#John Kim
#Use whatever you like



Help:=proc(): print(` AlgToSeq(F,P,x,K,f0)`): end:

Help10:=proc(): print(` GuessAlg(L,P,x,MaxDegP) `): 
print( ` GuessAlg1(L,P,x,DegP) , GuessAlg11(L,P,x,DegP,Degx) `):

end:

###stuff from C10.txt###########
#C10.txt: starting the Alg. Formal Power Series
#1+1!*x+2!*x^2+3!*x^3+ ......

#GuessAlg(L,P,x,MaxDegP): Inputs a sequence of numbers L
#and tries to guess  an alg. eq. of the from
#F(P,x)=0. where P is the symbol for the generating
#function of L, P(x)=add(L[i+1]*x^i,i=0..infinity)
#of max. degree MaxDegP
GuessAlg:=proc(L,P,x, MaxDegP) local g,DegP:

for DegP from 1 to MaxDegP do
 g:=GuessAlg1(L,P,x,DegP):
 if g<>FAIL then
     RETURN(g):
 fi:
od:

FAIL:

end:


#GuessAlg1(L,P,x,DegP): Inputs a sequence of numbers L
#and tries to guess  an alg. eq. of the from
#F(P,x)=0. where P is the symbol for the generating
#function of L, P(x)=add(L[i+1]*x^i,i=0..infinity)
#of degree DegP in P
GuessAlg1:=proc(L,P,x, DegP) local g,Degx:



for Degx from 0 while (1+Degx)*(1+DegP)<nops(L)-6 do

g:=GuessAlg11(L,P,x, DegP,Degx):

if g<>FAIL then
  RETURN(g):
fi:

od:

FAIL:

end:



#GuessAlg11(L,P,x,DegP, Degx): Inputs a sequence of numbers L
#and tries to guess  an alg. eq. of the from
#F(P,x)=0. where P is the symbol for the generating
#function of L, P(x)=add(L[i+1]*x^i,i=0..infinity)
#of degree DegP in P and Degx in x
GuessAlg11:=proc(L,P,x, DegP, Degx) local g, F,a, i,j,
var,P1,F1,eq, var1,v:

F:=add(add(a[i,j]*x^i*P^j,i=0..Degx),j=0..DegP):

var:={seq(seq(a[i,j],i=0..Degx),j=0..DegP)}:

P1:=add(L[i]*x^(i-1),i=1..nops(L)):

F1:=expand(subs(P=P1,F)):

eq:={seq( coeff(F1,x,i)=0, i=0..nops(L)-1)}:

var1:=solve(eq,var):

F:=subs(var1,F):

if F=0 then
 RETURN(FAIL):
else
F:=subs({seq(v=1, v  in var)},F):
 
F:=add(factor(coeff(F,P,i))*P^i,i=0..degree(F,P)): 
RETURN(F):
fi:



end:

##end stuff from C10.txt



#new stuff starts
#AlgToSeq(F,P,x,K,f0): INPUTS 
#(i) a polynomial F in the variables P and x
#denoting an algebraic relation satisfied by a formal
#power series P(x), such that
#F(P(x),x)=0 
#(ii) and (iii) variables P and as above
#(iv): a pos. integter K
#(v): f0 the constant term
#OUTPUTS:
#a sequence of length K+1 denoting the
#coeff. of x^0 through x^K of the  Maclaurin expansion
#of the power series of P(x).
#For example,
#AlgToSeq(P*(1-x)-1,P,x,4); should output
#[1,1,1,1,1];
#AlgToSeq(P-1-x*P^2,P,x,4); should yield:
#[1,1,2,5,14]
AlgToSeq:=proc(F,P,x,K,f0) local H,F1,i,f,L,m:

H:=coeff(F,P,1):

if subs(x=0,H)=0 then
 RETURN(FAIL):
fi:

