c++复数类库.txt

将数学运算中的复数运算抽象为C++类的源代码

#define G_COMPLEX_CTOR_OP_IMP(type,intype)     \
   G_COMPLEX(type)::G_COMPLEX(type)(const G_COMPLEX(intype) &v) : re((type)v.re) { im = (type)v.im; }

#define G_COMPLEX_OPEQ_DECL(type,from)      \
  G_COMPLEX(type)& operator += (from y) { re += (type) y; return *this; } \
  G_COMPLEX(type)& operator += (const G_COMPLEX(from) &y);   \
  G_COMPLEX(type)& operator -= (from y) { re -= (type) y; return *this; } \
  G_COMPLEX(type)& operator -= (const G_COMPLEX(from) &y);   \
  G_COMPLEX(type)& operator *= (from y) { re *=  (type) y; im *= (type) y; return *this;}\
  G_COMPLEX(type)& operator *= (const G_COMPLEX(from) &y);   \
  G_COMPLEX(type)& operator /= (from y);     \
  G_COMPLEX(type)& operator /= (const G_COMPLEX(from) &y);

#define G_COMPLEX_OPEQ_IMP(type,from)      \
G_COMPLEX(type)& G_COMPLEX(type)::operator += (const G_COMPLEX(from) &y) \
  { re += (type) y.re; im += (type) y.im; return *this; } \
          \
G_COMPLEX(type)& G_COMPLEX(type)::operator -= (const G_COMPLEX(from) &y) \
  { re -= (type) y.re; im -= (type) y.im; return *this; } \
          \
G_COMPLEX(type)& G_COMPLEX(type)::operator *= (const G_COMPLEX(from) &y) { \
   type r = re * (type) y.re - im * (type) y.im; \
   im = re * (type) y.im + im * (type) y.re;  \
   re = r;       \
   return *this;       \
        }        \
          \
G_COMPLEX(type)& G_COMPLEX(type)::operator /= (const G_COMPLEX(from)& y) \
{          \
  double den = fabs((double)y.re) + fabs((double)y.im);   \
  if (den == 0.0) error ("Attempted division by zero.");   \
  double xrden = (double) re / den;      \
  double xiden = (double) im / den;      \
  double yrden = (double) y.re / den;     \
  double yiden = (double) y.im / den;     \
  double nrm   = yrden * yrden + yiden * yiden;     \
  re = (type)((xrden * yrden + xiden * yiden) / nrm);    \
  im = (type)((xiden * yrden - xrden * yiden) / nrm);    \
  return *this;         \
}          \
          \
G_COMPLEX(type)& G_COMPLEX(type)::operator /= (from y)    \
{          \
  if (y == (from) 0) error ("Attempted division by zero.");   \
  re /= (type) y;  im /= (type) y;      \
  return *this;         \
}



class G_COMPLEX(double);
class G_COMPLEX(float);
class G_COMPLEX(int);

g_declare2(G_COMPLEX,double,G_COMPLEX_CTOR_OP_DECL(double,int) G_COMPLEX_CTOR_OP_DECL(double,float) G_COMPLEX_ASSIGN_OP_DECLS(double,double) G_COMPLEX_ASSIGN_OP_DECL(double,float) G_COMPLEX_ASSIGN_OP_DECL(double,int) G_COMPLEX_OPEQ_DECL(double,float) G_COMPLEX_OPEQ_DECL(double,int))
typedef G_COMPLEX(double) DComplex;
g_declare2(G_COMPLEX,float,G_COMPLEX_CAST_OP(double) G_COMPLEX_CTOR_OP_DECL(float,int) G_COMPLEX_CTOR_OP_DECL(float,double) G_COMPLEX_ASSIGN_OP_DECL(float,double) G_COMPLEX_ASSIGN_OP_DECL(float,int) G_COMPLEX_OPEQ_DECL(float,int) G_COMPLEX_OPEQ_DECL(float,double))
typedef G_COMPLEX(float) Complex;
g_declare2(G_COMPLEX,int,G_COMPLEX_CAST_OP(float) G_COMPLEX_CAST_OP(double) G_COMPLEX_CTOR_OP_DECL(int,float) G_COMPLEX_CTOR_OP_DECL(int,double) G_COMPLEX_ASSIGN_OP_DECL(int,double) G_COMPLEX_ASSIGN_OP_DECL(int,float) G_COMPLEX_OPEQ_DECL(int,float) G_COMPLEX_OPEQ_DECL(int,double))
typedef G_COMPLEX(int) IComplex;

