cx.cpp

复数运算库

#define WANT_MATH
#define WANT_STREAM

#include "cx.h"

#ifdef use_namespace
namespace RBD_COMPLEX {
#endif


ImaginaryUnit _I_;

Real pi, pi_times_2, pi_over_2, pi_over_4;


Real CX::cabs() const
{
   // reduce chances of floating point overflow

   Real x = ::fabs(X); Real y = ::fabs(Y);
   if (x > y) { Real r = y / x; return x * ::sqrt(1.0 + r * r); }
   else if (y == 0) return 0;
   else { Real r = x / y; return y * ::sqrt(1.0 + r * r); }
}

Real CX::arg() const { return atan2(Y, X); }

Real Imag::arg() const
   { return (Y > 0.0) ? pi_over_2 : (Y < 0.0) ? - pi_over_2 : 0; }


CX operator/(const CX& w1, const CX& w2)
{
   // reduce chances of floating point overflow
   if (::fabs(w2.X)>=::fabs(w2.Y))
   {
      Real r = w2.Y/w2.X; Real d = w2.X * (1.0 + r * r);
      return CX( (w1.X+w1.Y*r)/d, (w1.Y-w1.X*r)/d );
   }
   else
   {
      Real r = w2.X/w2.Y; Real d = w2.Y * (1.0 + r * r);
      return CX( (w1.X*r+w1.Y)/d, (w1.Y*r-w1.X)/d );
   }
}

CX operator/(Real x1, const CX& w)
{
   if (::fabs(w.X) >= ::fabs(w.Y))
   {
      Real r = w.Y/w.X; Real d = x1 / (w.X * (1.0 + r * r));
      return CX( d, -d*r );
   }
   else
   {
      Real r = w.X/w.Y; Real d = x1 / (w.Y * (1.0 + r * r));
      return CX( d*r, -d );
   }
}

CX operator/(Imag y1, const CX& w)
{
   if (::fabs(w.X) >= ::fabs(w.Y))
   {
      Real r = w.Y/w.X; Real d = y1.Y / (w.X * (1.0 + r * r));
      return CX( d*r, d );
   }
   else
   {
      Real r = w.X/w.Y; Real d = y1.Y / (w.Y * (1.0 + r * r));
      return CX( d, d*r );
   }
}

CX operator/(ImaginaryUnit, const CX& w)
{
   if (::fabs(w.X) >= ::fabs(w.Y))
   {
      Real r = w.Y/w.X; Real d = 1.0 / (w.X * (1.0 + r * r));
      return CX( d*r, d );
   }
   else
   {
      Real r = w.X/w.Y; Real d = 1.0 / (w.Y * (1.0 + r * r));
      return CX( d, d*r );
   }
}


CX exp(const CX& z) { return ::exp(z.X) * CX(::cos(z.Y), ::sin(z.Y)); }


CX log(const CX& z)
   { return CX(0.5 * ::log(z.X*z.X+z.Y*z.Y), ::atan2(z.Y,z.X)); }

CX sqrt(const CX& z)
{
   Real x = ::fabs(z.X); Real y = ::fabs(z.Y);
   if (x > y)
   {
      Real r = y / x; Real s = ::sqrt(1.0 + r * r) + 1.0;
      x = ::sqrt(x * s / 2); y = ::sqrt(y * r / (s * 2));
   }
   else if (y == 0) return CX(0, 0);
   else
   {
      Real r = x / y; Real s = ::sqrt(1.0 + r * r) + r;
      x = ::sqrt(y * s / 2); y = ::sqrt(y / (s * 2));
   }
   if (z.Y >= 0) { if (z.X >= 0) return CX(x, y); else return CX(y, x); }
   else { if (z.X >= 0) return CX(x, -y); else return CX(y, -x); }
}

CX sin(const CX& z)
   { return CX(::sin(z.X) * ::cosh(z.Y), ::cos(z.X) * ::sinh(z.Y)); }

CX cos(const CX& z)
   { return CX(::cos(z.X) * ::cosh(z.Y), - ::sin(z.X) * ::sinh(z.Y)); }

CX tan(const CX& z)
   { return CX(::tan(z.X), ::tanh(z.Y)) / CX(1, -::tan(z.X) * ::tanh(z.Y)); }

CX sinh(const CX& z) { return sin(_I_ * z) / _I_; }

CX cosh(const CX& z) { return cos(_I_ * z); }

CX tanh(const CX& z) { return tan(_I_ * z) / _I_; }



CX exp(Imag z) { return CX(::cos(z.Y), ::sin(z.Y)); }

CX log(Imag z)
{
   if (z.Y > 0.0) return CX(::log(z.Y), pi_over_2);
   else return CX(::log(-z.Y), -pi_over_2);
}

CX sqrt(Imag y)
{
   if (y.Y >= 0) { Real sy = ::sqrt(y.Y / 2.0); return CX(sy, sy); }
   else { Real sy = ::sqrt(-y.Y / 2.0); return CX(sy, -sy); }
}

CX pow(const CX& z, int n2)
{
   switch (n2)
   {
   case 0: return CX(1.0);
   case 1: return z;
   case 2: return square(z);
   case 3: return square(z) * z;
   case 4: { CX z2 = square(z); return square(z2); }
   case 5: { CX z2 = square(z); return square(z2) * z; }
   case 6: { CX z2 = square(z); return square(z2) * z2; }
   case 7: { CX z2 = square(z); return square(z2) * z2 * z; }
   case 8:
      { CX z2 = square(z); CX z4 = square(z2); return square(z4); }
   case 9:
      { CX z2 = square(z); CX z4 = square(z2); return square(z4) * z; }
   case 10:
      { CX z2 = square(z); CX z4 = square(z2); return square(z4) * z2; }
   case 11:
      { CX z2 = square(z); CX z4 = square(z2); return square(z4) * z2 * z; }
   case 12:
      { CX z2 = square(z); CX z4 = square(z2); return square(z4) * z4; }
   default:
      if (n2 < 0 && n2 >= -12)
      {
         if (z == 0)
            Throw(Runtime_error("pow: first arg 0, second (integer) arg < 0"));
         return pow(1.0 / z, -n2);
      }
      return pow(z, (Real)n2);
   }
}

CX pow(const CX& z1, Real r2)
{
   if (z1 == 0)
   {
      if (r2 > 0) return 0.0;
      else
         Throw(Runtime_error("pow: first arg 0, second arg real part <= 0"));
   }
   return exp(r2 * log(z1));
}

CX pow(const CX& z1, Imag y2)
{
   if (z1 == 0)
      Throw(Runtime_error("pow: first arg 0, second arg imaginary"));
   return exp(y2 * log(z1));
}

CX pow(const CX& z1, const CX& z2)
{
   if (z1 == 0)
   {
      if (z2.X > 0) return 0.0;
      else
         Throw(Runtime_error("pow: first arg 0, second arg real part <= 0"));
   }
   return exp(z2 * log(z1));
}

CX pow(Real r1, const CX& z2)
{
   if (r1 == 0)
   {
      if (z2.X > 0) return 0.0;
      else
         Throw(Runtime_error("pow: first arg 0, second arg real part <= 0"));
   }
   else if (r1 > 0) return exp(z2 * ::log(r1));
   return exp(z2 * (_I_ * pi + ::log(-r1)));
}

CX pow(Imag y1, const CX& z2)
{
   if (y1 == 0)
   {
      if (z2.X > 0) return 0.0;
      else
         Throw(Runtime_error("pow: first arg 0, second arg real part <= 0"));
   }
   return exp(z2 * log(y1));
}

CX pow(Imag y1, int n2)
{
   Real R = ipow(y1.Y, n2);
   int RA = (n2 >= 0) ? n2 & 3 : (4 - ((-n2) & 3)) & 3;
   switch (RA)
   {
   case 0: return CX(R);
   case 1: return CX(0,R);
   case 2: return CX(-R);
   case 3: return CX(0,-R);
   }
   return 0;
}


CX pow(Imag y1, Real r2)
{
   if (y1 == 0)
   {
      if (r2 > 0) return 0.0;
      else
         Throw(Runtime_error("pow: first arg 0, second arg real part <= 0"));
   }
   return exp(r2 * log(y1));
}

CX pow(Imag y1, Imag y2)
{
   if (y1 == 0)
      Throw(Runtime_error("pow: first arg 0, second arg imaginary"));
   return exp(y2 * log(y1));
}

CX pow(Real r1, Imag y2)
{
   if (r1 == 0)
      Throw(Runtime_error("pow: first arg 0, second arg imaginary"));
   else if (r1 > 0) return exp(y2 * ::log(r1));
   return exp(y2 * (_I_ * pi + ::log(-r1)));
}



// won't need this when integer version of pow is available

Real ipow(Real x, int n)
{
   switch (n)
   {
   case 0:  return 1.0;
   case 1:  return x;
   case 2:  return x*x;
   case 3:  return x*x*x;
   case 4:  { Real x2 = x*x; return x2*x2; }
   case 5:  { Real x2 = x*x; return x2*x2 * x; }
   case 6:  { Real x2 = x*x; return x2*x2 * x2; }
   case 7:  { Real x2 = x*x; return x2*x2 * x2 * x; }
   case 8:  { Real x2 = x*x; Real x4 = x2*x2; return x4*x4; }
   case 9:  { Real x2 = x*x; Real x4 = x2*x2; return x4*x4 * x; }
   case 10: { Real x2 = x*x; Real x4 = x2*x2; return x4*x4 * x2; }
   case 11: { Real x2 = x*x; Real x4 = x2*x2; return x4*x4 * x2 * x; }
   case 12: { Real x2 = x*x; Real x4 = x2*x2; return x4*x4 * x4; }
   default:
      if (n < 0 && n >= -12) return ipow(1.0 / x, -n);
      return ::pow(x, (Real)n);
   }
}

ComplexPackageInitialise::ComplexPackageInitialise()
{
   if ( !pi)
   {
      pi = 3.141592653589793238462643;
      pi_times_2 = pi * 2;
      pi_over_2 = pi / 2;
      pi_over_4 = pi / 4;
   }
}


#ifdef use_namespace
}
#endif