cxtest2.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 array2<T1>& A, const array2<T2>& B)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   Real maxdiff = 0.0;
   for (int i = 0; i < nr; i++) for (int j = 0; j < nc; j++)
   {
      Real diff = fabs(A(i,j) - B(i,j));
      if (maxdiff < diff) maxdiff = diff;
      if (diff > 0.0000001)
      {
         cout << i << " " << j << ": ";
         cout << real(A(i,j)) << " " << imag(A(i,j)) << "; ";
         cout << real(B(i,j)) << " " << imag(B(i,j)) << " -> " << diff << endl;
      }
   }
   cout << maxdiff << endl;
   if (OverallMax < maxdiff) OverallMax = maxdiff;
}

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

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

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

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

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

template <class T1>
void operator+=(array2<T1>& A, ImaginaryUnit)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   for (int i = 0; i < nr; i++) for (int j = 0; j < nc; j++) A(i,j) += _I_;
}

template <class T1>
void operator-=(array2<T1>& A, ImaginaryUnit)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   for (int i = 0; i < nr; i++) for (int j = 0; j < nc; j++) A(i,j) -= _I_;
}

template <class T1>
void operator*=(array2<T1>& A, ImaginaryUnit)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   for (int i = 0; i < nr; i++) for (int j = 0; j < nc; j++) A(i,j) *= _I_;
}

template <class T1>
void operator/=(array2<T1>& A, ImaginaryUnit)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   for (int i = 0; i < nr; i++) for (int j = 0; j < nc; j++) A(i,j) /= _I_;
}

void AssertIsValid(const array2<Polar>& p)
{
   int nr = p.N_rows(); int nc = p.N_cols();
   for (int i = 0; i < nr; i++) for (int j = 0; j < nc; j++)
      p(i,j).AssertIsValid();
}


template <class T1, class T2>
void TestCompare(const array2<T1>& A, const array2<T2>& B)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   bool OnDiag = true; bool OffDiag = false;
   int eq = 0; int neq = 0;

   for (int i = 0; i < nr; i++) for (int j = 0; j < nc; j++)
   {
      CX a(A(i,j)), b(B(i,j));
      if (a == b)                 // complex versions equal
      {
         if ( A(i,j) == B(i,j) ) {} else OnDiag = false;
         if ( A(i,j) != B(i,j) ) OnDiag = false;
         eq++;
      }
      else
      {
         if ( A(i,j) != B(i,j) ) {} else OffDiag = true;
         if ( A(i,j) == B(i,j) ) OffDiag = true;
         neq++;
      }
   }

   if (OnDiag && !OffDiag) cout << "test passes" << endl;
   else { cout << "test fails" << endl; OverallMax = 100; }
   cout << "eq = " << eq << ";  neq = " << neq << endl;
}

