c++复数类库.txt

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

          \
inline int operator < (const G_COMPLEX(type) x, const G_COMPLEX(type) y) { \
  return norm(x) < norm(y);       \
}          \
          \
inline int operator < (const G_COMPLEX(type) x, const type y) {   \
  return norm(x) < y*y;        \
}          \
          \
inline int operator < (const type x, const G_COMPLEX(type) y) {   \
  return x*x < norm(y);        \
}          \
inline G_COMPLEX(type) min (const G_COMPLEX(type) x, const G_COMPLEX(type) y)   \
{          \
  return x <= y ? x : y;                                          \
}          \
          \
inline G_COMPLEX(type) max (const G_COMPLEX(type) x, const G_COMPLEX(type) y)   \
{          \
  return x >= y ? x : y;                                          \
}          \
          \
          \
/*          \
* NEGATION OPERATORS        \
*/          \
inline G_COMPLEX(type) operator - (const G_COMPLEX(type) x)   \
{          \
  return G_COMPLEX(type)(-x.re, -x.im);     \
}          \
          \
inline G_COMPLEX(type) conj(const G_COMPLEX(type) x)    \
{          \
  return G_COMPLEX(type)(x.re, -x.im);     \
}          \
          \
/*          \
* SIMPLE NUMERIC OPERATIONS       \
*/          \
inline G_COMPLEX(type) operator + (const G_COMPLEX(type) x, const G_COMPLEX(type) y)\
{          \
  return G_COMPLEX(type)(x.re + y.re, x.im + y.im);  \
}          \
          \
inline G_COMPLEX(type) operator + (const G_COMPLEX(type) x, type y)  \
{          \
  return G_COMPLEX(type)(x.re + y, x.im);    \
}          \
          \
inline G_COMPLEX(type) operator + (type x, const G_COMPLEX(type) y)  \
{          \
  return G_COMPLEX(type)(x + y.re, y.im);    \
}          \
          \
inline G_COMPLEX(type) operator - (const G_COMPLEX(type) x, const G_COMPLEX(type) y)\
{          \
  return G_COMPLEX(type)(x.re - y.re, x.im - y.im);  \
}          \
          \
inline G_COMPLEX(type) operator - (const G_COMPLEX(type) x, type y)  \
{          \
  return G_COMPLEX(type)(x.re - y, x.im);    \
}          \
          \
inline G_COMPLEX(type) operator - (type x, const G_COMPLEX(type) y)  \
{          \
  return G_COMPLEX(type)(x - y.re, -y.im);    \
}          \
          \
inline G_COMPLEX(type) operator * (const G_COMPLEX(type) x, const G_COMPLEX(type) y)\
{          \
  return G_COMPLEX(type)(x.re * y.re - x.im * y.im,   \
                 x.re * y.im + x.im * y.re);   \
}          \
          \
inline G_COMPLEX(type) operator * (const G_COMPLEX(type) x, type y)  \
{          \
  return G_COMPLEX(type)(x.re * y, x.im * y);    \
}          \
          \
inline G_COMPLEX(type) operator * (type x, const G_COMPLEX(type) y)  \
{          \
  return G_COMPLEX(type)(x * y.re, x * y.im);    \
}

// </group>

/*
* --   --   --   --   --   -- Implementation --   --   --   --   --   --
*/


#define G_COMPLEXimplement2(type,CTOR)      \
          \
CTOR          \
          \
void  G_COMPLEX(type)::error(const char* msg) const    \
{          \
  (*Complex_error_handler)(msg);      \
}          \
          \
/* from [e-mail] */      \
G_COMPLEX(type) /* const */ operator / (const G_COMPLEX(type) x, const G_COMPLEX(type) 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(type)((type)((xrden * yrden + xiden * yiden) / nrm),  \
                 (type)((xiden * yrden - xrden * yiden) / nrm));  \
}          \
          \
G_COMPLEX(type)& G_COMPLEX(type)::operator /= (const G_COMPLEX(type) y)  \
{          \
  double den = fabs((double)y.re) + fabs((double)y.im);   \
  if (den == 0.0) error ("Attempted division by zero.");   \
  double xrden = re / den;       \
  double xiden = im / den;       \
  double yrden = y.re / den;      \
  double yiden = 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) /* const */ operator / (type x, const G_COMPLEX(type) y) \
{          \
  double den = norm(y);        \
  if (den == 0.0) y.error ("Attempted division by zero.");   \
  return G_COMPLEX(type)((type)((x * y.re) / den), -(type)((x * y.im) / den));\
}          \
          \
G_COMPLEX(type) /* const */ operator / (const G_COMPLEX(type) x, type y) \
{          \
  if (y == 0.0) x.error ("Attempted division by zero.");   \
  return G_COMPLEX(type)((type)(x.re / y), (type)(x.im / y));  \
}          \
          \
          \
G_COMPLEX(type)& G_COMPLEX(type)::operator /= (type y)    \
{          \
  if (y == 0.0) error ("Attempted division by zero.");    \
  re /= y;  im /= y;        \
  return *this;         \
}          \
          \
          \
G_COMPLEX(type) /* const */ exp(const G_COMPLEX(type) x)   \
{          \
  double r = exp((double)x.re);      \
  return G_COMPLEX(type)((type)(r * cos((double)x.im)),    \
                 (type)(r * sin((double)x.im)));    \
}          \
          \
G_COMPLEX(type) /* const */ cosh(const G_COMPLEX(type) x)   \
{          \
  return G_COMPLEX(type)((type)(cos((double)x.im) * cosh((double)x.re)),\
                 (type)(sin((double)x.im) * sinh((double)x.re))); \
}          \
          \
