cx_polar.cpp

复数运算库

#define WANT_MATH
#define WANT_STREAM
#include "include.h"
#include "myexcept.h"

#include "cx.h"

#ifdef use_namespace
namespace RBD_COMPLEX {
#endif



// Polar functions


//                      RA = 1
//                  .     |     .
//                    .   |   .
//                      . | .
//    RA = 2  ------------.------------   RA = 0
//                      . | .
//                    .   |   .
//                  .     |     .
//                      RA = 3


Polar::Polar(Real x) : Theta(0)
   { if (x < 0) { R = -x; RA = 2; }  else { R = x; RA = 0; } }


Polar::Polar(Real r, Real theta) : R(r)
{
   Real rra = floor(theta / pi_over_2 + 0.5);
   Theta = theta - rra * (pi_over_2);
   RA = (int)(rra - 4 * floor(rra / 4));
   if (r < 0) { RA.pluseq2(); R = -R; }
}

Polar::Polar(const CX& z)
{
   R = z.cabs();
   if (z.Y > z.X)
   {
      if (z.Y > -z.X) { RA = 1; Theta = ::atan(-z.X / z.Y); }
      else { RA = 2; Theta = ::atan(z.Y / z.X); }
   }
   else
   {
      if (z.Y < -z.X) { RA = 3; Theta = ::atan(-z.X / z.Y); }
      else { RA = 0; Theta = (z.X != 0) ? ::atan(z.Y / z.X) : 0.0; }
   }
}

Polar::Polar(Imag y) : Theta(0)
{
   if (y.Y > 0) { RA = 1; R = y.Y; }
   else if (y.Y < 0) { RA = 3; R = -y.Y; }
   else { RA = 0; R = 0; }
}

void Polar::operator=(Real x)

{
   if (x < 0) { R = -x; RA = 2; Theta = 0; }
   else { R = x; RA = 0; Theta = 0; } }

void Polar::operator=(Imag y)
{
   if (y.Y > 0) { RA = 1; R = y.Y; Theta = 0; }
   else if (y.Y < 0) { RA = 3; R = -y.Y; Theta = 0; }
   else { RA = 0; R = 0; Theta = 0; }
}

Real Polar::real() const
{
   switch (RA)
   {
   case 0: return R * ::cos(Theta);
   case 1: return - R * ::sin(Theta);
   case 2: return - R * ::cos(Theta);
   case 3: return R * ::sin(Theta);
   default: Throw(Runtime_error("real: invalid polar variable")); return 0;
   }
}

Real Polar::imag() const
{
   switch (RA)
   {
   case 0: return R * ::sin(Theta);
   case 1: return R * ::cos(Theta);
   case 2: return - R * ::sin(Theta);
   case 3: return - R * ::cos(Theta);
   default: Throw(Runtime_error("imag: invalid polar variable")); return 0;
   }
}

Real Polar::arg() const
{
   switch (RA)
   {
   case 0: return Theta;
   case 1: return Theta + pi_over_2;
   case 2: return (Theta <= 0) ? Theta + pi : Theta  - pi;
   case 3: return Theta - pi_over_2;
   default: Throw(Runtime_error("arg: invalid polar variable")); return 0;
   }
}

Polar operator*(const Polar& p1, const Polar& p2)
{
   Quadrant RA = p1.RA + p2.RA;
   Real R = p1.R * p2.R;
   Real Theta = p1.Theta + p2.Theta;
   if (Theta > pi_over_4) { Theta -= pi_over_2; ++RA; }
   else if (Theta < -pi_over_4) { Theta += pi_over_2; --RA; }
   return Polar(RA, R, Theta);
}

void CX::operator*=(const Polar& p1) { operator*=(CX(p1)); }

void CX::operator/=(const Polar& p1) { operator/=(CX(p1)); }

void CX::operator+=(const Polar& p1) { operator+=(CX(p1)); }

void CX::operator-=(const Polar& p1) { operator-=(CX(p1)); }



Polar polar_exp(const CX& z) { return Polar( ::exp(z.X), z.Y); }


Polar operator/(const Polar& p1, const Polar& p2)
{
   Quadrant RA = p1.RA - p2.RA;
   Real R = p1.R / p2.R;
   Real Theta = p1.Theta - p2.Theta;
   if (Theta > pi_over_4) { Theta -= pi_over_2; ++RA; }
   else if (Theta < -pi_over_4) { Theta += pi_over_2; --RA; }
   return Polar(RA, R, Theta);
}

CX operator+(const Polar& p1, const Polar& p2) { return CX(p1) + CX(p2); }

CX operator-(const Polar& p1, const Polar& p2) { return CX(p1) - CX(p2); }

CX::CX(const Polar& p)
{
   switch (p.RA)
   {
   case 0: X = p.R * ::cos(p.Theta); Y = p.R * ::sin(p.Theta); return;
   case 1: X = - p.R * ::sin(p.Theta); Y = p.R * ::cos(p.Theta); return;
   case 2: X = - p.R * ::cos(p.Theta); Y = - p.R * ::sin(p.Theta); return;
   case 3: X = p.R * ::sin(p.Theta); Y = - p.R * ::cos(p.Theta); return;
   default: Throw(Runtime_error("CX: invalid polar variable"));
   }
}

void CX::operator=(const Polar& p)
{
   switch (p.RA)
   {
   case 0: X = p.R * ::cos(p.Theta); Y = p.R * ::sin(p.Theta); return;
   case 1: X = - p.R * ::sin(p.Theta); Y = p.R * ::cos(p.Theta); return;
   case 2: X = - p.R * ::cos(p.Theta); Y = - p.R * ::sin(p.Theta); return;
   case 3: X = p.R * ::sin(p.Theta); Y = - p.R * ::cos(p.Theta); return;
   default: Throw(Runtime_error("=: invalid polar variable"));
   }
}

CX operator+(const Polar& p1, const CX& z2) { return CX(p1) + z2; }
CX operator-(const Polar& p1, const CX& z2) { return CX(p1) - z2; }
CX operator+(const CX& z1, const Polar& p2) { return z1 + CX(p2); }
CX operator-(const CX& z1, const Polar& p2) { return z1 - CX(p2); }


Polar operator*(const Polar& p1, Real x2)
{
   if (x2 > 0) return Polar(p1.RA, p1.R * x2, p1.Theta);
   else if (x2 < 0) return Polar(p1.RA.plus2(), - p1.R * x2, p1.Theta);
   else return Polar(0, 0, 0);
}

Polar operator/(const Polar& p1, Real x2)
{
   if (x2 >= 0) return Polar(p1.RA, p1.R / x2, p1.Theta);
   else return Polar(p1.RA.plus2(), - p1.R / x2, p1.Theta);
}

Polar operator/(Real x1, const Polar& p2)
{
   if (x1 >= 0) return Polar(p2.RA.h_reflect(), x1 / p2.R, -p2.Theta);
   else return Polar(p2.RA.v_reflect(), - x1 / p2.R, -p2.Theta);
}

CX operator+(const Polar& p1, Real x2) { return CX(p1) + x2; }
CX operator-(const Polar& p1, Real x2) { return CX(p1) - x2; }
CX operator+(Real x1, const Polar& p2) { return x1 + CX(p2); }
CX operator-(Real x1, const Polar& p2) { return x1 - CX(p2); }

Polar operator*(const Polar& p1, Imag y2)
{
   if (y2.Y > 0) return Polar(p1.RA.plus1(), p1.R * y2.Y, p1.Theta);
   else if (y2.Y < 0) return Polar(p1.RA.plus3(), - p1.R * y2.Y, p1.Theta);
   else return Polar(0, 0, 0);
}

Polar operator/(const Polar& p1, Imag y2)
{
   if (y2.Y >= 0) return Polar(p1.RA.plus3(), p1.R / y2.Y, p1.Theta);
   else return Polar(p1.RA.plus1(), - p1.R / y2.Y, p1.Theta);
}

CX operator+(const Polar& p1, Imag y2) { return CX(p1) + y2; }
CX operator-(const Polar& p1, Imag y2) { return CX(p1) - y2; }

Polar operator/(Imag y1, const Polar& p2)
{
   if (y1.Y >= 0) return Polar(p2.RA.h_reflect_plus1(), y1.Y / p2.R, -p2.Theta);
   else return Polar(p2.RA.h_reflect_minus1(), - y1.Y / p2.R, -p2.Theta);
}

CX operator+(Imag y1, const Polar& p2) { return y1 + CX(p2); }
CX operator-(Imag y1, const Polar& p2) { return y1 - CX(p2); }

CX operator+(const Polar& p1, ImaginaryUnit) { return CX(p1) + _I_; }
CX operator-(const Polar& p1, ImaginaryUnit) { return CX(p1) - _I_; }
CX operator+(ImaginaryUnit, const Polar& p2) { return _I_ + CX(p2); }
CX operator-(ImaginaryUnit, const Polar& p2) { return _I_ - CX(p2); }


void Polar::operator*=(const Polar& p)
{
   RA += p.RA; R *= p.R; Theta += p.Theta;
   if (Theta > pi_over_4) { Theta -= pi_over_2; ++RA; }
   else if (Theta < -pi_over_4) { Theta += pi_over_2; --RA; }
}

void Polar::operator/=(const Polar& p)
{
   RA -= p.RA; R /= p.R; Theta -= p.Theta;
   if (Theta > pi_over_4) { Theta -= pi_over_2; ++RA; }
   else if (Theta < -pi_over_4) { Theta += pi_over_2; --RA; }
}

void Polar::operator+=(const Polar& p) { *this = Polar(*this + p); }
void Polar::operator-=(const Polar& p) { *this = Polar(*this - p); }

void Polar::operator*=(const CX& z) { operator*=(Polar(z)); }
void Polar::operator/=(const CX& z) { operator/=(Polar(z)); }
void Polar::operator+=(const CX& z) { *this = Polar(CX(*this) + z); }
void Polar::operator-=(const CX& z) { *this = Polar(CX(*this) - z); }

void Polar::operator*=(Real x)
   { if (x >= 0) { R *= x; } else { R *= (-x); RA.pluseq2(); } }

void Polar::operator/=(Real x)
   { if (x >= 0) { R /= x; } else { R /= (-x); RA.pluseq2(); } }

void Polar::operator+=(Real x) { *this = Polar(CX(*this) + x); }
void Polar::operator-=(Real x) { *this = Polar(CX(*this) - x); }

void Polar::operator*=(Imag y)
   { if (y.Y >= 0) { R *= y.Y; ++RA; } else { R *= (-(y.Y)); --RA; } }

void Polar::operator/=(Imag y)
   { if (y.Y >= 0) { R /= y.Y; --RA; } else { R /= (-(y.Y)); ++RA; } }

void Polar::operator+=(Imag y) { *this = Polar(CX(*this) + y); }
void Polar::operator-=(Imag y) { *this = Polar(CX(*this) - y); }

void Polar::operator+=(ImaginaryUnit) { *this = Polar(CX(*this) + _I_); }
void Polar::operator-=(ImaginaryUnit) { *this = Polar(CX(*this) - _I_); }

void Polar::AssertIsValid() const
{
   if (RA < 0 || RA > 3 || R < 0 || ::fabs(Theta) > pi_over_4)
      { Throw(Runtime_error("Assert fails: invalid polar variable")); }
}




CX log(const Polar& p)
{
   Real lpr = ::log(p.R);
   switch (p.RA)
   {
   case 0: return CX(lpr, p.Theta);
   case 1: return CX(lpr, p.Theta + pi_over_2);
   case 2: return CX(lpr, (p.Theta > 0) ? p.Theta - pi : p.Theta + pi);
   case 3: return CX(lpr, p.Theta - pi_over_2);
   default: Throw(Runtime_error("log: invalid polar variable")); return 0;
   }
}

Polar sqrt(const Polar& p)
{
   Real spr = ::sqrt(p.R); Real pt2 = p.Theta / 2;
   switch (p.RA)
   {
   case 0: return Polar(0, spr, pt2);
   case 1:
      if (pt2 <= 0) return Polar(0, spr, pt2 + pi_over_4);
      else return Polar(1, spr, pt2 - pi_over_4);
   case 2:
      if (pt2 <= 0) return Polar(1, spr, pt2);
      else return Polar(3, spr, pt2);
   case 3:
      if (pt2 <= 0) return Polar(3, spr, pt2 + pi_over_4);
      else return Polar(0, spr, pt2 - pi_over_4);
   default: Throw(Runtime_error("sqrt: invalid polar variable")); return 0;
   }
}


Polar pow(const Polar& p1, int n2)
{
   Real Theta = p1.Theta * n2; int p1RA = p1.RA;
   int RA = (n2 >= 0) ? (n2 & 3) * p1RA : (4 - ((-n2) & 3)) * p1RA;
   Real rra = floor(Theta / pi_over_2 + 0.5);
   Theta -= rra * (pi_over_2);
   RA += (int)(rra - 4 * floor(rra / 4));
   return Polar(RA & 3, ipow(p1.R, n2), Theta);
}



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

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

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


// equality operators - must allow for zero values when Theta can be anything

bool operator==(const Polar& p1, const Polar& p2)
{
   return (p1.R == 0.0 && p2.R == 0.0)
     || (p1.RA == p2.RA && p1.R == p2.R && p1.Theta == p2.Theta);
}

bool operator==(const Polar& p1, Real r2)
{
   return (r2 == 0.0 && p1.R == 0.0)
     || (r2 > 0.0 && p1.RA == 0 && p1.R == r2 && p1.Theta == 0.0)
     || (r2 < 0.0 && p1.RA == 2 && p1.R == -r2 && p1.Theta == 0.0);
}

bool operator==(const Polar& p1, Imag y2)
{
   return (y2.Y == 0.0 && p1.R == 0.0)
     || (y2.Y > 0.0 && p1.RA == 1 && p1.R == y2.Y && p1.Theta == 0.0)
     || (y2.Y < 0.0 && p1.RA == 3 && p1.R == -y2.Y && p1.Theta == 0.0);
}

#ifdef use_namespace
}
#endif