typedef Complex float_complex;
typedef DComplex double_complex;

/*
* POLAR
*/
inline G_COMPLEX(double)  polar(double r, double t)
{
  return G_COMPLEX(double)(r * cos(t), r * sin(t));
}

/*
* near() functions
*/

Bool near(G_COMPLEX(float) val1, G_COMPLEX(float) val2, double tol=1.0e-5);
Bool near(G_COMPLEX(double) val1, G_COMPLEX(double) val2, double tol=1.0e-13);
Bool nearAbs(G_COMPLEX(float) val1, G_COMPLEX(float) val2, double tol=1.0e-5);
Bool nearAbs(G_COMPLEX(double) val1, G_COMPLEX(double) val2, double tol=1.0e-13);

inline Bool allNear(G_COMPLEX(float) val1, G_COMPLEX(float) val2, 
      double tol=1.0e-5)
    { return near(val1, val2, tol); }
inline Bool allNear(G_COMPLEX(double) val1, G_COMPLEX(double) val2, 
      double tol=1.0e-13)
    { return near(val1, val2, tol); }
inline Bool allNearAbs(G_COMPLEX(float) val1, G_COMPLEX(float) val2, 
         double tol=1.0e-5)
    { return nearAbs(val1, val2, tol); }
inline Bool allNearAbs(G_COMPLEX(double) val1, G_COMPLEX(double) val2, 
         double tol=1.0e-13)
    { return nearAbs(val1, val2, tol); }

#if defined(COMPLEX_STUPID_COMPILER)
//#***************************************************************************
//#********************* Should not be needed, casts *************************
//#********************* should take care of these   *************************
//#***************************************************************************
#define G_COMPLEX_BIN_OP_DECL(ret,left,right)     \
/*          \
* non-inline functions        \
*/          \
G_COMPLEX(ret) operator /  (const G_COMPLEX(left)& x, const G_COMPLEX(right)& y);\
G_COMPLEX(ret) operator /  (const G_COMPLEX(left)& x, right y);   \
G_COMPLEX(ret) operator /  (left x, const G_COMPLEX(right)& y);   \
          \
/*          \
* EQUALITY OPERATORS        \
*/          \
inline int  operator == (const G_COMPLEX(left)& x, const G_COMPLEX(right)& y) \
{          \
  return (ret) x.re == (ret) y.re && (ret) x.im == (ret) y.im; \
}          \
          \
inline int  operator == (const G_COMPLEX(left)& x, right y)   \
{          \
  return (ret) x.im == (ret) 0.0 && (ret) x.re == (ret) y;  \
}          \
          \
inline int  operator != (const G_COMPLEX(left)& x, const G_COMPLEX(right)& y) \
{          \
  return (ret) x.re != (ret) y.re || (ret) x.im != (ret) y.im; \
}          \
          \
inline int  operator != (const G_COMPLEX(left)& x, right y)   \
{          \
  return (ret) x.im != (ret) 0.0 || (ret) x.re != (ret) y;  \
}          \
          \
/*          \
* COMPARISON OPERATORS (using norm)      \
*/          \
inline int operator >= (const G_COMPLEX(left) &x, const G_COMPLEX(right) &y) { \
  return norm(x) >= norm(y);       \
}          \
          \
inline int operator >= (const G_COMPLEX(left) &x, const right y) {  \
  return norm(x) >= (double) y * (double) y;     \
}          \
          \
inline int operator >= (const left x, const G_COMPLEX(right) &y) {  \
  return (double) x * (double) x >= norm(y);     \
}          \
          \
inline int operator > (const G_COMPLEX(left) &x, const G_COMPLEX(right) &y) { \
  return norm(x) > norm(y);       \
}          \
          \
inline int operator > (const G_COMPLEX(left) &x, const right y) {  \
  return norm(x) > (double) y * (double) y;     \
}          \
          \
inline int operator > (const left x, const G_COMPLEX(right) &y) {  \
  return (double) x * (double) x > norm(y);     \
}          \
          \
inline int operator <= (const G_COMPLEX(left) &x, const G_COMPLEX(right) &y) { \
  return norm(x) <= norm(y);       \
}          \
          \
inline int operator <= (const G_COMPLEX(left) &x, const right y) {  \
  return norm(x) <= (double) y * (double)y;     \
}          \
          \
inline int operator <= (const left x, const G_COMPLEX(right) &y) {  \
  return (double) x * (double) x <= norm(y);     \
}          \
          \
inline int operator < (const G_COMPLEX(left) &x, const G_COMPLEX(right) &y) { \
  return norm(x) < norm(y);       \
}          \
          \
inline int operator < (const G_COMPLEX(left) &x, const right y) {  \
  return norm(x) < (double) y * (double) y;     \
}          \
          \
inline int operator < (const left x, const G_COMPLEX(right) &y) {  \
  return (double) x * (double) x < norm(y);     \
}          \
          \
/*          \
* SIMPLE NUMERIC OPERATIONS       \
*/          \
inline G_COMPLEX(ret) operator + (const G_COMPLEX(left)& x, const G_COMPLEX(right)& y)\
{          \
  return G_COMPLEX(ret)((ret) x.re + (ret) y.re, (ret) x.im + (ret) y.im);\
}          \
          \
inline G_COMPLEX(ret) operator + (const G_COMPLEX(left)& x, right y)  \
{          \
  return G_COMPLEX(ret)((ret) x.re + (ret) y, (ret) x.im);  \
}          \
          \
inline G_COMPLEX(ret) operator + (left x, const G_COMPLEX(right)& y)  \
{          \
  return G_COMPLEX(ret)((ret) x + (ret) y.re, (ret) y.im);  \
}          \
          \
inline G_COMPLEX(ret) operator - (const G_COMPLEX(left)& x, const G_COMPLEX(right)& y)\
{          \
  return G_COMPLEX(ret)((ret) x.re - (ret) y.re, (ret) x.im - (ret) y.im);\
}          \
          \
inline G_COMPLEX(ret) operator - (const G_COMPLEX(left)& x, right y)  \
{          \
  return G_COMPLEX(ret)((ret) x.re - (ret) y, (ret) x.im);  \
}          \
          \
inline G_COMPLEX(ret) operator - (left x, const G_COMPLEX(right)& y)  \
{          \
  return G_COMPLEX(ret)((ret) x - (ret) y.re, - (ret) y.im);  \
}          \
          \
inline G_COMPLEX(ret) operator * (const G_COMPLEX(left)& x, const G_COMPLEX(right)& y)\
{          \
  return G_COMPLEX(ret)((ret) x.re * (ret) y.re - (ret) x.im * (ret) y.im,\
                 (ret) x.re * (ret) y.im + (ret) x.im * (ret) y.re);\
}          \
          \
inline G_COMPLEX(ret) operator * (const G_COMPLEX(left)& x, right y)  \
{          \
  return G_COMPLEX(ret)((ret) x.re * (ret) y, (ret) x.im * (ret) y); \
}          \
          \
inline G_COMPLEX(ret) operator * (left x, const G_COMPLEX(right)& y)  \
{          \
  return G_COMPLEX(ret)((ret) x * (ret) y.re, (ret) x * (ret) y.im); \
}          \
          \
G_COMPLEX(ret) pow(const G_COMPLEX(left) &x, const G_COMPLEX(right) &p);


#define G_COMPLEX_BIN_OP_IMP(ret,left,right)     \
          \
