/* test_complex.c    using complex.h and complex.c */

#include "complex.h"
#include <stdio.h>

int main(int argc, char *argv[])
{
  complex z,zc,zn,znc;
  double x, y, nrm;
  double mag, ang;
  complex t1, t2, t3,  adet, aadet;
  complex A[4][4], X[4], Y[4];
  complex AA[4][4], XX[4], TT[4], YY[4];
  complex BB[4][4], CC[4][4];
  int i, j, n;

  printf("test_complex.c \n");
  z = cmplx(1.5,-3.5);
  zc = cxconj(z);
  zn = cxneg(z);
  znc = cxconj(cxneg(z));
  printf("z=(%g,%g) \n", z.re, z.im);
  printf("z=(%g,%g) \n", z);         /* complex is two doubles */

  printf("real(z) 1.5 =%g \n",        real(z));
  printf("imag(z) -3.5 =%g \n",        imag(z));
  printf("cxconj(z) (1.5,3.5) =(%g,%g) \n", zc);
  printf("cxneg(z) (-1.5,3.5) =(%g,%g) \n", zn);
  printf("cxconj(cxneg(z)) (-1.5,-3.5) =(%g,%g) \n", znc);
  printf("cxabs(z) =%g \n",       cxabs(z));
  printf("cxarg(z) =%g \n",       cxarg(z));

  printf("cxadd(z,z) =(%g.%g) \n", cxadd(z,z));    
  printf("cxsub(z,z/2.0) =(%g,%g) \n", cxsub(z,cdivd(z,2.0)));  
  printf("cxmul(z,z) = (%g,%g) \n", cxmul(z,z));      
  printf("cxdiv(z*z,z) = (%g,%g) \n", cxdiv(cxmul(z,z),z));  
  printf("\n");
  printf("cxsqrt(z*z) = (%g,%g) \n", cxsqrt(cxmul(z,z)));
  printf("cxsqrt(%g,%g)=(%g,%g), squared=(%g,%g)\n",
         z,cxsqrt(z),cxmul(cxsqrt(z),cxsqrt(z)));
  printf("cxsqrt(%g,%g)=(%g,%g), squared=(%g,%g)\n",
         zc,cxsqrt(zc),cxmul(cxsqrt(zc),cxsqrt(zc)));
  printf("cxsqrt(%g,%g)=(%g,%g), squared=(%g,%g)\n",
         zn,cxsqrt(zn),cxmul(cxsqrt(zn),cxsqrt(zn)));
  printf("cxsqrt(%g,%g)=(%g,%g), squared=(%g,%g)\n",
         znc,cxsqrt(znc),cxmul(cxsqrt(znc),cxsqrt(znc)));
  printf("\n");
  t1 = cxexp(z);
  printf("cxexp(z)   =  (%g,%g)\n", t1);
  t2 = cxlog(z);
  printf("cxlog(z)   =  (%g,%g)\n", t2);
  printf("z=(%g,%g), cxlog(cxexp(z))   =  (%g,%g)\n", z, cxlog(t1));
  printf("z=(%g,%g), cxexp(cxlog(z))   =  (%g,%g)\n", z, cxexp(t2));
  t1 = cxexp(zc);
  printf("cxexp(z)   =  (%g,%g)\n", t1);
  t2 = cxlog(zc);
  printf("cxlog(z)   =  (%g,%g)\n", t2);
  printf("z=(%g,%g), cxlog(cxexp(z))   =  (%g,%g)\n", zc, cxlog(t1));
  printf("z=(%g,%g), cxexp(cxlog(z))   =  (%g,%g)\n", zc, cxexp(t2));
  t1 = cxexp(zn);
  printf("cxexp(z)   =  (%g,%g)\n", t1);
  t2 = cxlog(zn);
  printf("cxlog(z)   =  (%g,%g)\n", t2);
  printf("z=(%g,%g), cxlog(cxexp(z))   =  (%g,%g)\n", zn, cxlog(t1));
  printf("z=(%g,%g), cxexp(cxlog(z))   =  (%g,%g)\n", zn, cxexp(t2));
  t1 = cxexp(znc);
  printf("cxexp(z)   =  (%g,%g)\n", t1);
  t2 = cxlog(znc);
  printf("cxlog(z)   =  (%g,%g)\n", t2);
  printf("z=(%g,%g), cxlog(cxexp(z))   =  (%g,%g)\n", znc, cxlog(t1));
  printf("z=(%g,%g), cxexp(cxlog(z))   =  (%g,%g)\n", znc, cxexp(t2));
  
  printf("\n");
  printf("cxsin(z)   =  (%g,%g)\n", cxsin(z));   
  printf("cxcos(z)   =  (%g,%g)\n", cxcos(z));   
  printf("cxtan(z)   =  (%g,%g)\n", cxtan(z));   
  printf("cxatan(z)  =  (%g,%g)\n", cxatan(z));   
  printf("cxsinh(z)  =  (%g,%g)\n", cxsinh(z));  
  printf("cxcosh(z)  =  (%g,%g)\n", cxcosh(z));  
  printf("cxtanh(z)  =  (%g,%g)\n", cxtanh(z));  
  printf("cxatanh(z) =  (%g,%g)\n", cxatanh(z));  
  printf("\n");

  n = 4;
  cxmran(n,A);
  printf("random matrix A\n");
  cxmprint(n,A);
  printf("\n");
  cxvran(n,X);
  printf("random vector X\n");
  cxvprint(n,X);
  cxmvmul(n,A,X,Y);
  printf("\n");
  printf("Y = A*X  cxmvmul(n,A,X,Y) \n");
  cxvprint(n,Y);
  printf("\n");

  cxmcopy(n,A,AA); /* save A */
  printf("cxsimeq(n,A,Y,XX) solve for XX=X\n");
  cxsimeq(n, AA, Y, XX);
  printf("check computed solution XX\n");
  cxvprint(n,XX);
  cxvsub(n,XX,X,TT);
  printf("TT=XX-X should be near zero \n");
  cxvprint(n,TT);
  nrm = cxvnorm(n,TT);
  printf("norm of all zero vector near zero %g\n",nrm);
  printf("\n");
  printf("YY=A*X  cxmvmul(n,A,XX,YY) \n");
  cxmvmul(n,A,XX,YY);
  printf("YY is\n");
  cxvprint(n,YY);
  printf("check computed solution YY-Y should be near zero\n");
  cxvsub(n,YY,Y,TT);
  cxvprint(n,TT);
  nrm = cxvnorm(n,TT);
  printf("norm of all zero vector near zero %g\n",nrm);
  printf("\n");

  printf("cxinv(n,A,AA) is AA\n");
  cxinv(n,A,AA);
  cxmprint(n,AA);
  printf("\n");
  cxmmul(n,A,AA,BB);
  printf("should be identity matrix BB\n");
  cxmprint(n,BB);
  printf("\n");  
  cxmident(n,CC);
  cxmsub(n,BB,CC,AA);
  printf("should be all zero matrix inv(A)*A-I\n");
  cxmprint(n,AA);
  nrm = cxmnorm(n,AA);
  printf("norm of all zero matrix near zero %g\n",nrm);
  printf("\n");  

  adet = cxdet(n,A);
  printf("cxdet(n,A)=(%g,%g)\n",adet);
  cxinv(n,A,AA);
  aadet = cxdet(n,AA);
  printf("cxdet(inv(A))=(%g,%g)  and product (%g,%g)\n",
         aadet, cxmul(adet,aadet));
  printf("\n");  

  printf("end test_complex.c\n");
  return 0;
} /* end test_complex.c */