template <class T1>
void TestCompareWithI(const array2<T1>& A)
{
   int nr = A.N_rows(); int nc = A.N_cols();
   bool OnDiag = true; bool OffDiag = false;
   int eq = 0; int neq = 0;

   for (int i = 0; i < nr; i++) for (int j = 0; j < nc; j++)
   {
      CX a(A(i,j));
      if (a == CX(0,1))                 // complex versions equal
      {
         if ( A(i,j) == _I_ && _I_ == A(i,j) ) {} else OnDiag = false;
         if ( A(i,j) != _I_ || _I_ != A(i,j) ) OnDiag = false;
         eq++;
      }
      else
      {
         if ( A(i,j) != _I_ && _I_ != A(i,j) ) {} else OffDiag = true;
         if ( A(i,j) == _I_ || _I_ == A(i,j) ) OffDiag = true;
         neq++;
      }
   }

   if (OnDiag && !OffDiag) cout << "test passes" << endl;
   else { cout << "test fails" << endl; OverallMax = 100; }
   cout << "eq = " << eq << ";  neq = " << neq << endl;
}
Real cxtest2()
{
   cout << "begin binary tests" << endl;

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

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

   for (i = 0; i < n; i++)
   {
      int j, k;
      CA(0, i) = 0.0;
      CA(1, i) = 1.0;
      CA(2, i) = -1.0;
      CA(3, i) = 1.9;
      CA(4, i) = -1.9;
      CA(5, i) = 0.75;
      CA(6, i) = -0.75;
      CA(7, i) = 2.35;
      CA(8, i) = -2.35;
      CA(9, i) = 4.33;
      CA(10, i) = -4.33;
      for (j = 1; j <= 10; j++) CA(j + 10, i) = Imag(CA(j, i).real());
      for (j = 1; j <= 10; j++) for (k = 1; k <= 10; k++)
         CA(j + 10 * k + 10, i) = CA(j, i) + CA(k + 10, i);
   }

   for (i = 0; i < n; i++)
   { for (int j = 0; j < n; j++) { CB(i, j) = CA(j, i); }  }

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

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      RA(i,j) = real(CA(i,j)); IA(i,j) = Imag(imag(CA(i,j))); PA(i,j) = CA(i,j);
      RB(i,j) = real(CB(i,j)); IB(i,j) = Imag(imag(CB(i,j))); PB(i,j) = CB(i,j);
      RC(i,j) = imag(CA(i,j)); RD(i,j) = imag(CB(i,j));
      PA(i,j).AssertIsValid(); PB(i,j).AssertIsValid();
   }

   PrintMaxDiff(PA, CA);  PrintMaxDiff(PB, CB);

   cout << "check complex op complex" << endl;

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) + CB(i,j);
      CD(i,j) = CX(RA(i,j)+RB(i,j), RC(i,j)+RD(i,j));
   }
   PrintMaxDiff(CC,CD);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) - CB(i,j);
      CD(i,j) = CX(RA(i,j)-RB(i,j), RC(i,j)-RD(i,j));
   }
   PrintMaxDiff(CC,CD);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) * CB(i,j);
      CD(i,j) = CX(RA(i,j)*RB(i,j)-RC(i,j)*RD(i,j),
                   RA(i,j)*RD(i,j)+RB(i,j)*RC(i,j));
   }
   PrintMaxDiff(CC,CD);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      if (j != 0)
      {
         Real f = cabs(CB(i,j)); f = f * f;
         CC(i,j) = CA(i,j) / CB(i,j);
         CD(i,j) = CX((RA(i,j)*RB(i,j)+RC(i,j)*RD(i,j))/f,
                      (RB(i,j)*RC(i,j)-RA(i,j)*RD(i,j))/f);
      }
      else { CC(i,j) = 0.0; CD(i,j) = 0.0; }
   }
   PrintMaxDiff(CC,CD);

   cout << "check complex op real, real op complex" << endl;

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) + RB(i,j);
      CD(i,j) = RB(i,j) + CA(i,j);
      CE(i,j) = CA(i,j) + CX(RB(i,j));
   }
   PrintMaxDiff(CC,CE);  PrintMaxDiff(CD,CE);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) - RB(i,j);
      CD(i,j) = RB(i,j) - CA(i,j);
      CE(i,j) = CA(i,j) - CX(RB(i,j));
      CF(i,j) = CX(RB(i,j)) - CA(i,j);
   }
   PrintMaxDiff(CC,CE);  PrintMaxDiff(CD,CF);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) * RB(i,j);
      CD(i,j) = RB(i,j) * CA(i,j);
      CE(i,j) = CA(i,j) * CX(RB(i,j));
   }
   PrintMaxDiff(CC,CE);  PrintMaxDiff(CD,CE);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      if (RB(i,j) != 0.0)
      {
         CC(i,j) = CA(i,j) / RB(i,j);
         CE(i,j) = CA(i,j) / CX(RB(i,j));
      }
      else { CC(i,j) = 0; CE(i,j) = 0; }
      if (CA(i,j) != 0.0)
      {
         CF(i,j) = CX(RB(i,j)) / CA(i,j);
         CD(i,j) = RB(i,j) / CA(i,j);
      }
      else { CD(i,j) = 0; CF(i,j) = 0; }
   }
   PrintMaxDiff(CC,CE);  PrintMaxDiff(CD,CF);

   cout << "check complex op imag, imag op complex" << endl;

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) + IB(i,j);
      CD(i,j) = IB(i,j) + CA(i,j);
      CE(i,j) = CA(i,j) + CX(IB(i,j));
   }
   PrintMaxDiff(CC,CE);  PrintMaxDiff(CD,CE);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) - IB(i,j);
      CD(i,j) = IB(i,j) - CA(i,j);
      CE(i,j) = CA(i,j) - CX(IB(i,j));
      CF(i,j) = CX(IB(i,j)) - CA(i,j);
   }
   PrintMaxDiff(CC,CE);  PrintMaxDiff(CD,CF);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) * IB(i,j);
      CD(i,j) = IB(i,j) * CA(i,j);
      CE(i,j) = CA(i,j) * CX(IB(i,j));
   }
   PrintMaxDiff(CC,CE);  PrintMaxDiff(CD,CE);

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      if (IB(i,j) != 0.0)
      {
         CC(i,j) = CA(i,j) / IB(i,j);
         CE(i,j) = CA(i,j) / CX(IB(i,j));
      }
      else { CC(i,j) = 0; CE(i,j) = 0; }
      if (CA(i,j) != 0.0)
      {
         CF(i,j) = CX(IB(i,j)) / CA(i,j);
         CD(i,j) = IB(i,j) / CA(i,j);
      }
      else { CD(i,j) = 0; CF(i,j) = 0; }
   }
   PrintMaxDiff(CC,CE);  PrintMaxDiff(CD,CF);

   cout << "check complex op _I_, _I_ op complex" << endl;

   for (i = 0; i < n; i++) for (j = 0; j < n; j++)
   {
      CC(i,j) = CA(i,j) + _I_;
      CD(i,j) = _I_ + CA(i,j);