cxtest1.cpp

复数运算库

#define WANT_MATH
#define WANT_STREAM

#include "include.h"
#include "myexcept.h"
#include "array1.h"
#include "cx.h"
#include "cxtest.h"

#ifdef use_namespace
using namespace RBD_COMMON;             // needed for VC++
using namespace RBD_COMPLEX;
#endif




static Real real(Real x) { return x; }
static Real imag(Real) { return 0.0; }

static Real OverallMax = 0.0;



template <class T1, class T2>
void PrintMaxDiff(const array1<T1>& A, const array1<T2>& B)
{
   int n = A.size();
   Real maxdiff = 0.0;
   for (int i = 0; i < n; i++)
   {
      Real diff = fabs(A(i) - B(i));
      if (maxdiff < diff) maxdiff = diff;
      if (diff > 0.0000001)
      {
         cout << i << ": ";
         cout << real(A(i)) << " " << imag(A(i)) << "; ";
         cout << real(B(i)) << " " << imag(B(i)) << " -> " << diff << endl;
      }
   }
   cout << maxdiff << endl;
   if (OverallMax < maxdiff) OverallMax = maxdiff;
}

template <class T1, class T2>
void PrintMaxRelDiff(const array1<T1>& A, const array1<T2>& B)
{
   int n = A.size();
   Real maxdiff = 0.0;
   for (int i = 0; i < n; i++)
   {
      Real diff = fabs(A(i) - B(i));
      Real sum = fabs(A(i)) + fabs(B(i));
      if (sum != 0.0) diff /= sum;
      if (maxdiff < diff) maxdiff = diff;
      if (diff > 0.0000001)
      {
         cout << i << ": ";
         cout << real(A(i)) << " " << imag(A(i)) << "; ";
         cout << real(B(i)) << " " << imag(B(i)) << " -> " << diff << endl;
      }
   }
   cout << maxdiff << endl;
   if (OverallMax < maxdiff) OverallMax = maxdiff;
}

template <class T1, class T2>
void operator+=(array1<T1>& A, const array1<T2>& B)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   for (int j = 0; j < nc; j++) A(j) += B(j);
}

template <class T1, class T2>
void operator-=(array1<T1>& A, const array1<T2>& B)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   for (int j = 0; j < nc; j++) A(j) -= B(j);
}

template <class T1, class T2>
void operator*=(array1<T1>& A, const array1<T2>& B)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   for (int j = 0; j < nc; j++) A(j) *= B(j);
}

template <class T1, class T2>
void operator/=(array1<T1>& A, const array1<T2>& B)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   for (int j = 0; j < nc; j++) A(j) /= B(j);
}


void AssertIsValid(const array1<Polar>& p)
{
   int n = p.size();
   for (int i = 0; i < n; i++) p(i).AssertIsValid();
}





