#define WANT_MATH#define WANT_STREAM#include "include.h"#include "myexcept.h"#include "array1.h"#include "cx.h"#include "cxtest.h"#ifdef use_namespaceusing namespace RBD_COMMON;             // needed for VC++using namespace RBD_COMPLEX;#endifstatic 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);