G_COMPLEX(type) /* const */ sinh(const G_COMPLEX(type) x)   \
{          \
  return G_COMPLEX(type)((type)(cos((double)x.im) * sinh((double)x.re)),\
                 (type)(sin((double)x.im) * cosh((double)x.re))); \
}          \
          \
G_COMPLEX(type) /* const */ cos(const G_COMPLEX(type) x)   \
{          \
  return G_COMPLEX(type)((type)(cos((double)x.re) * cosh((double)x.im)),\
                 (type)(-sin((double)x.re) * sinh((double)x.im))); \
}          \
          \
G_COMPLEX(type) /* const */ sin(const G_COMPLEX(type) x)   \
{          \
  return G_COMPLEX(type)((type)(sin((double)x.re) * cosh((double)x.im)),\
                 (type)(cos((double)x.re) * sinh((double)x.im))); \
}          \
          \
G_COMPLEX(type) /* const */ log(const G_COMPLEX(type) x)   \
{          \
  double h = hypot((double)x.re, (double)x.im);    \
  if (h <= 0.0) x.error("attempted log of zero magnitude number.");  \
  return G_COMPLEX(type)((type)log(h), (type)(atan2((double)x.im, (double) x.re)));\
}          \
          \
/* Corrections based on reports from: [e-mail] & [e-mail] */ \
G_COMPLEX(type) /* const */ pow(const G_COMPLEX(type) x, const G_COMPLEX(type) 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(type)((type)(lr * cos(li)), (type)(lr * sin(li)));  \
}          \
          \
G_COMPLEX(type) /* const */ pow(const G_COMPLEX(type) x, double p)  \
{          \
  double h = hypot((double)x.re, (double)x.im);    \
  if (h <= 0.0) x.error("attempted power of zero magnitude number.");  \
  double lr = pow(h, p);       \
  double a = atan2((double)x.im, (double)x.re);    \
  double li = p * a;        \
  return G_COMPLEX(type)((type)(lr * cos(li)), (type)(lr * sin(li)));  \
}          \
          \
          \
G_COMPLEX(type) /* const */ sqrt(const G_COMPLEX(type) x)   \
{          \
  if (x.re == 0.0 && x.im == 0.0)     \
    return G_COMPLEX(type)((type)0, (type)0);     \
  else          \
  {          \
    double s = sqrt((fabs((double)x.re) + hypot((double)x.re, (double)x.im)) * 0.5);\
    double d = (x.im / s) * 0.5;      \
    if (x.re > (type)0)       \
      return G_COMPLEX(type)((type) s, (type) d);    \
    else if (x.im >= (type) 0)      \
      return G_COMPLEX(type)((type)d, (type)s);     \
    else         \
      return G_COMPLEX(type)((type)-d, (type)-s);    \
  }          \
}          \
          \
G_COMPLEX(type) log10(const G_COMPLEX(type) x)     \
{          \
    return log(x)*(type)C::log10e;      \
}          \
          \
          \
G_COMPLEX(type) /* const */ pow(const G_COMPLEX(type) x, int p)   \
{          \
  if (p == 0)         \
    return G_COMPLEX(type)((type)1, (type)0);     \
  else if (x == (type) 0)       \
    return G_COMPLEX(type)((type)0, (type)0);     \
  else          \
  {          \
    G_COMPLEX(type) res((type)1, (type)0);     \
    G_COMPLEX(type) b = x;       \
    if (p < 0)         \
    {          \
      p = -p;         \
      b = (type) 1 / b;        \
    }          \
          \
    do {         \
      if ( p & 1 )        \
        res *= b;        \
    } while ( ( p >>= 1 ) != 0 && ((b *= b),1) );    \
    return res;         \
          \
  }          \
}          \
          \
G_COMPLEX(type) cube(const G_COMPLEX(type) val)     \
{          \
    G_COMPLEX(type) retval(val);      \
    retval *= val;        \
    retval *= val;        \
    return retval;        \
}          \
          \
ostream& operator << (ostream& s, const G_COMPLEX(type) x)   \
{          \
  return s << "(" << x.re << "," << x.im << ")" ;   \
}          \
          \
istream& operator >> (istream& s, G_COMPLEX(type)& x)    \
{          \
  if (!s.good())        \
  {          \
    return s;         \
  }          \
  type r, i;         \
  char ch;         \
  s >> ws;         \
  s.get(ch);         \
  if (ch == '(')        \
  {          \
    s >> r;         \
    s >> ws;         \
    s.get(ch);         \
    if (ch == ',')        \
    {          \
      s >> i;         \
      s >> ws;         \
      s .get(ch);        \
    }          \
    else         \
      i = 0;         \
    if (ch != ')')        \
      s.clear(ios::failbit);       \
  }          \
  else          \
  {          \
    s.putback(ch);        \
    s >> r;         \
    i = 0;         \
  }          \
  x = G_COMPLEX(type)(r, i);       \
  return s;         \
}

#define G_COMPLEX_CAST_OP(type)       \
   operator G_COMPLEX(type)() const {return G_COMPLEX(type)((type)re,(type)im);}

#define G_COMPLEX_ASSIGN_OP_DECL(type,from)     \
  G_COMPLEX(type)& operator=(from y) { re = (type) y; im = (type) 0; return *this;} \
  G_COMPLEX(type)& operator=(const G_COMPLEX(from)& y);
#define G_COMPLEX_ASSIGN_OP_DECLS(type,from)     \
  G_COMPLEX(type)& operator=(from y) { re = (type) y; im = (type) 0; return *this;}
#define G_COMPLEX_ASSIGN_OP_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;}


#define G_COMPLEX_CTOR_OP_DECL(type,intype)     \
   G_COMPLEX(type)(const G_COMPLEX(intype) &v);