Real cxtest1()
{
   cout << "begin unary tests" << endl;

   int n = 441;
   array1<CX>    CA(n), CB(n), CC(n), CD(n);
   array1<CX>    CE(n), CF(n), CG(n), CH(n);
   array1<Polar> PA(n), PB(n), PC(n), PD(n);
   array1<Real>  RA(n), RB(n), RC(n), RD(n);
   array1<Real>  RE(n), RF(n), RG(n), RH(n);
   array1<Imag>  IA(n), IB(n), IC(n), ID(n);
   array1<Imag>  IE(n), IF(n), IG(n), IH(n);


   // load CA so we get a lot of combinations of values
   // including pure real and pure imaginary and equal real and imaginary parts

      int j, k;
      CA(0) = 0.0;
      CA(1) = 1.0;
      CA(2) = -1.0;
      CA(3) = 1.9;
      CA(4) = -1.9;
      CA(5) = 0.75;
      CA(6) = -0.75;
      CA(7) = 2.35;
      CA(8) = -2.35;
      CA(9) = 4.33;
      CA(10) = -4.33;
      CA(11) = 0.001;
      CA(12) = -0.001;
      CA(13) = 0.25;
      CA(14) = -0.25;
      CA(15) = 0.5;
      CA(16) = -0.5;
      CA(17) = 2;
      CA(18) = -2;
      CA(19) = 7;
      CA(20) = -7;
      for (j = 1; j <= 20; j++) CA(j + 20) = Imag(CA(j).real());
      for (j = 1; j <= 20; j++) for (k = 1; k <= 20; k++)
         CA(j + 20 * k + 20) = CA(j) + CA(k + 20);


   // get Real part and Imaginary part and make Polar version

   for (j = 0; j < n; j++)
   {
      RA(j) = real(CA(j)); IA(j) = Imag(imag(CA(j))); PA(j) = CA(j);
      RB(j) = imag(CA(j));
   }

   PrintMaxDiff(PA, CA); AssertIsValid(PA);



   cout << "check pow(CX, int) against pow(CX, Real)" << endl;
   for (k = -14; k <= 14; k++)
   {
      for (j = 0; j < n; j++)
      {
         if (CA(j) != 0 || k > 0)
            { CC(j) = pow(CA(j),k); CD(j) = pow(CA(j),(Real)k); }
         else { CC(j) = 0.0; CD(j) = 0.0; }
      }
      PrintMaxRelDiff(CC,CD);
   }

   cout << "check pow(Imag, int) against pow(Imag, Real)" << endl;
   for (k = -14; k <= 14; k++)
   {
      for (j = 0; j < n; j++)
      {
         if (IA(j) != 0 || k > 0)
            { CC(j) = pow(IA(j),k); CD(j) = pow(IA(j),(Real)k); }
         else { CC(j) = 0.0; CD(j) = 0.0; }
      }
      PrintMaxRelDiff(CC,CD);
   }

   cout << "check pow(Polar, int) against pow(Polar, Real)" << endl;
   for (k = -14; k <= 14; k++)
   {
      for (j = 0; j < n; j++)
      {
         if (PA(j) != 0 || k > 0)
            { PC(j) = pow(PA(j),k); PD(j) = pow(PA(j),(Real)k); }
         else { PC(j) = 0.0; PD(j) = 0.0; }
      }
      PrintMaxRelDiff(PC,PD);
   }

   cout << "spot check pow(CX, Real)" << endl;
   for (j = 0; j < n; j++) { CC(j) = pow(CA(j), 0.5); CD(j) = sqrt(CA(j)); }
   PrintMaxRelDiff(CC,CD);

   for (j = 0; j < n; j++) { CC(j) = pow(CA(j), 2.0); CD(j) = square(CA(j)); }
   PrintMaxRelDiff(CC,CD);

   for (j = 0; j < n; j++)
   {
      if (CA(j) != 0) { CC(j) = pow(CA(j),-1.0); CD(j) = 1.0 / CA(j); }
      else { CC(j) = 0.0; CD(j) = 0.0; }
   }
   PrintMaxRelDiff(CC,CD);

   cout << "check ipow(Real, int) against pow(Real, Real)" << endl;
   for (k = -14; k <= 14; k++)
   {
      for (j = 0; j < n; j++)
      {
         if (RA(j) != 0 || k > 0)
            { RE(j) = ipow(RA(j),k); RF(j) = pow(RA(j),(Real)k); }
         else { RE(j) = 0.0; RF(j) = 0.0; }
      }
      PrintMaxRelDiff(RE,RF);
   }

   cout << "check exponential for complex" << endl;
   for (j = 0; j < n; j++)
   {
      CC(j) = exp(CA(j));
      CD(j) = exp(real(CA(j))) * CX( cos(imag(CA(j))), sin(imag(CA(j))) );
      PC(j) = polar_exp(CA(j));
   }
   PrintMaxRelDiff(CC,CD); PrintMaxRelDiff(CC,PC); AssertIsValid(PC);

   cout << "check exponential for imaginary" << endl;
   for (j = 0; j < n; j++)
   {
      CC(j) = exp(IA(j));
      CD(j) = exp(CX(IA(j)));
      PC(j) = polar_exp(IA(j));
   }
   PrintMaxRelDiff(CC,CD); PrintMaxRelDiff(CC,PC); AssertIsValid(PC);

   cout << "check log for complex and polar" << endl;
   for (j = 0; j < n; j++)
   {
      if (CA(j) != 0)
      {
         CC(j) = log(CA(j));
         if (fabs(imag(CC(j))) > pi) cout << "log outlier" << endl;
         CE(j) = log(PA(j));
         if (fabs(imag(CE(j))) > pi) cout << "log outlier" << endl;
         CD(j) = exp(CC(j));
      }
      else { CC(j) = 0; CE(j) = 0; CD(j) = 0; }
   }
   PrintMaxDiff(CC,CE); PrintMaxRelDiff(CA,CD);

   cout << "check log for imaginary" << endl;
   for (j = 0; j < n; j++)
   {
      if (IA(j) != 0)
      {
         CC(j) = log(IA(j));
         if (fabs(imag(CC(j))) > pi) cout << "log outlier" << endl;
         CD(j) = exp(CC(j));
      }
      else { CC(j) = 0; CD(j) = 0; }
   }
   PrintMaxRelDiff(IA,CD);

   cout << "check sin, cos, tan for complex" << endl;
   for (j = 0; j < n; j++)
   {
      CB(j) = exp(_I_ * CA(j));
      CC(j) = sin(CA(j)); CD(j) = cos(CA(j)); CE(j) = tan(CA(j));
      CF(j) = (CB(j) - 1 / CB(j)) / (2 * _I_);
      CG(j) = (CB(j) + 1 / CB(j)) / 2;
      CH(j) = CF(j) / CG(j);
   }
   PrintMaxRelDiff(CC,CF); PrintMaxRelDiff(CD,CG);  PrintMaxRelDiff(CE,CH);

   cout << "check sin, cos, tan for imaginary" << endl;
   for (j = 0; j < n; j++)
   {
      RE(j) = exp(_I_ * IA(j));
      IE(j) = sin(IA(j)); RF(j) = cos(IA(j)); IF(j) = tan(IA(j));
      CF(j) = (RE(j) - 1 / RE(j)) / (2 * _I_);
      CG(j) = (RE(j) + 1 / RE(j)) / 2;
      CH(j) = CF(j) / CG(j);
   }
   PrintMaxRelDiff(IE,CF); PrintMaxRelDiff(RF,CG);  PrintMaxRelDiff(IF,CH);

   cout << "check sinh, cosh, tanh for complex" << endl;
   for (j = 0; j < n; j++)
   {
      CB(j) = exp(CA(j));
      CC(j) = sinh(CA(j)); CD(j) = cosh(CA(j)); CE(j) = tanh(CA(j));
      CF(j) = (CB(j) - 1 / CB(j)) / 2;
      CG(j) = (CB(j) + 1 / CB(j)) / 2;
      CH(j) = CF(j) / CG(j);
   }
   PrintMaxRelDiff(CC,CF); PrintMaxRelDiff(CD,CG);  PrintMaxRelDiff(CE,CH);

   cout << "check sinh, cosh, tanh for imaginary" << endl;
   for (j = 0; j < n; j++)
   {
      CB(j) = exp(IA(j));
      IE(j) = sinh(IA(j)); RE(j) = cosh(IA(j)); IF(j) = tanh(IA(j));
      CF(j) = (CB(j) - 1 / CB(j)) / 2;
      CG(j) = (CB(j) + 1 / CB(j)) / 2;
      CH(j) = CF(j) / CG(j);
   }
   PrintMaxRelDiff(IE,CF); PrintMaxRelDiff(RE,CG);  PrintMaxRelDiff(IF,CH);

   cout << "check sqrt for complex and polar" << endl;
   for (j = 0; j < n; j++)
   {
      CB(j) = sqrt(CA(j));
      if (real(CB(j)) < 0) cout << "CX sqrt outlier" << endl;
      CC(j) = CB(j) * CB(j);
      PB(j) = sqrt(PA(j));
      if (arg(PB(j)) > pi_over_2 && arg(PB(j)) < 3 * pi_over_2)
         cout << "Polar sqrt outlier" << endl;
      PC(j) = PB(j) * PB(j);
   }
   PrintMaxRelDiff(CC,CA);  PrintMaxRelDiff(PC,PA);
   AssertIsValid(PB); AssertIsValid(PC);

   cout << "check sqrt for imag" << endl;
   for (j = 0; j < n; j++)
   {
      CB(j) = sqrt(IA(j));
      if (real(CB(j)) < 0) cout << "sqrt outlier" << endl;
      CC(j) = CB(j) * CB(j);
   }
   PrintMaxRelDiff(CC,IA);

   cout << "check square for complex and polar" << endl;
   for (j = 0; j < n; j++)
   {
      CB(j) = square(CA(j));
      CC(j) = CA(j) * CA(j);
      PB(j) = square(PA(j));
      PC(j) = PA(j) * PA(j);
   }
   PrintMaxRelDiff(CC,CB);  PrintMaxRelDiff(PC,PB);  PrintMaxRelDiff(CB,PB);
   AssertIsValid(PB); AssertIsValid(PC);

   cout << "check sum_square for complex and polar" << endl;
   for (j = 0; j < n; j++)
   {
      RE(j) = sum_square(CA(j));
      CC(j) = CA(j) * conj(CA(j));
      RG(j) = sum_square(PA(j));
      PC(j) = PA(j) * conj(PA(j));
   }
   PrintMaxRelDiff(CC,RE);  PrintMaxRelDiff(PC,RG);  PrintMaxRelDiff(RE,RG);

   cout << "check square and sum_square for imaginary" << endl;
   for (j = 0; j < n; j++)
   {
      RE(j) = square(IA(j));
      RF(j) = IA(j) * IA(j);
      RG(j) = -sum_square(IA(j));
   }
   PrintMaxRelDiff(RF,RE); PrintMaxRelDiff(RF,RG);

   cout << "check operator+ and operator- for complex and polar" << endl;
   for (j = 0; j < n; j++)
   {
      CB(j) = +CA(j);
      CC(j) = -CA(j);
      CD(j) = 0.0 - CA(j);
      PB(j) = +PA(j);
      PC(j) = -PA(j);
      PD(j) = 0.0 - PA(j);
   }
   PrintMaxRelDiff(CA,CB);  PrintMaxRelDiff(CC,CD);  PrintMaxRelDiff(PC,CC);
   PrintMaxRelDiff(PA,PB);  PrintMaxRelDiff(PC,PD);
   AssertIsValid(PB); AssertIsValid(PC); AssertIsValid(PD);

   cout << "check real, imag, conj for complex and polar" << endl;
   for (j = 0; j < n; j++)
   {
      RE(j) = real(CA(j));
      RF(j) = imag(CA(j));
      CB(j).real() = RE(j);
      CB(j).imag() = RF(j);
      CC(j).real() = RE(j);
      CC(j).imag() = -RF(j);
      CD(j) = conj(CA(j));
      RG(j) = real(PA(j));
      RH(j) = imag(PA(j));
   }
   PrintMaxRelDiff(CA,CB);  PrintMaxRelDiff(CC,CD);
   PrintMaxRelDiff(RE,RG);  PrintMaxRelDiff(RF,RH);

   cout << "check real, imag, conj for imaginary" << endl;
   for (j = 0; j < n; j++)
   {
      RE(j) = real(IA(j));
      RF(j) = imag(IA(j));
      IB(j) = Imag(RF(j));
      IC(j) = -IB(j);
      ID(j) = conj(IA(j));
      RG(j) = 0;
   }
   PrintMaxRelDiff(IA,IB);  PrintMaxRelDiff(IC,ID);
   PrintMaxRelDiff(RE,RG);

   cout << "check arg for complex and polar" << endl;
   for (j = 0; j < n; j++)
   {
      if (CA(j) != 0.0)
      {
         RE(j) = arg(CA(j));
         RF(j) = atan2(imag(CA(j)), real(CA(j)));
         RG(j) = arg(PA(j));
         if (fabs(RE(j)) > pi) cout << "outlier error in arg" << endl;
      }
      else { RE(j) = 0; RF(j) = 0; RG(j) = 0; }
   }
   PrintMaxRelDiff(RE,RF);  PrintMaxRelDiff(RF,RG);

   cout << "check arg for imaginary" << endl;
   for (j = 0; j < n; j++)
   {
      if (IA(j) != 0.0)
      {
         RE(j) = arg(IA(j));
         RF(j) = atan2(imag(IA(j)), 0.0);
         if (fabs(RE(j)) > pi) cout << "outlier error in arg" << endl;
      }
      else { RE(j) = 0; RF(j) = 0; }
   }
   PrintMaxRelDiff(RE,RF);

   cout << "check fabs and cabs for complex and polar" << endl;
   for (j = 0; j < n; j++)
   {
      RE(j) = fabs(CA(j));
      RF(j) = cabs(CA(j));
      RG(j) = fabs(PA(j));
      RH(j) = cabs(PA(j));
      RD(j) = sqrt(CA(j) * CA(j).conj()).real();
   }
   PrintMaxDiff(RE,RF);  PrintMaxDiff(RF,RG);
   PrintMaxDiff(RG,RH);  PrintMaxDiff(RE,RD);

   cout << "check fabs and cabs for imaginary" << endl;
   for (j = 0; j < n; j++)
   {
      RE(j) = fabs(IA(j));
      RF(j) = cabs(IA(j));
      RD(j) = fabs(IA(j).imag());
   }
   PrintMaxDiff(RE,RF);  PrintMaxDiff(RE,RD);

   cout << "check norm1 for complex" << endl;
   for (j = 0; j < n; j++)
   {
      RE(j) = norm1(CA(j));
      RF(j) = fabs(RA(j)) + fabs(IA(j));
   }
   PrintMaxDiff(RE,RF);

   cout << "check polar constructor" << endl;
   for (j = 0; j < n; j++)
   {
      PB(j) = Polar(RA(j), RB(j));
      CB(j) = RA(j) * exp(_I_ * RB(j));
      PC(j) = Polar(RA(j), 10.0 * RB(j));
      CC(j) = RA(j) * exp(10.0 * _I_ * RB(j));
   }
   AssertIsValid(PB); PrintMaxRelDiff(PB, CB);
   AssertIsValid(PC); PrintMaxRelDiff(PC, CC);

   cout << "check polar variables access" << endl;
   for (j = 0; j < n; j++)
   {
      PC(j).theta() = PA(j).theta();
      PC(j).quadrant() = PA(j).quadrant();
      PC(j).cabs() = PA(j).cabs();
   }
   AssertIsValid(PC); PrintMaxDiff(PA, PC);

   cout << "check imag_as_imag" << endl;
   for (j = 0; j < n; j++)
   {
      IE(j) = CA(j).imag_as_imag();
      IF(j) = IA(j).imag_as_imag();
      IG(j) = PA(j).imag_as_imag();
   }
   PrintMaxDiff(IA, IE); PrintMaxDiff(IA, IF); PrintMaxDiff(IA, IG);

   cout << "check imag_as_imag global function" << endl;
   for (j = 0; j < n; j++)
   {
      IE(j) = imag_as_imag(CA(j));
      IF(j) = imag_as_imag(IA(j));
      IG(j) = imag_as_imag(PA(j));
   }
   PrintMaxDiff(IA, IE); PrintMaxDiff(IA, IF); PrintMaxDiff(IA, IG);

   cout << "end unary tests" << endl;

   return OverallMax;
}