### Homework 7.  Patrick Devlin.  February 11, 2012  ###
### I have no preferrences about keeping my work private ###


### Programs from classes 7 and 6 at the bottom ###

Help:=proc():
  print(`From homework 7:  Gsomosk(Ini, k, n),  GenHomPol(X,d,a,co:=1), GenNonHomPol(X,d,a,co:=1)`):
  print(`From class 7:  Gsomos4(Ini,n),  mysteryFib(n),  Guess2NLR(L,f,n,maxD),  GenPol(X,d,a,co:=1)`):
  print(`From class 6:  GuessRec1(L,r),  GuessRec(L) `):
end proc:




#############
# Problem 1 #
#############


#Gsomosk(Ini,k,n): outputs the n-th term of the Somos-k recurrence with initial conditions Ini, i.e. sk[1]=Ini[1], etc.

Gsomosk:=proc(Ini,k,n,symb:=false) option remember: local i,x:
 if (n<=k) then if((n<=nops(Ini)) and (n>0)) then return Ini[n]: else if(symb) then return x[n]: else return 1: fi: fi: fi:
 if(symb) then
   return normal((add(Gsomosk(Ini,k,n-i, symb)*Gsomosk(Ini,k,n-k+i, symb), i=1..(k/2)))/Gsomosk(Ini,k,n-k, symb)):
 else
   return (add(Gsomosk(Ini,k,n-i, symb)*Gsomosk(Ini,k,n-k+i, symb), i=1..(k/2)))/Gsomosk(Ini,k,n-k, symb):
 fi:
end proc:


#############
# Problem 2 #
#############

# Running a function such as
# 
# test3:=(k,stop)-> {seq(denom(Gsomosk([1$k], k, i)), i=1..stop)}
# 
# Seems to reveal that Gsomosk([1$k], k, i) is always an integer
# if and only if 2 <= k <= 7.
# 
# If k is one of those values, the experimental data obviously cannot
# rigorously prove that Gsomosk(...)  is always an integer.
# However, we have found non-integer values for   8 <= k <= 35,
# which suggests the conjecture.



#############
# Problem 3 #
#############

# As before, by running a function such as
#  
# test4:=(k,stop)-> {seq(evalb( op(0,expand(denom(Gsomosk([], k, i, true)))) <> `+`), i=1..stop)}:
# 
# seems to reveal that Gsomosk([], k, i, true) always produces Laurent
# polynomials if and only if  2 <= k <= 7.
# 
# However, this is harder to test since maple is slower with the purely symbolic manipulations



#############
# Problem 4 #
#############

# Takes a list of variables X, the total degree d, a list of coefficients a,
# and an optional parameter co, which keeps track of depth in recursive calls.
# 
# A modified version of GenPol.  It returns a polynomial in the variables
# X[1], ... , X[nops(X)],  where each term is a monomial of total degree d.

GenHomPol:=proc(X,d,a,co:=1) local P, var, co1,m,i,x,X1,P1: option remember:
 m:=nops(X):
 x:=X[m]:

 if m=1 then
  RETURN(x^d * a[co], {a[co]},co+1):
 fi:

 co1:=co:
 X1:=[op(1..m-1,X)]:
 P:=0:
 var:={}:

 for i from 0 to d do
   P1:=GenHomPol(X1,d-i,a,co1):
   P:=expand(P+x^i * P1[1]):
   var:=var union P1[2]:
   co1:=P1[3]:
 od:

 return (P,var,co1):
end:




#############
# Problem 5 #
#############

# Takes a list of variables X, a maximum total degree d, a variable name for
# the list of coefficients, a,  and an option argument co, indicating which
# index of the list of variables should be the first used for the coefficients.
# 
# Outputs a general polynomial in X, whose maximum total degree is exactly d

GenNonHomPol:= proc(X,d,a, co:=1) local P, var, i, tempAns, co1:
  P:=0:  var:={}:  co1:=co:
  for i from 0 to d do
    tempAns:=GenHomPol(X,i,a,co1):
    P:=expand(P+tempAns[1]):
    var:=var union tempAns[2]:
    co1:=tempAns[3]:
  od:
  return (P, var, co):
end proc:



######################
# Stuff from class 7 #
######################


#Gsomos4(Ini,n): outputs the n-th term of the Somos4 recurrence with initial conditions Ini, i.e. s[1]=Ini[1], etc.
# For example, Gsomos4([1,1,1,1,20);
Gsomos4:=proc(Ini,n) option remember:
 if n<=4 then return Ini[n]: fi:
 return normal((Gsomos4(Ini,n-1)*Gsomos4(Ini,n-3)+Gsomos4(Ini,n-2)^2)/Gsomos4(Ini,n-4)):
end proc:


#Guess2NLR(L,f,n, MaxD) inputs a sequene L and symbols f and n and tries to guess
# a NON-LINEAR recurrene of 2nd order of the form P( f(n-2),f(n-1),f(n))=0 of degree MaxD

Guess2NLR:=proc(L,f,n, MaxD) local d,ans:
 for d from 1 to MaxD do
   ans:=Guess2NLRd(L,f,n,d):
   if ans<>FAIL then return ans: fi:
 od:
 return FAIL:
end:

## GenPol(X,d,a): 
#
## Inputs 
#   (i)    a list of variables X e.g. [x,y,z] (let m=nops(X), the number of variables)
#   (ii)   a pos. integer d, for the degree
#   (iii)  a symbol a by which the coeff. (undetermined) of the pol. are expressed
#   (iv)   co, starting counter
#
## Outputs: 
#   (i)    a generic polynomial in the variables X of degree d with coeff. expressed as a[co],ac[co+1], ...,  a[co+AsNeeded]
#   (ii)   the set of coefficients
#   (iii)  the new value of the counter
#
## For example, GenPol([x],1,a,0) --outputs-->   a[0]+a[1]*x, {a[0],a[1]}, 2

GenPol:=proc(X,d,a,co:=1) local P, var, co1,m,i,x,X1,P1: option remember:
 m:=nops(X):
 x:=X[1]:

 if m=1 then
  RETURN(add( a[co+i]*x^i,i=0..d), {seq(a[co+i],i=0..d)},co+d+1):
 fi:

 co1:=co:
 X1:=[op(2..m,X)]:
 P:=0:
 var:={}:

 for i from 0 to d do
   P1:=GenPol(X1,d,a,co1):
   P:=expand(P+P1[1]*x^i):
   var:=var union P1[2]:
   co1:=P1[3]:
 od:

 return (P,var,co1):
end:



######################
# Stuff from Class 6 #
######################

# GuessRec1(L,r):  given a list of numbers and a positive integer, r,
#  it tries to guess a linear recurrence relation of order r
#  that fits the list

# Looking for a recurrence of the form   L[n] + c_1 L[n-1] + ... + c_r L[n-r] = 0   for all relevant n
# We need r data points to determine this system of (c_1, c_2, ... , c_{r})  unknowns

# So we'll output [ [c_1, c_2, c_3, ... , c_r], [initial_conditions]].
# (i.e., Fibonacci would be  F[n] - F[n-1] - F[n-2] = 0   -> [[-1, -1], [1,1]]

GuessRec1:=proc(L,r) local var, c, i, n, eqns:
 if (nops(L) <= 2*r + 6) then return FAIL; fi: # The list isn't long enough, so return fail
 var:={seq(c[i], i=1..r)}:
 eqns:=[ seq( (L[n] + sum(c[i] * L[n-i],i=1..r))=0, n=(r+1)..nops(L))]:
 var:=solve(eqns,var);
 if(var = NULL) then return FAIL: fi:
 return [[seq(eval(c[i],var),i=1..r)],L[1..r]];
end proc:

# GuessRec(L) returns a linear recurrence relation of minimal order that fits L

GuessRec:=proc(L) local r, ans:
 for r from 1 to (nops(L)-6)/2 do
  ans:=GuessRec1(L,r):
  if (ans <> FAIL) then return ans: fi:
 od:
 return FAIL:
end proc: