#define WANT_MATH#define WANT_STREAM#include "cx.h"#ifdef use_namespacenamespace RBD_COMPLEX {#endifImaginaryUnit _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 availableReal 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