#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const int precision=8;
const int zmaxsize=3;
const int cmaxsize=40;

#include "genprec.h"
#include "cgenprec.h"
#include "apps.h"

int n;
int nmax;

complexmp p,cwz; 
complexmp xmad;           
complexmp a[cmaxsizepad];
complexmp b[cmaxsizepad];
complexmp c[cmaxsizepad];
cpademp d;
complexmp dnum[cmaxsizepad],dden[cmaxsizepad];
complexmp myrootsofden[cmaxsizepad];
complexmp myrootsofnum[cmaxsizepad];

int main()
  {int i;
   zsuppi();
   zsupe();                 //run this if e is needed
   zsuplogsp();             //run this if log base e of 10 is needed
   n=20;       
              
  cout<<"The example is the 10,10 Pade Approximant to"<<'\n'; 
  cout<<"sqrt(1+(2+i)x)/sqrt(1+(1+i)x)"<<'\n';
  cout<<"The zeros of the numerator and denominator are calculated, and the"<<'\n';
  cout<<"polynomial solver prints out rrr... as it finds the roots"<<'\n';
  cout<<"at the end the residues at the poles of the denominator are printed"<<'\n';
  zmessage="The calculation takes a few minutes on a 486DX2 50mhz PC";
  zbreak();
  cout<<'\n';

   a[0]=cinttomp(1,0);       //compute series 
   a[1]=cinttomp(2,1);
   for (i=2;i<n+1;i++)       //remaining terms in a are zero
     {a[i]=cinttomp(0,0);}            
   p=cinttomp(1,0);
   p=p/2;
   a<cstonth(n,a,p);            //sqrt(a)

   for (i=0;i<n+1;i++) printmp(a[i]);
   zmessage="coeffs in sqrt(1+(2+i)x)";
   zbreak();

   b[0]=cinttomp(1,0);         //b=1+(1+i)x
   b[1]=cinttomp(1,1);
   for (i=2;i<n+1;i++)
     {b[i]=cinttomp(0,0);}
   p=-p;
   b<cstonth(n,b,p);            //1/sqrt(b)

   for (i=0;i<n+1;i++) printmp(b[i]);
   zmessage="coeffs in 1/sqrt(1+(1+i)x)";
   zbreak();

   c<csmult(n,a,b);
   for (i=0;i<n+1;i++) printmp(c[i]);
   zmessage="the series for sqrt(1+(2+i)x)/sqrt(1+(1+i)x)";
   zbreak();

   d<cpadeapprox(n/2,n/2,c);     //form pade approximant
   dnum<cextractnum(n/2,d);
   dden<cextractden(n/2,d);

   for (i=0;i<n/2+1;i++)         //print numerator
     {printmp(dnum[i]);}
   zmessage="numerator";
   zbreak();
   cout<<'\n';

   for (i=0;i<n/2+1;i++)        //print denominator
     {printmp(dden[i]);}
   zmessage="denominator";
   zbreak();
   cout<<'\n';

   rotationangle=inttomp(2)/10;   //real polys.txt

   myrootsofnum<csolvpoly(n/2,dnum,rotationangle);  //get roots of 
   for (i=1;i<n/2+1;i++)                            //numerator
     {printmp(myrootsofnum[i]);}
   zmessage="roots of numerator";
   zbreak();
   cout<<'\n';
   myrootsofden<csolvpoly(n/2,dden,rotationangle);   //and denominator
   for (i=1;i<n/2+1;i++)
     {printmp(myrootsofden[i]);}
   zmessage="roots of denominator";
   zbreak();
   dden<csder(n/2,dden);
   for (i=0;i<n/2-1;i++) printmp(dden[i]);
   zmessage="terms in derivative of denominator";
   zbreak();
   for (i=1;i<n/2+1;i++)
     {xmad<cseval(n/2-1,dden,myrootsofden[i]);
      printmp(xmad);
      zmessage="derivative of denominator at ith root of denominator, ";
      cerr<< zmessage << "i="<<i<<'\n'; 
      cwz<cseval(n/2,dnum,myrootsofden[i]);
      printmp(cwz);
      zmessage="numerator at ith root of denominator";
      cerr<<zmessage<<'\n';   //zbreak();
      cwz<cwz/xmad;
      printmp(cwz);
      zmessage="residue=N/(dD/dx) at ith root of denominator";
      zbreak();
     }

   return 0;
  }