#a0(x)+H(x)*P(x)+a2(x)*P(x)^2+...=0
#P(x)=-(a0(x)+a2(x)*P(x)^2+....)/H(x)

F1:=(P*H-F)/H:

#start at f(x)=1 and iterate f(x)->F1(f(x))

f:=f0:
L:=[f0]:

for i from 1 to K do
f:=subs(P=f,F1):

f:=taylor(f,x=0,i+1):



L:=[op(L),coeff(f,x,i)]:

f:=add(coeff(f,x,m)*x^m,m=0..i):

od:
L:
end:




#GuessDE(L,x,Di,Max): inputs a list of L, representing the first nops(L) terms of a formal power series
#P(x)=L[1]+L[2]*x+ ...+L[nops(L)]*x^(L-1)+ ...
#and outputs a differential operator with polynomial coefficients
#OPER(x,Di)
#(where Di is differentiation with respect to x) such that the ORDER ("degree" in Di) and DEGREE (degree in x) satisfy
#DEGREE+ORDER<=Max, but with DEGREE+ORDER as small as possible, such that
#when you apply OPER(x,Di) to P(x) it is (conjenturally) 0 
#(and surely the first L coefficients of the Maclaurin expansion of OPER(x,Di) should all vanish. 
#Otherwise it is to return FAIL.
#For example,
#GuessDE([seq(1/i!,i=0..20)],x,Di,1);
#should return Di-1 

GuessDE:=proc(L,x,Di,Max) local M,oper:

for M from 1 to Max do:
 oper:=GuessDE1(L,x,Di,M):
 if oper<>FAIL then
  RETURN(oper):
 fi:
od:

FAIL:

end:



#GuessDE1(L,x,Di,M): inputs a list of L, representing the first nops(L) terms of a formal power series
#P(x)=L[1]+L[2]*x+ ...+L[nops(L)]*x^(L-1)+ ...
#and outputs a differential operator with polynomial coefficients
#OPER(x,Di)
#(where Di is differentiation with respect to x) such that the ORDER ("degree" in Di) and DEGREE (degree in x) satisfy
#DEGREE+ORDER=Max, such that
#when you apply OPER(x,Di) to P(x) it is (conjenturally) 0 
#(and surely the first L coefficients of the Maclaurin expansion of OPER(x,Di) should all vanish. 
#Otherwise it is to return FAIL.
#For example,
#GuessDE1([seq(1/i!,i=0..20)],x,Di,1);
#should return Di-1  

GuessDE1:=proc(L,x,Di,M) local f,oper,dif,i,j,a,var,var1,eq:

f:=add(L[i]*x^(i-1),i=1..nops(L)):
oper:=add(add(a[i,j]*x^i*Di^j,j=0..M-i),i=0..M):
dif:=add(a[i,0]*x^i*f,i=0..M)+add(add(a[i,j]*x^i*diff(f,x$j),j=1..M-i),i=0..M):

var:={seq(seq(a[i,j],j=0..M-i),i=0..M)}:
eq:={seq( coeff(dif,x,i)=0, i=0..nops(L)-M-1)}:

var1:=solve(eq,var):
oper:=subs(var1,oper):

if oper=0 then
 RETURN(FAIL):
else
 oper:=subs({seq(v=1, v  in var)},oper):
 RETURN(oper):
fi:

end:



#AlgToDE(F,P,x,Di,K) that inputs a polynomial F(P,x) in the variables P and x, representing an algebraic equation
#F(P(x),x)=0 satisfied by a formal power series P(x), and a positive integer K, and conjenctures a differential 
#operator "annihilating" P(x), by using the first K+1 coefficients of the Maclaurin expansion of P(x).
#For example, AlgToDE(P*(1-x)-1,P,x,Di,30);
#should return (1-x)*Di-1. 

AlgToDE:=proc(F,P,x,Di,K) local L:

L:=AlgToSeq(F,P,x,K,1):
GuessDE(L,x,Di,nops(L)):

end: