#C10.txt Feb. 25, 2019 Math640(RU) Spring 2019
Help:=proc(): print(`SeqToC(L) , CtoSeq(C,N), CtoRat(C,x), RatToSeq(R,x,N),  `):
 print(` GuessAlg1(L,F,x,d), GuessAlg(L,F,x) `): end:


##from C8.txt
#CtoSeq(C,N): inputs a C-finite sequence C=[INI,Rec] and a pos. integer N and
  #outputs the first N terms starting at n=0.
  CtoSeq:=proc(C,N) local L,INI,R,nt,i:
   INI:=C[1]: R:=C[2]:
    if not (type(C,list) and nops(C)=2 and type(C[1],list) and type(C[2],list) and
       nops(C[1])=nops(C[2])) then
     print(`bad input`):
     RETURN(FAIL):
    fi:
   L:=INI:
   while nops(L)<=N do
    nt:=add(R[i]*L[-i],i=1..nops(R)):
    L:=[op(L),nt]:
   od:

   L:
  end:

#end from C8.txt
#From hw8.txt (thanks to Victoria Chayes)
#SeqToC1(L,d) inputs a list of integers of at least length 2d+4 and and looks for a linear recurrence with constant coefficients of order d, outputting it in C-finite description C=[INI,Rec]. Initial conditions are assumed to be the first d terms of the list.
SeqToC1:=proc(L,d) local a,i,n,eq, Rec,j:
if nops(L)<2*d+4 then
  print(`not enough data`):
  return(FAIL):
fi:

Rec:={seq(a[i],i=1..d)}:
eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d)=0,n=(d+1)..nops(L))}:
Rec:=solve(eq, Rec):

if Rec=NULL then
  return(FAIL):
fi:

[[op(1..d,L)],[seq(subs(Rec,a[j]),j=1..d)]]:

end:

SeqToC:=proc(L) local i,d,C:
d:=floor((nops(L)-4)/2):
for i from 1 to d do
  C:=SeqToC1(L,i):
  if C<>FAIL then
    return(C):
  fi:
od:
end:

#end from hw8.txt

 #CtoRat(C,x): inputs a C-finite sequence [INI,Rec] and a symbol x
 #outputs the rational generating function of the sequence
 CtoRat:=proc(C,x) local INI,Rec,P,Q,d,i,f,L:
  if not type(x,symbol) then
    print(x, `must be a symbol . `):
    RETURN(FAIL):
  fi: 

  if not (type(C,list) and nops(C)=2 and type(C[1],list) 
     and type(C[2],list) and nops(C[1])=nops(C[2]) ) then
     print(`Bad input`):
      RETURN(FAIL):
  fi:
   INI:=C[1]: Rec:=C[2]:
    d:=nops(INI):
  #f(x)=Sum(a[n]*x^n,n=0..infinity)
  #Sum((a[n]-Rec[1]*a[n-1]- Rec[2]*a[n-2]- ... - Rec[d]*a[n-d])*x^n,n=0..infinity)
  #f(x)(1-Rec[1]*x -Rec[2]*x^2-...-Rec[d]*x^d)=Poly of degree -1
   Q:=1-add(Rec[i]*x^i,i=1..d):
   L:=CtoSeq(C,d):
   f:=add(L[i]*x^(i-1),i=1..nops(L)):
   f:=expand(Q*f):
    P:=add(coeff(f,x,i)*x^i,i=0..d-1):
    factor(P/Q):

   end:

  #RatToSeq(R,x,N): the first N+1 terms of the sequence of coefficients of the rational
   #function R of x starting at n=0.
  RatToSeq:=proc(R,x,N) local R1,f:
   R1:=normal(R):
    if subs(x=0,denom(R1))=0 then
       RETURN(FAIL):
    fi:

    f:=taylor(R1,x=0,N+1):
    [seq(coeff(f,x,i),i=0..N)]:

  end:

  #Def: a sequence a(n) (n=0..infinity) is an ALGEBRAIC sequence if the
  #generating function f(x)=Sum(a(i)*x^i,i=0..infinity) satifies an algebraic equation 
  #P(x,F)==0  (F=f(x))
    
   #GuessAlg1(L,F,x,d): inputs a list of numbers L (starting at n=0)
   #symbols F and x and a pos. number d, outputs a polynomial P, in F and x
   #of degree d such that  P(x,Sum(L[i+1]*x^i,i=0..nops(L)-1) had degree nops(L)
   #improved after class to get rid of annoying factors
   GuessAlg1:=proc(L,F,x,d) local P,i,j,a,eq,var,f,tem,c:
   P:=add(add(a[i,j]*x^i*F^j,    j=0..d-i), i=0..d):
   f:=add(L[i+1]*x^i,i=0..nops(L)-1):
   var:={seq(seq(a[i,j],j=0..d-i),i=0..d)}:
   #nops(var)=(d+1)*(d+2)/2
    if nops(L)<=(d+1)*(d+2)/2+3 then
     print(`Make the list bigger`):
      RETURN(FAIL):
     fi:
   tem:=subs(F=f,P):
    eq:={seq(coeff(tem,x,i)=0,i=0..nops(L)-1)}:
    var:=solve(eq,var):

   

    P:=subs(var,P):
    if P=0 then
      RETURN(FAIL):
     else
     
      P:=normal(P):
      c:=lcoeff(lcoeff(P,F),x):
     P:=normal(P/c):
    fi:

    end:

 
   #GuessAlg(L,F,x): inputs a list of numbers L (starting at n=0)
   #symbols F and x and a pos. number d, outputs a polynomial P, in F and x
   #of degree d such that  P(x,Sum(L[i+1]*x^i,i=0..nops(L)-1) had degree nops(L)
     GuessAlg:=proc(L,F,x) local P,d:

       for d from 1 while nops(L)>=(d+1)*(d+2)/2+3 do
         P:=GuessAlg1(L,F,x,d):
         if P<>FAIL then
                RETURN(P):
         fi:
       od:
 FAIL:
end: