//This computes Pade Approximants to sqr((1+2x)/(1+x)) and evaluates them
//at infinity (see Baker,Essentials of Pade Approximants, p.5). The number
//of significant figures (32) and maximum order of approximation (19) could
//be increased easily.

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

const int precision = 8;    //note const           
const int zmaxsize=40;
const int cmaxsize=0;

//#include "genprecd.h"    //old style
//#include "seriesd.h"          
//#include "genlind.h"
//#include "paded.h"

#include "genprec.h"      //new style
#include "cgenprec.h"
#include "series&p.h"

int n;
int nmax;

realscalar p_ram;unsigned int *p=p_ram.s;
realscalar xuv_ram;unsigned int *xuv=xuv_ram.s;
realvector a_ram;unsigned int * *a=a_ram.v;
realvector b_ram;unsigned int * *b=b_ram.v;
realvector c_ram;unsigned int * *c=c_ram.v;

int main()
  {int i;
   zsupe();                 //run this if e is needed
   zsuplogsp();             //run this if log base e of 10 is needed
   n=38;
   zconinttomp(1);       //compute series for sqrt[(1+2x)/(1+x)]
   zstore(a[0]);         //a=1 + 2x,b+1+x
   zstore(b[0]);
   zstore(b[1]);
   zconinttomp(2);
   zstore(a[1]);
   for (i=2;i<n+1;i++)       //remaining terms in a,b are zero
     {zzero(a[i]);
      zzero(b[i]);}            //end set up
   zconinttomp(1);
   zchs(z);
   zstore(p);
   stonth(n,b,p);            //b^-1
   sstorit(n,sout,c);
   smult(n,a,c);          //a*b^-1
   sstorit(n,sout,c);             //=c
   zconinttomp(1);
   zdivint(z,2);
   zstore(p);
   stonth(n,c,p);           //c^(1/2)
   sstorit(n,sout,c);
   for (i=1;i<n/2;i++)
     {cout<<i<<"th Pade Approximant at infinity"<<'\n';
      padeapprox(i,i,c);                                     
      zdiv(pnum[i],pden[i]);
      zstore(xuv);
      zprnt(xuv);
      zmessage=" ";
      zbreak();
     }
   return 0;
  }