//: complex.h
#ifndef _COMPLEX_H_
#define _COMPLEX_H_
namespace numeric {
//////////////////////////////////////////////////////////////////////////
#include <math.h>
#define LOG10 2.30258512496948
//////////////////////////////////////////////////////////////////////////
//template<typename _tp>
//class complex;
//
//template<class _arg1, class _arg2, class _oper>
//struct _binary_expr_cc;
//////////////////////////////////////////////////////////////////////////
template<typename _tp1, typename _tp2>
struct _return_type
{
typedef double value_type;
};
template<typename _tp1>
struct _return_type<_tp1, long double>
{
typedef long double value_type;
};
template<typename _tp2>
struct _return_type<long double, _tp2>
{
typedef long double value_type;
};
template<>
struct _return_type<float, float>
{
typedef float value_type;
};
//////////////////////////////////////////////////////////////////////////
template<typename _tp>
struct constnt
{
typedef typename _tp value_type;
public:
const _tp& _m_cnst;
public:
constnt(const _tp& __t) : _m_cnst(__t) {}
//constnt(const constnt& __c) : _m_cnst(__c._m_cnst) { }
public:
inline value_type real() const
{ return _m_cnst; }
inline value_type imag() const
{ return 0; }
};
//////////////////////////////////////////////////////////////////////////
struct _abs_cc
{
template<class _arg>
static inline typename _arg::value_type
evaluate(const _arg& __a)
{ return sqrt(__a.real()*__a.real()+__a.imag()*__a.imag()); }
};
struct _arg_cc
{
template<class _arg>
static inline typename _arg::value_type
evaluate(const _arg& __a)
{ return atan2(__a.imag(), __a.real()); }
};
struct _norm_cc
{
template<class _arg>
static inline typename _arg::value_type
evaluate(const _arg& __a)
{ return (__a.real()*__a.real()+__a.imag()*__a.imag()); }
};
//////////////////////////////////////////////////////////////////////////
struct _conjg_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{ return __a.real(); }
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{ return -__a.imag(); }
};
struct _cos_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{ return cos(__a.real())*cosh(__a.imag()); }
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{ return -sin(__a.real())*sinh(__a.imag()); }
};
struct _cosh_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{ return cosh(__a.real())*cos(__a.imag()); }
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{ return sinh(__a.real())*sin(__a.imag()); }
};
struct _exp_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{
typedef typename _arg::value_type value_type;
//value_type tmp = exp(__a.real());
return exp(__a.real())*cos(__a.imag());
}
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{
typedef typename _arg::value_type value_type;
//value_type tmp = exp(__a.real());
return exp(__a.real())*sin(__a.imag());
}
};
struct _log_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{
typedef typename _arg::value_type value_type;
value_type tmp = __a.real()*__a.real()+__a.imag()*__a.imag();
return log(tmp)*0.5;
}
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{
return atan2(__a.imag(), __a.real());
}
};
struct _log10_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{
typedef typename _arg::value_type value_type;
value_type tmp = __a.real()*__a.real()+__a.imag()*__a.imag();
return log10(tmp)*0.5;
}
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{
return arg(__a)/LOG10;
}
};
struct _sin_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{ return sin(__a.real())*cosh(__a.imag()); }
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{ return cos(__a.real())*sinh(__a.imag()); }
};
struct _sinh_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{ return sinh(__a.real())*cos(__a.imag()); }
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{ return cosh(__a.real())*sin(__a.imag()); }
};
//////////////////////////////////////////////////////////////////////////
struct _pow_cc
{
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
real(const _arg1& __a, const _arg2& __b)
{
//typedef typename
//_return_type<typename _arg1::value_type,
// typename _arg2::value_type
// >::value_type value_type;
//value_type t1 = log(__a.real()*__a.real()+__a.imag()*__a.imag())*0.5;
//value_type t2 = arg(__a);
//return exp(t1*__b.real()-t2*__b.imag())*cos(t1*__b.imag()+t2*__b.real());
return exp(__b*log(__a)).real();
}
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
imag(const _arg1& __a, const _arg2& __b)
{
//typedef typename
//_return_type<typename _arg1::value_type,
// typename _arg2::value_type
// >::value_type value_type;
//value_type t1 = log(__a.real()*__a.real()+__a.imag()*__a.imag())*0.5;
//value_type t2 = arg(__a);
//return exp(t1*__b.real()-t2*__b.imag())*sin(t1*__b.imag()+t2*__b.real());
return exp(__b*log(__a)).imag();
}
};
//////////////////////////////////////////////////////////////////////////
struct _unary_add_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{ return __a.real(); }
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{ return __a.imag(); }
};
//////////////////////////////////////////////////////////////////////////
struct _unary_sub_cc
{
template<class _arg>
static inline typename _arg::value_type
real(const _arg& __a)
{ return -__a.real(); }
template<class _arg>
static inline typename _arg::value_type
imag(const _arg& __a)
{ return -__a.imag(); }
};
//////////////////////////////////////////////////////////////////////////
struct _binary_add_cc
{
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
real(const _arg1& __a, const _arg2& __b)
{ return __a.real() + __b.real(); }
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
imag(const _arg1& __a, const _arg2& __b)
{ return __a.imag() + __b.imag(); }
};
//////////////////////////////////////////////////////////////////////////
struct _binary_sub_cc
{
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
real(const _arg1& __a, const _arg2& __b)
{ return __a.real() - __b.real(); }
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
imag(const _arg1& __a, const _arg2& __b)
{ return __a.imag() - __b.imag(); }
};
//////////////////////////////////////////////////////////////////////////
struct _binary_mul_cc
{
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
real(const _arg1& __a, const _arg2& __b)
{ return __a.real()*__b.real() - __a.imag()*__b.imag(); }
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
imag(const _arg1& __a, const _arg2& __b)
{ return __a.imag()*__b.real() + __a.real()*__b.imag(); }
};
//////////////////////////////////////////////////////////////////////////
struct _binary_div_cc
{
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
real(const _arg1& __a, const _arg2& __b)
{ return (__a.real()*__b.real() + __a.imag()*__b.imag())/
(__b.real()*__b.real() + __b.imag()*__b.imag()); }
template<class _arg1, class _arg2>
static inline typename _return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
imag(const _arg1& __a, const _arg2& __b)
{ return (__a.imag()*__b.real() - __a.real()*__b.imag())/
(__b.real()*__b.real() + __b.imag()*__b.imag()); }
};
//////////////////////////////////////////////////////////////////////////
template<class _arg, class _oper>
struct _unary_expr_cc
{
typedef typename _arg::value_type value_type;
public:
const _arg& _m_arg;
public:
_unary_expr_cc(const _arg& __a) : _m_arg(__a)
{ }
inline value_type real() const
{ return _oper::real(_m_arg); }
inline value_type imag() const
{ return _oper::imag(_m_arg); }
public:
_unary_expr_cc<_unary_expr_cc<_arg,_oper>,_unary_add_cc>
operator+() const
{
typedef _unary_expr_cc<_unary_expr_cc<_arg,_oper>,
_unary_add_cc> _expr;
return _expr(*this);
}
_unary_expr_cc<_unary_expr_cc<_arg,_oper>,_unary_sub_cc>
operator-() const
{
typedef _unary_expr_cc<_unary_expr_cc<_arg,_oper>,
_unary_sub_cc> _expr;
return _expr(*this);
}
};
//////////////////////////////////////////////////////////////////////////
template<class _arg1, class _arg2, class _oper>
struct _binary_expr_cc
{
typedef typename
_return_type<typename _arg1::value_type,
typename _arg2::value_type>::value_type
value_type;
public:
const _arg1& _m_arg1;
const _arg2& _m_arg2;
public:
_binary_expr_cc(const _arg1& __a, const _arg2& __b)
: _m_arg1(__a), _m_arg2(__b)
{ }
public:
inline value_type real() const
{ return _oper::real(_m_arg1, _m_arg2); }
inline value_type imag() const
{ return _oper::imag(_m_arg1, _m_arg2); }
public:
_unary_expr_cc<_binary_expr_cc<_arg1,_arg2,_oper>,_unary_add_cc>
operator+() const
{
typedef _unary_expr_cc<_binary_expr_cc<_arg1,_arg2,_oper>,
_unary_add_cc
> _expr;
return _expr(*this);
}
_unary_expr_cc<_binary_expr_cc<_arg1,_arg2,_oper>,_unary_sub_cc>
operator-() const
{
typedef _unary_expr_cc<_binary_expr_cc<_arg1,_arg2,_oper>,
_unary_sub_cc
> _expr;
return _expr(*this);
}
};
//////////////////////////////////////////////////////////////////////////
#ifdef _DEBUG
// 在debug下，由于编译器对生成的临时变量严格按照C++标准来处理(也就是临时
// 变量生存期问题)，所以需要针对常量类型进行特化。但是现在的编译器在
// release下都能对代码作出很好的优化，临时变量被优化掉了，所以没有特化
// 也能得到正确的结果!
template<typename _tp, class _arg2, class _oper>
struct _binary_expr_cc<constnt<_tp>,_arg2,_oper>
{
typedef typename
_return_type<typename _tp,
typename _arg2::value_type>::value_type
value_type;
public:
const constnt<_tp> _m_arg1;
const _arg2& _m_arg2;
public:
_binary_expr_cc(const constnt<_tp>& __a, const _arg2& __b)
: _m_arg1(__a), _m_arg2(__b)
{ //cout << "hello" << endl; }
public:
inline value_type real() const
{ return _oper::real(_m_arg1, _m_arg2); }
inline value_type imag() const
{ return _oper::imag(_m_arg1, _m_arg2); }
public:
_unary_expr_cc<_binary_expr_cc<constnt<_tp>,_arg2,_oper>,_unary_add_cc>
operator+() const
{
typedef _unary_expr_cc<_binary_expr_cc<constnt<_tp>,_arg2,_oper>,
_unary_add_cc
> _expr;
return _expr(*this);
}
};
//////////////////////////////////////////////////////////////////////////
template<class _arg1, class _tp, class _oper>
struct _binary_expr_cc<_arg1,constnt<_tp>,_oper>
{
typedef typename
_return_type<typename _arg1::value_type,
typename _tp
>::value_type value_type;
public:
const _arg1& _m_arg1;
const constnt<_tp> _m_arg2;
public:
_binary_expr_cc(const _arg1& __a, const constnt<_tp>& __b)
: _m_arg1(__a), _m_arg2(__b)
{ }
public:
inline value_type real() const
{ return _oper::real(_m_arg1, _m_arg2); }
inline value_type imag() const
{ return _oper::imag(_m_arg1, _m_arg2); }
public:
_unary_expr_cc<_binary_expr_cc<_arg1,constnt<_tp>,_oper>,
_unary_add_cc
>
operator+() const
{
typedef _unary_expr_cc<_binary_expr_cc<_arg1,constnt<_tp>,_oper>,
_unary_add_cc
> _expr;
return _expr(*this);
}
_unary_expr_cc<_binary_expr_cc<_arg1,constnt<_tp>,_oper>,
_unary_sub_cc
>
operator-() const
{
typedef _unary_expr_cc<_binary_expr_cc<_arg1,constnt<_tp>,_oper>,
_unary_sub_cc
> _expr;
return _expr(*this);
}
};
#endif // _DEBUG
//////////////////////////////////////////////////////////////////////////
template<typename _tp>
class complex
{
public:
typedef typename _tp value_type;
protected:
value_type _m_real;
value_type _m_imag;
public:
complex(const _tp& __real = _tp(), const _tp& __imag = _tp())
: _m_real(__real), _m_imag(__imag)
{
#ifdef _DEBUG
if ((typeid(value_type) != typeid(float)) &&
(typeid(value_type) != typeid(double)) &&
(typeid(value_type) != typeid(long double))) {
cout << "complex模板参数类型只支持浮点数!" << endl;
exit(0);
}
#endif // _DEBUG
}
complex(const complex& __c)
: _m_real(__c._m_real), _m_imag(__c._m_imag)
{ }
public:
template<class _arg, class _oper>
complex(const _unary_expr_cc<_arg,_oper>& __e)
: _m_real(__e.real()), _m_imag(__e.imag())
{ }
template<class _arg1, class _arg2, class _oper>
complex(const _binary_expr_cc<_arg1,_arg2,_oper>& __e)
: _m_real(__e.real()), _m_imag(__e.imag())
{ }
public:
inline complex& operator=(const complex& __c)
{
_m_real = __c.real();
_m_imag = __c.imag();
return *this;
}
template<class _arg, class _oper>
inline complex& operator=(const _unary_expr_cc<_arg,_oper>& __e)
{
_m_real = __e.real();
_m_imag = __e.imag();
return *this;
}
template<class _arg1, class _arg2, class _oper>
inline complex& operator=(const _binary_expr_cc<_arg1,_arg2,_oper>& __e)
{
_m_real = __e.real();
_m_imag = __e.imag();
return *this;
}
public:
inline complex& operator+=(const complex& __c)
{
_m_real += __c.real();
_m_imag += __c.imag();
return *this;
}
template<class _arg, class _oper>
inline complex& operator+=(const _unary_expr_cc<_arg,_oper>& __e)
{
_m_real += __e.real();
_m_imag += __e.imag();
return *this;
}
template<class _arg1, class _arg2, class _oper>
inline complex& operator+=(const _binary_expr_cc<_arg1,_arg2,_oper>& __e)
{
_m_real += __e.real();
_m_imag += __e.imag();
return *this;
}
inline complex& operator-=(const complex& __c)
{
_m_real -= __c.real();
_m_imag -= __c.imag();
return *this;
}
template<class _arg, class _oper>
inline complex& operator-=(const _unary_expr_cc<_arg,_oper>& __e)
{
_m_real -= __e.real();
_m_imag -= __e.imag();
return *this;
}
template<class _arg1, class _arg2, class _oper>
inline complex& operator-=(const _binary_expr_cc<_arg1,_arg2,_oper>& __e)
{
_m_real -= __e.real();
_m_imag -= __e.imag();
return *this;
}
inline complex& operator*=(const complex& __c)
{
static _tp rx, iy;
rx = _m_real*__c.real() - _m_imag*__c.imag();
iy = _m_real*__c.imag() + _m_imag*__c.real();
_m_real = rx;
_m_imag = iy;
return *this;
}
template<class _arg, class _oper>
inline complex& operator*=(const _unary_expr_cc<_arg,_oper>& __e)
{
static _tp rx, iy, ex, ey;
ex = __e.real();
ey = __e.imag();
rx = _m_real*ex - _m_imag*ey;
iy = _m_real*ey + _m_imag*ex;
_m_real = rx;
_m_imag = iy;
return *this;
}
template<class _arg1, class _arg2, class _oper>
inline complex& operator*=(const _binary_expr_cc<_arg1,_arg2,_oper>& __e)
{
static _tp rx, iy;
rx = _m_real*__e.real() - _m_imag*__e.imag();
iy = _m_real*__e.imag() + _m_imag*__e.real();
_m_real = rx;
_m_imag = iy;
return *this;
}
inline complex& operator/=(const complex& __c)
{
static _tp rx, iy, cx, cy, rt;
cx = __c.real();
cy = __c.imag();
rt = 1/(cx*cx+cy*cy);
rx = _m_real*cx + _m_imag*cy;
iy = -_m_real*cy + _m_imag*cx;
_m_real = rx*rt;
_m_imag = iy*rt;
return *this;
}
template<class _arg, class _oper>
inline complex& operator/=(const _unary_expr_cc<_arg,_oper>& __e)
{
static _tp rx, iy, ex, ey, rt;
ex = __e.real();
ey = __e.imag();
rt = 1/(ex*ex+ey*ey);
rx = _m_real*ex + _m_imag*ey;
iy = -_m_real*ey + _m_imag*ex;
_m_real = rx*rt;
_m_imag = iy*rt;
return *this;
}
template<class _arg1, class _arg2, class _oper>
inline complex& operator/=(const _binary_expr_cc<_arg1,_arg2,_oper>& __e)
{
static _tp rx, iy, ex, ey, rt;
ex = __e.real();
ey = __e.imag();
rt = 1/(ex*ex+ey*ey);
rx = _m_real*ex + _m_imag*ey;
iy = -_m_real*ey + _m_imag*ex;
_m_real = rx*rt;
_m_imag = iy*rt;
return *this;
}
public:
inline value_type real() const { return _m_real; }
inline value_type imag() const { return _m_imag; }
inline value_type& real() { return _m_real; }
inline value_type& imag() { return _m_imag; }
public:
_unary_expr_cc<complex<_tp>,_unary_add_cc> operator+() const
{ return _unary_expr_cc<complex<_tp>,_unary_add_cc>(*this); }
_unary_expr_cc<complex<_tp>,_unary_sub_cc> operator-() const
{ return _unary_expr_cc<complex<_tp>,_unary_sub_cc>(*this); }
};
//////////////////////////////////////////////////////////////////////////
//template<class _arg1, class _arg2>
//_binary_expr_cc<_arg1,_arg2,_add_cc>
//operator+(const _arg1& a, const _arg2& b)
//{
// typedef _binary_expr_cc<_arg1,_arg2,_add_cc> _expr;
// return _expr(a,b);
//}
#define _DEFINE_COMPLEX_BINARY_FUNC(_func, _name) template<class _arg11,class _arg12,class _oper1, class _arg21,class _arg22,class _oper2> inline _binary_expr_cc<_binary_expr_cc<_arg11,_arg12,_oper1>, _binary_expr_cc<_arg21,_arg22,_oper2>, _name > _func(const _binary_expr_cc<_arg11,_arg12,_oper1>& a, const _binary_expr_cc<_arg21,_arg22,_oper2>& b) { typedef _binary_expr_cc<_binary_expr_cc<_arg11,_arg12,_oper1>, _binary_expr_cc<_arg21,_arg22,_oper2>, _name > _expr; return _expr(a,b); } template<class _arg1,class _oper1, class _arg21,class _arg22,class _oper2> inline _binary_expr_cc<_unary_expr_cc<_arg1,_oper1>, _binary_expr_cc<_arg21,_arg22,_oper2>, _name > _func(const _unary_expr_cc<_arg1,_oper1>& a, const _binary_expr_cc<_arg21,_arg22,_oper2>& b) { typedef _binary_expr_cc<_unary_expr_cc<_arg1,_oper1>, _binary_expr_cc<_arg21,_arg22,_oper2>, _name > _expr; return _expr(a,b); } template<class _arg11,class _arg12,class _oper1, class _arg2,class _oper2> inline _binary_expr_cc<_binary_expr_cc<_arg11,_arg12,_oper1>, _unary_expr_cc<_arg2,_oper2>, _name > _func(const _binary_expr_cc<_arg11,_arg12,_oper1>& a, const _unary_expr_cc<_arg2,_oper2>& b) { typedef _binary_expr_cc<_binary_expr_cc<_arg11,_arg12,_oper1>, _unary_expr_cc<_arg2,_oper2>, _name > _expr; return _expr(a,b); } template<class _arg1,class _oper1, class _arg2,class _oper2> inline _binary_expr_cc<_unary_expr_cc<_arg1,_oper1>, _unary_expr_cc<_arg2,_oper2>, _name > _func(const _unary_expr_cc<_arg1,_oper1>& a, const _unary_expr_cc<_arg2,_oper2>& b) { typedef _binary_expr_cc<_unary_expr_cc<_arg1,_oper1>, _unary_expr_cc<_arg2,_oper2>, _name > _expr; return _expr(a,b); } template<class _arg11, class _arg12, class _oper, class _tp> inline _binary_expr_cc<_binary_expr_cc<_arg11,_arg12,_oper>, complex<_tp>, _name > _func(const _binary_expr_cc<_arg11,_arg12,_oper>& a, const complex<_tp>& b) { typedef _binary_expr_cc<_binary_expr_cc<_arg11,_arg12,_oper>, complex<_tp>, _name > _expr; return _expr(a,b); } template<class _tp, class _arg21, class _arg22, class _oper> inline _binary_expr_cc<complex<_tp>, _binary_expr_cc<_arg21,_arg22,_oper>, _name > _func(const complex<_tp>& a, const _binary_expr_cc<_arg21,_arg22,_oper>& b) { typedef _binary_expr_cc<complex<_tp>, _binary_expr_cc<_arg21,_arg22,_oper>, _name > _expr; return _expr(a,b); } template<class _arg11, class _arg12, class _oper, class _tp> inline _binary_expr_cc<_binary_expr_cc<_arg11,_arg12,_oper>, constnt<_tp>, _name > _func(const _binary_expr_cc<_arg11,_arg12,_oper>& a, const _tp& b) { typedef _binary_expr_cc<_binary_expr_cc<_arg11,_arg12,_oper>, constnt<_tp>, _name > _expr; return _expr(a,b); } template<class _tp, class _arg21, class _arg22, class _oper> inline _binary_expr_cc<constnt<_tp>, _binary_expr_cc<_arg21,_arg22,_oper>, _name > _func(const _tp& a, const _binary_expr_cc<_arg21,_arg22,_oper>& b) { typedef _binary_expr_cc<constnt<_tp>, _binary_expr_cc<_arg21,_arg22,_oper>, _name > _expr; return _expr(a,b); } template<class _arg,class _oper,class _tp> inline _binary_expr_cc<_unary_expr_cc<_arg,_oper>, complex<_tp>, _name > _func(const _unary_expr_cc<_arg,_oper>& a, const complex<_tp>& b) { typedef _binary_expr_cc<_unary_expr_cc<_arg,_oper>, complex<_tp>, _name > _expr; return _expr(a,b); } template<class _tp, class _arg,class _oper> inline _binary_expr_cc<complex<_tp>, _unary_expr_cc<_arg,_oper>, _name > _func(const complex<_tp>& a, const _unary_expr_cc<_arg,_oper>& b) { typedef _binary_expr_cc<complex<_tp>, _unary_expr_cc<_arg,_oper>, _name > _expr; return _expr(a,b); } template<class _arg,class _oper,class _tp> inline _binary_expr_cc<_unary_expr_cc<_arg,_oper>, constnt<_tp>, _name > _func(const _unary_expr_cc<_arg,_oper>& a, const constnt<_tp>& b) { /*这个函数不常用，只是用作备用*/ typedef _binary_expr_cc<_unary_expr_cc<_arg,_oper>, constnt<_tp>, _name > _expr; return _expr(a,b); } template<class _tp, class _arg,class _oper> inline _binary_expr_cc<constnt<_tp>, _unary_expr_cc<_arg,_oper>, _name > _func(const constnt<_tp>& a, const _unary_expr_cc<_arg,_oper>& b) { typedef _binary_expr_cc<constnt<_tp>, _unary_expr_cc<_arg,_oper>, _name > _expr; return _expr(a,b); } template<class _arg,class _oper,class _tp> inline _binary_expr_cc<_unary_expr_cc<_arg,_oper>, constnt<_tp>, _name > _func(const _unary_expr_cc<_arg,_oper>& a, const _tp& b) { typedef _binary_expr_cc<_unary_expr_cc<_arg,_oper>, constnt<_tp>, _name > _expr; return _expr(a,b); } template<class _tp, class _arg,class _oper> inline _binary_expr_cc<constnt<_tp>, _unary_expr_cc<_arg,_oper>, _name > _func(const _tp& a, const _unary_expr_cc<_arg,_oper>& b) { typedef _binary_expr_cc<constnt<_tp>, _unary_expr_cc<_arg,_oper>, _name > _expr; return _expr(a,b); } template<class _tp1, class _tp2> inline _binary_expr_cc<complex<_tp1>,complex<_tp2>,_name> _func(const complex<_tp1>& __a, const complex<_tp2>& __b) { typedef _binary_expr_cc<complex<_tp1>, complex<_tp2>, _name > _expr; return _expr(__a,__b); } template<class _tp1, class _tp2> inline _binary_expr_cc<complex<_tp1>, constnt<_tp2>, _name > _func(const complex<_tp1>& __a, const _tp2& __b) { typedef _binary_expr_cc<complex<_tp1>, constnt<_tp2>, _name > _expr; return _expr(__a,__b); } template<class _tp1, class _tp2> inline _binary_expr_cc<constnt<_tp1>, complex<_tp2>, _name > _func(const _tp1& a, const complex<_tp2>& b) { typedef _binary_expr_cc<constnt<_tp1>, complex<_tp2>, _name > _expr; return _expr(a,b); }
_DEFINE_COMPLEX_BINARY_FUNC(operator+, _binary_add_cc)
_DEFINE_COMPLEX_BINARY_FUNC(operator-, _binary_sub_cc)
_DEFINE_COMPLEX_BINARY_FUNC(operator*, _binary_mul_cc)
_DEFINE_COMPLEX_BINARY_FUNC(operator/, _binary_div_cc)
_DEFINE_COMPLEX_BINARY_FUNC( pow, _pow_cc)
//////////////////////////////////////////////////////////////////////////
#define _DEFINE_COMPLEX_UNARY_FUNC(_func, _name) template<class _arg1, class _arg2, class _oper> _unary_expr_cc<_binary_expr_cc<_arg1,_arg2,_oper>, _name > _func(const _binary_expr_cc<_arg1,_arg2,_oper>& __e) { typedef _unary_expr_cc<_binary_expr_cc<_arg1,_arg2,_oper>, _name > _expr; return _expr(__e); } template<class _arg, class _oper> _unary_expr_cc<_unary_expr_cc<_arg,_oper>, _name > _func(const _unary_expr_cc<_arg,_oper>& __e) { typedef _unary_expr_cc<_unary_expr_cc<_arg,_oper>, _name > _expr; return _expr(__e); } template<class _tp> _unary_expr_cc<complex<_tp>, _name > _func(const complex<_tp>& __c) { typedef _unary_expr_cc<complex<_tp>, _name > _expr; return _expr(__c); }
_DEFINE_COMPLEX_UNARY_FUNC( cos, _cos_cc)
_DEFINE_COMPLEX_UNARY_FUNC( cosh, _cosh_cc)
_DEFINE_COMPLEX_UNARY_FUNC( exp, _exp_cc)
_DEFINE_COMPLEX_UNARY_FUNC( log, _log_cc)
_DEFINE_COMPLEX_UNARY_FUNC(log10, _log10_cc)
_DEFINE_COMPLEX_UNARY_FUNC( sin, _sin_cc)
_DEFINE_COMPLEX_UNARY_FUNC( sinh, _sinh_cc)
//////////////////////////////////////////////////////////////////////////
#define _DEFINE_COMPLEX_UNARY_FUNC_EX(_func, _name) template<class _arg1, class _arg2, class _oper> typename _return_type<typename _arg1::value_type, typename _arg2::value_type >::value_type _func(const _binary_expr_cc<_arg1,_arg2,_oper>& __e) { return _name::evaluate(__e); } template<class _arg, class _oper> typename _arg::value_type _func(const _unary_expr_cc<_arg,_oper>& __e) { return _name::evaluate(__e); } template<typename _tp> typename _tp _func(const complex<_tp>& __c) { return _name::evaluate(__c); }
_DEFINE_COMPLEX_UNARY_FUNC_EX( abs, _abs_cc)
_DEFINE_COMPLEX_UNARY_FUNC_EX( arg, _arg_cc)
_DEFINE_COMPLEX_UNARY_FUNC_EX(norm, _norm_cc)
//////////////////////////////////////////////////////////////////////////
}
#endif // _COMPLEX_H_
// test.cpp
#include <ctime>
#include <iostream>
using namespace std;
//#include <complex>
//using namespace std;
#include "complex.h"
using namespace numeric;
void main()
{
complex<double> c1(1,1), c2(2,2), c3(3,3), c4(4,4), c5(5,5), cc(0,0);
double d = 1.2, rx = 0, iy = 0;
clock_t time = clock();
for (int i=0; i<10000; ++i) {
for (int j=0; j<10000; ++j) {
// 由于被编译器的循环优化会优化=,所以采用@=
cc *= cos(d+(c1+c2)+c3);
}
}
time = clock() - time;
cout << (double)time/CLK_TCK << endl;
cout << '(' << rx
<< ',' << iy
<< ')' << endl;
cout << '(' << cc.real()
<< ',' << cc.imag()
<< ')' << endl;
}
/*************************************************************************
程序分别用intel c++ compiler 8.0, vc7.1和gcc 3.3.3编译。其中icc80和vc71
顺利编译通过，gcc编译出错。
首先用intel的编译器来测试程序：
(1) 简单运算(+,-,*,/,=,+=,-=,*=,/=)
通过测试发现，它在+、-及+=、-=上的效率约是compaq visual fortran 6.5
的complex类型的两倍多。而在稍微复杂的*、/、*=、/=上，它的效率要高出
cvf好多。
(2) 函数计算(cos, sin, exp, log)
这些函数计算与cvf相比仍然约是两倍多(fortran没有log10、sinh和cosh)。
(3) pow函数
在测试的时候发现一个有趣的问题：intel编译器在处理pow函数时对数据非常
敏感。比如申明
complex<double> c1(1,1),c2(2,2),c3(3,3),cc(0,0);
double d = 1.2;
然后计算下面语句10000*10000次
cc += pow(d+c1+c2,c3);
如果将c2(2,2)的虚部改为任何不为2的数，对比两次计算耗时，可以发现两者
差别居然有近70倍。
fortan没有pow函数，但是pow(x,y)可以表达成exp(y*log(x))。
使用vs71编译器来测试程序：
(1) 简单运算(+,-,*,/,=,+=,-=,*=,/=)
编译出的程序运行效率和visual fortran差不多，+=和-=稍微低于fortran，而
*=则约是fortran的1/3，但是/=居然有fortran的两倍。
(2) 函数计算(cos, sin, exp, sinh, cosh)
使用@=(@可以是+、-、*、/)运算符时，计算的效率和stl的complex没什么差别。
(3) 函数(log, log10)
这两个函数与(2)中的函数表现差不多，但是如果将@=(+=或-=)中的一个赋值
表达式注释掉，效率就会高出很多，但这个方法对(2)中的函数不管用。
(4) 函数(pow)
这个函数如果直接使用@=的话，也和stl的complex类差不多。但是(3)中的方法对
它也有效，更有意思的地方是——如果在循环中直接赋值：
cc.real() += pow(d+c1+c2,c3).real();
cc.imag() += pow(d+c1+c2,c3).imag();
这样得到的效率居然是fortran的两倍多。但是这种方法对(3)(2)都不管用。
*************************************************************************/
