(*PADE' RATIONAL APPROXIMATION ALGORITHM 8.1
*
*  To obtain the rational approximation
*    r(x) = p(x) / q(x)
*      	  = (p0 + p1*x + ... + pn*x^n)/(q0 + q1*x + ... + qm*x^m)
*   for a given function f(x):
*
*   INPUT: nonnegative integers m and n.
* 
*  OUTPUT: coefficients qo, q1, ... , qm, p0, p1, ... , pn.
*
*  The coefficients of the Maclaurin poplynomial a0, a1, ... could
*   be calculated instead of input as is assumed in this program.
*)
OK=0;
While[OK == 0,
   LM=Input["This is Pade' Approximation.\n
   \n
   Input m\n"];
   LN=Input["Input n\n"];
   BN=LM+LN;
   If[LM >= 0 && LN >= 0,
      OK = 1,
      Input["m and n must both be nonnegative\n
      \n
      Press 1 [enter] to continue\n"];
   ];
   If[LM == 0 && LN == 0,
      OK = 0;
      Input["Not both m and n cannot be zero\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
OK = 0;
While[OK == 0,
   FLAG=Input["The McLaurin coefficients a[0], a[1], ... a[n]\n
   are to be input.\n
   Choice of input method:\n
   1. Input entry by entry from keyboard\n
   2. Input data from a text file\n
   Choose 1 or 2 please\n"];
   If[FLAG == 1 || FLAG == 2,
      OK = 1;
   ];
];
If[FLAG == 1,
   For[i=0,
      i <= BN,
      i++,
      AA[i]=Input["Input A["<>ToString[i]<>"]\n"];
   ];
];
If[FLAG == 2,
   AAA=InputString["As many entries as desired can be placed\n
   on each line of the file separated by a blank.\n
   \n
   Has such a text file been created?\n
   Enter 'yes' or 'no'\n"];
   If[AAA=="Y" || AAA=="y" || AAA=="yes",
      NAME=InputString["Input the file name in the form - \n
      drive:\ name.ext      for example\n
         A:\\DATA.DTA\n"];
      INP=OpenRead[NAME];
      For[i = 0,
         i <= BN,
         i++,
         AA[i]=Read[INP,Number];
      ];
      Close[INP],
      Input["Please create the input file.\n
      This program will end so the input file can\n
      be created.\n
      \n
      Press 1 [enter] to continue\n"];
      OK = 0;
   ];
];
If[OK == 1,
   (* Step 1 *)
   n=BN;
   M=n+1;
   (* Step 2 - performed in input *)
   For[i=1,
      i <= n,
      i++,
      NROW[i-1]=i;
   ];
   (* initialize row pointer for linear system *)
   NN=n-1;
   (* Step 3 *)
   q[0]=1.0;
   P[0]=AA[0];
   (* Step 4 - Set up a linear system, but use A[i,j] 
      instead of B[i,j]. *)
   For[i=1,
      i <= n,
      i++,
      (* Step 5 *)
      For[J=1,
         J <= i-1,
	 J++,
   	 If[J <= LN,
	    A[i-1,J-1]=0;
	 ];
      ];
      (* Step 6 *)
      If[i <= LN,
         A[i-1,i-1]=1.0;
      ];
      (* Step 7 *)
      For[J=i+1,
	 J <= n,
	 J++,
 	 A[i-1,J-1]=0;
      ];
      (* Step 8 *)
      For[J = 1,
	 J <= i,
	 J++,
	 If[J <= LM,
	    A[i-1,LN+J-1]=-AA[i-J];
	 ];
      ];
      (* Step 9 *)
      For[J = LN+i+1,
	 J <= n,
	 J++,
	 A[i-1,J-1]=0;
      ];
      (* Step 10 *)
      A[i-1,n]=AA[i];
   ];
   (* Solve the linear system using partial pivoting *)
   i=LN+1;
   (* Step 11 *)
   While[OK == 1 && i <= NN,
      (* Step 12 *)
      IMAX=NROW[i-1];
      AMAX=Abs[A[IMAX-1,i-1]];
      IMAX=i;
      JJ=i+1;
      For[IP = JJ,
	 IP <= n,
	 IP++,
	 JP=NROW[IP-1];
	 If[Abs[A[JP-1,i-1]] > AMAX,
	    AMAX=Abs[A[JP-1,i-1]];
	    IMAX=IP;
	 ];
      ];
      (* Step 13 *)
      If[AMAX <= 10^-20,
	 OK = 0,
         (* Step 14 - Simulate row interchange *)
	 If[NROW[i-1] != NROW[IMAX-1],
	    NCOPY=NROW[i-1];
            NROW[i-1]=NROW[IMAX-1];
	    NROW[IMAX-1]=NCOPY;
	 ];
	 I1=NROW[i-1];
         (* Step 15 - Perform elimination *)
	 For[J=JJ,
	    J <= M,
	    J++,
	    J1=NROW[J-1];
            (* Step 16 *)
	    XM=A[J1-1,i-1]/A[I1-1,i-1];
            (* Step 17 *)
	    For[K = JJ,
	       K <= M,
	       K++,
	       A[J1-1,K-1]=A[J1-1,K-1]-XM*A[I1-1,K-1];
	    ];
            (* Step 18 *)
	    A[J1-1,i-1]=0;
	 ];
      ];
      i=i+1;
   ];
   If[OK == 1,
      (* Step 19 *)
      N1=NROW[n-1];
      If[Abs[A[N1-1,n-1]] <= 10^-20,
	 OK = 0,
         (* System has no unique solution *)
         (* Step 20 - Start backward substitution *)
  	 If[LM > 0,
	    q[LM]=A[N1-1,M-1]/A[N1-1,n-1];
	    A[N1-1,M-1]=q[LM];
	 ];
	 PP=1;
         (* Step 21 *)
	 For[K=LN+1,
	    K <= NN,
	    K++,
	    i=NN-K+LN+1;
	    JJ=i+1;
	    N2=NROW[i-1];
	    SUM=A[N2-1,n];
	    For[KK=JJ,
	       KK <= n,
	       KK++,
	       LL=NROW[KK-1];
	       SUM=SUM-A[N2-1,KK-1]*A[LL-1,M-1];
	    ];
	    A[N2-1,M-1]=SUM/A[N2-1,i-1];
	    q[LM-PP]=A[N2-1,M-1];
	    PP=PP+1;
 	 ];
         (* Step 22 *)
	 For[K=1,
	    K <= LN,
	    K++,
	    i=LN-K+1;
	    N2=NROW[i-1];
	    SUM=A[N2-1,n];
 	    For[KK=LN+1,
	       KK <= n,
	       KK++,
	       LL=NROW[KK-1];
	       SUM=SUM-A[N2-1,KK-1]*A[LL-1,M-1];
	    ];
	    A[N2-1,M-1]=SUM;
	    P[LN-K+1]=A[N2-1,M-1];   
	 ];
         (* Step 23 - Procedure completed successfully *)
	 FLAG = Input["Select output destination\n
                 1. Screen\n
                 2. Text file\n
                 Enter 1 or 2\n"];
         If[FLAG == 2,
            NAME = InputString["Input the file name\n
                    For example:   output.dta\n"];
            OUP = OpenWrite[NAME,FormatType->OutputForm],
            OUP = "stdout"
         ];
	 Write[OUP,"PADE' RATIONAL APPROXIMATION\n"];
	 Write[OUP,"\n"];
	 Write[OUP,"Denominator coefficients Q[0], ..., Q[M]\n"];
	 For[i =0,
	    i <= LM,
	    i++,
	    Write[OUP,SetPrecision[q[i],9]];
	 ]; 
	 Write[OUP,"\n"];
	 Write[OUP,"Numerator coefficients P[0], ..., P[n]\n"];
	 For[i =0,
	    i <= LN,
	    i++,
	    Write[OUP,SetPrecision[P[i],9]];
	 ];
         Write[OUP,"\n"];
         If[OUP == "OutputStream[",NAME," 3]",
            Print["Output file: ",NAME," created successfully\n"];
            Close[OUP];
         ];
      ];
   ];
   If[OK == 0,
      Print["System has no unique solution\n"];
   ];
];
Quit[];