/* FILE: myclx4.c begins                               */
/*  gcc myclx4.c -L../lib -o a -lm -lnr  */
   
#define NRANSI
#include <stdio.h>
#include <math.h>
#include "../lib/nr.h"
#include "../lib/nrutil.h"
#include "../lib/nrutil.c"
#include "../lib/complex.h"  


#include "../util/matutl6.c"
#define NP 20
#define MP 20
#define MAXSTR 80

int main()
{
  int j;
  long int ma,na;
  float **a,**b;
  fcomplex z, z1, w, w1;
  char pmdm;    
  ma=20; na=4; 
  printf("Test program for complex arithmetic. Enter 1,2,3 or 4:\n");
  printf("1: Addition/subtraction \n");
  printf("2: Multiplication/division \n");
  printf("3: Root\n");
  printf("4: Break\n");
  scanf("%s",&pmdm );
  switch (pmdm) {
  case '1': {printf("Case 1: \n");  break;};
  case '2': {printf("Case 2: \n"); break;};
  case '3': {printf("Case 3: \n"); na=2;  break;};
  case '4' : return 0;}
  a=ranmat2(ma,na,-0.98,0.99);
  b=matrix(1,ma,1,na+3);
  printf("\nThe matrix you generated is: \n");
  showmat(a,ma,na);
  getchar();
  printf("%15s\n","Error estimate");
  for (j=1;j<=ma;j++) {
    z.r=a[j][1]; z.i=a[j][2];
    b[j][1]=a[j][1]; b[j][2]=a[j][2];
    switch (pmdm) {
    case '1': {
      w.r=a[j][3]; w.i=a[j][4];
                 b[j][3]=a[j][3]; b[j][4]=a[j][4];
                 z1=Cadd(z,w);
                 w1=Csub(z1,w); w1=Csub(w1,z) ; /* w1=(z+w)-w-z */
                 break; };
    case '2': {
      w.r=a[j][3]; w.i=a[j][4];
      b[j][3]=a[j][3]; b[j][4]=a[j][4];
      z1=Cmul(z, w);
      w1=Cdiv(z1,w); w1=Csub(w1,z); /* w1=(z*w)/w-z */
      break; };
    case '3': {    
      z1=Csqrt(z);
      w=Cmul(z1,z1);
      w1=Csub(w,z);  /* w1=sqrt(z)^2-z */
      break; };
    default : break;}
    b[j][na+1]=z1.r; b[j][na+2] =z1.i;
    b[j][na+3]=Cabs(w1); 
    printf("%15.4e \n",b[j][na+3]);
  }
  switch (pmdm) {
  case '1': 
    {printf("\n%10s %10s %10s %10s %10s %10s %10s\n", "Re{z}",
	    "Im{z}", "Re{w}","Im{w}", "Re{z+w}","Im{z+w}","Error");
    break;}
  case '2': 
    {printf("\n%10s %10s %10s %10s %10s %10s %10s\n", "Re{z}",
	    "Im{z}", "Re{w}","Im{w}", "Re{z*w}","Im{z*w}","Error");
    break;}
  case '3':
    {printf("\n%10s %10s %12s %12s %10s\n", "Re{z}",
	    "Im{z}", "Re{sqrt(z)}","Im{sqrt(z)}" ,"Error"); 
    break;}
  default: break; 
  }
  showmat(b,ma,na+3);
  for (j=1;j<=na+3;j++) printf("%12.4e \n", b[j][na+3] );
  switch (pmdm) {
  case '1' : case '2' : break;
  case '3': 
             {printf("\nNext z= exp(i theta), theta=2*pi*(j-1)/20,
j=1,2,...,20 \n");
	     printf("\n%10s %10s %12s %12s  \n", "Re{z}",
		    "Im{z}", "Re{sqrt(z)}","Im{sqrt(z)}" ); 
	     for (j=1;j<=20;j++) {
	       z.r=cos(2*M_PI*(j-1)/20);
	       z.i=sin(2*M_PI*(j-1)/20);
	       w=Csqrt(z);
	       printf("%11.5f %11.5f %11.5f %11.5f \n", 
		      z.r, z.i, w.r, w.i);} 
             break;}
  default: break;
  }
  free_matrix(a,1,ma,1,na);
  free_matrix(b,1,ma,1,na+3);
  return 0;
}
#undef NRANSI
/* FILE: myclx4.c ends                                */