/* from [e-mail] */      \
G_COMPLEX(ret) /* const */ operator / (const G_COMPLEX(left)& x, const G_COMPLEX(right)& y)\
{          \
  double den = fabs((double)y.re) + fabs((double)y.im);   \
  if (den == 0.0) x.error ("Attempted division by zero.");   \
  double xrden = x.re / den;      \
  double xiden = x.im / den;      \
  double yrden = y.re / den;      \
  double yiden = y.im / den;      \
  double nrm   = yrden * yrden + yiden * yiden;     \
  return G_COMPLEX(ret)((ret)((xrden * yrden + xiden * yiden) / nrm),  \
                 (ret)((xiden * yrden - xrden * yiden) / nrm));   \
}          \
          \
G_COMPLEX(ret) /* const */ operator / (left x, const G_COMPLEX(right)& y) \
{          \
  double den = norm(y);        \
  if (den == 0.0) y.error ("Attempted division by zero.");   \
  return G_COMPLEX(ret)((ret)((x * y.re) / den), -(ret)((x * y.im) / den));\
}          \
          \
G_COMPLEX(ret) /* const */ operator / (const G_COMPLEX(left)& x, right y) \
{          \
  if (y == 0.0) x.error ("Attempted division by zero.");   \
  return G_COMPLEX(ret)((ret)(x.re / y), (ret)(x.im / y));  \
}          \
          \
/* Corrections based on reports from: [e-mail] & [e-mail] */ \
G_COMPLEX(ret) /* const */ pow(const G_COMPLEX(left) &x, const G_COMPLEX(right) &p)\
{          \
  double h = hypot((double)x.re, (double)x.im);    \
  if (h <= 0.0) x.error("attempted power of zero magnitude number.");  \
          \
  double a = atan2((double)x.im, (double)x.re);    \
  double lr = pow(h, (double)p.re);      \
  double li = (double)p.re * a;      \
  if (p.im != 0.0)        \
  {          \
    lr /= exp((double)p.im * a);      \
    li += p.im * log(h);       \
  }          \
  return G_COMPLEX(ret)((ret)(lr * cos(li)), (ret)(lr * sin(li)));  \
}          \
          \


#define G_COMPLEX_DO_BIN_OP_DECL  \
/*G_COMPLEX_BIN_OP_DECL(double,double,double)*/ \
G_COMPLEX_BIN_OP_DECL(double,double,float) \
G_COMPLEX_BIN_OP_DECL(double,float,double) \
G_COMPLEX_BIN_OP_DECL(double,double,int) \
G_COMPLEX_BIN_OP_DECL(double,int,double) \
/*G_COMPLEX_BIN_OP_DECL(float,float,float)*/ \
G_COMPLEX_BIN_OP_DECL(float,float,int)  \
G_COMPLEX_BIN_OP_DECL(float,int,float)  \
/*G_COMPLEX_BIN_OP_DECL(int,int,int)*/

#define G_COMPLEX_DO_BIN_OP_IMP   \
/*G_COMPLEX_BIN_OP_IMP(double,double,double)*/ \
G_COMPLEX_BIN_OP_IMP(double,double,float) \
G_COMPLEX_BIN_OP_IMP(double,float,double) \
G_COMPLEX_BIN_OP_IMP(double,double,int)  \
G_COMPLEX_BIN_OP_IMP(double,int,double)  \
/*G_COMPLEX_BIN_OP_IMP(float,float,float)*/ \
G_COMPLEX_BIN_OP_IMP(float,float,int)  \
G_COMPLEX_BIN_OP_IMP(float,int,float)  \
/*G_COMPLEX_BIN_OP_IMP(int,int,int)*/

G_COMPLEX_DO_BIN_OP_DECL
//#***************************************************************************
//#********************* Should not be needed, casts *************************
//#********************* should take care of these   *************************
//#***************************************************************************
#endif

// Functions to set and test for IEEE NaN's. A complex number is a NaN if either
// the real or imaginary part is a NaN. setNaN sets both the real and imaginary
// parts to NaN.
// <group>
Bool isNaN(const Complex &val);
Bool isNaN(const DComplex &val);
void setNaN(Complex &val);
void setNaN(DComplex &val);
// </group>


#endif