######################################################################
##PADE: Save this file as PADE. To use it, stay in the               #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read PADE : <Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg@math.Rutgers.edu.                                          # 
#######################################################################

#Created: Feb. 16, 2002
#This version: Feb. 16, 2002
#PADE: A Maple package to study the PADE Appx.
#Please report bugs to zeilberg@math.Rutgers.edu



print(`Created: Feb. 16, 2002.`):
print(`This version: Feb. 16, 2002`):
lprint(``):
print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`):
lprint(``):
print(`Please report bugs to zeilberg@math.rutgers.edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/`):
 print(`For a list of the procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(``):
ezra:=proc()
if args=NULL then
 print(`Contains the following procedures: delt, DeltaSeq,`):
 print(` DiagSeq, Pade1, Pade2 , PolTab `):
fi:

if nops([args])=1 and op(1,[args])=delt then
print(`delt(a,c): Given a rational number a and a constant`):
print(` c, finds the delta such that |a-c|=1/denom(a)^(1+delta) .`):
print(` For example, try delta(22/7,evalf(Pi)): `):
fi:

if nops([args])=1 and op(1,[args])=DeltaSeq then
print(`DeltaSeq(f,x,M,x0): the first M terms of the`):
print(`sequence of deltas for the rational numbers obtained`):
print(`by plugging in x=x0 in the diagonal of the Pade table for`):
print(`f(x) : For example, try:DeltaSeq(exp(x),x,10,1); .`):
fi:

if nops([args])=1 and op(1,[args])=DiagSeq then
print(`DiagSeq(f,x,M,x0): the first M terms of the`):
print(`sequence of rational numbers obtained`):
print(`by plugging in x=x0 in the diagonal of the Pade table for`):
print(`f(x) : For example, try:DiagSeq(exp(x),x,10,1); .`):
fi:

if nops([args])=1 and op(1,[args])=Pade1 then
print(`Pade1(f,x,M,N): Given an expression f(x) and two integers, M, and N`):
print(`finds the best rational appx. with numerator of degree M and `):
print(`denominator of degree N for example, try: Pade1(exp(x),x,3,5);`):
fi:

if nops([args])=1 and op(1,[args])=Pade2 then
print(`Pade2(f,x,y,lam,mu,nu): Given an expression f(x,y) `):
print(` and two partitions  `):
print(`finds the best rational appx. with numerator `):
print(`of degree M and denominator of degree N for example, try: `):
print(` Pade2(log(1+x+y),x,y,[1,1],[1,1],[3,3]);`):
fi:

if nops([args])=1 and op(1,[args])=PolTab then
print(`PolTab(lam,x,y,a): Given a  sequence of positive integers`):
print(`(usually a tableau), and variables, x,y and a letter a`):
print(`returns a generic polynomial p_0(x)+p_1(x)*y+...+p_r(x)*y^r`):
print(`such that the degree of p_i(x) in x  is lam[i+1],`):
print(`followed by the set of coefficients, for example,`):
print(`try PolTab([3,2],x,y,a);`):
fi:

end:


#Pade1(f,x,M,N): Given an expression f(x) and two integers, M, and N
#finds the best rational appx. with numerator of degree M and denominator
#of degree N for example, try: Pade1(exp(x),x,3,5);
Pade1:=proc(f,x,M,N)
local P,Q,a,b,eq,var,i,f1,f2:
f1:=taylor(f,x=0,M+N+3):
f2:=0:
for i from 0 to M+N+2 do
f2:=f2+coeff(f1,x,i)*x^i:
od:
eq:={}:
P:=0:
Q:=1:
var:={}:
for i from 0 to M do
var:=var union {a[i]}:
P:=P+a[i]*x^i:
od:

for i from 1 to N do
var:=var union {b[i]}:
Q:=Q+b[i]*x^i:
od:

f2:=expand(f2*Q-P):

eq:={}:

for i from 0 to M+N do
eq:=eq union {coeff(f2,x,i)}:
od:

var:=solve(eq,var):

P:=subs(var,P):Q:=subs(var,Q):
normal(P/Q):

end:



#delt(a,c): Given a rational number a and a constant
#c, finds the delta such that |a-c|=1/denom(a)^(1+delta)
#For example, try delta(22/7,evalf(Pi)):
delt:=proc(a,c): evalf(-log(abs(c-a))/log(denom(a)))-1: end:

#DiagSeq(f,x,M,x0): the first M terms of the
#sequence of rational numbers obtained
#by plugging in x=x0 in the diagonal of the Pade table for
#f(x) : For example, try:DiagSeq(exp(x),x,10,1); .
DiagSeq:=proc(f,x,M,x0)
local gu,i,mu:
gu:=[]:
for i from 1 to M do
mu:=Pade1(f,x,i,i):
gu:=[op(gu),subs(x=x0,mu)]:
od:
gu:
end:
 
#PolTab(lam,x,y,a): Given a  sequence of positive integers
#(usually a tableau), and variables, x,y and a letter a
#returns a generic polynomial p_0(x)+p_1(x)*y+...+p_r(x)*y^r
#such that the degree of p_i(x) in x  is lam[i+1],
#followed by the set of coefficients, for example,
#try PolTab([3,2],x,y,a);
PolTab:=proc(lam,x,y,a)
local pol,gu,i,j,co:
pol:=0:
gu:={}:
co:=0:
for i from 1 to nops(lam) do
for j from 0 to lam[i] do
co:=co+1:
pol:=pol+a[co]*y^(i-1)*x^j:
gu:=gu union {a[co]}:
od:
od:
pol,gu:
end:

#godel(lam): the "size" of the array lam, i.e. the sums of
#its entries plus nops(lam)
godel:=proc(lam) convert(lam,`+`)+nops(lam):end:

#Pade2(f,x,y,lam,mu,nu): Given an expression f(x,y) 
#and two partitions 
#finds the best rational appx. with numerator of degree M and denominator
#of degree N for example, try: 
#Pade2(log(1+x+y),x,y,[1,1],[1,1],[3,3]);
Pade2:=proc(f,x,y,lam,mu,nu)
local M,N,gu,F,F1,i,j,TOP1,BOT1,var,eq,a,b,lu,var1:
if godel(lam)+godel(mu)<>godel(nu)+1 then
ERROR(`The sizes of lam and mu should add up to nu's+1`):
fi:
N:=max(nops(lam),nops(mu),nops(nu))+2:
M:=max(op(lam),op(mu),op(nu))+2:

gu:=0:
F:=taylor(f,x=0,M+1):
for i from 0 to M do
F1:=coeff(F,x,i):
F1:=taylor(F1,y=0,N+1):
for j from 0 to N do
gu:=gu+coeff(F1,y,j)*x^i*y^j:
od:
od:

TOP1:=PolTab(lam,x,y,a):
var:=TOP1[2]:
TOP1:=TOP1[1]:
BOT1:=PolTab(mu,x,y,b):
var:=var union BOT1[2]:
BOT1:=BOT1[1]:

gu:=expand(BOT1*gu-TOP1):
eq:={}:

for i from 1 to nops(nu) do
 lu:=coeff(gu,y,i-1):

 for j from 0 to nu[i] do
  eq:=eq union {coeff(lu,x,j)}:
 od:
od:


var1:=solve(eq,var):

#normal(subs(var1,TOP1/BOT1)):
subs(var1,TOP1),subs(var1,BOT1):



end:




#DeltaSeq(f,x,M,x0): the first M terms of the
#delta's for the sequence of rational numbers obtained
#by plugging in x=x0 in the diagonal of the Pade table for
#f(x) : For example, try:DeltaSeq(exp(x),x,10,1); .
DeltaSeq:=proc(f,x,M,x0)
local gu,i,mu:

gu:=[]:
for i from 4 to M do
mu:=Pade1(f,x,i,i):
#print(subs(x=x0,mu)):
#print(delt(subs(x=x0,mu),subs(x=x0,f))):
gu:=[op(gu),delt(subs(x=x0,mu),subs(x=x0,f))]:
od:
evalf(gu,6):
end: