#Homework by Ben Miles for Experimental math assigned 2/9 OK to post
#Gsomosk(Ini,k,n) outputs the n-th term of the Somos-k recurrence, where Somos-k is the natural k-th order analog of the Somos-4 recurrence

Gsomosk:=proc(Ini,k,n) local j,m: option remember:
m:=floor(k/2):
if k<>nops(Ini) then
 RETURN("Wrong number of Initial Conditions"):
elif n<(nops(Ini)+1) then
 Ini[n]:
else:
 add(Gsomosk(Ini,k,n-j)*Gsomosk(Ini,k,n-(k-j)),j=1..m)/Gsomosk(Ini,k,n-k);
fi:
end:

#It seemed to me that only the numbers up to 7 produced just integers. I recorded the first few that broke down (and when) and got the following sequence:
#17, 19, 20, 22, 24, 27 Which gave three results in OEIS including A030127 which confired my feelings about where the break down occurs.
#Obviously k>7 produces non monomial denominators for otherwise there wouldn't be fractions.
#Testing k=7 seemed to imply that k=7 gives Laurent Polynomials. Further research suggested that it is still an open question for k>5

GenPol:=proc(X,d,a,co) 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:
P,var,co1:
end:

GenHomPol:=proc(X,d,a,co) local P,i,Q,b,co1, vars: option remember:
P:=GenPol(X,d,b,co)[1];
Q:=0:
co1:=co:
vars:={}:
if d=0 then
 return a[co],{a[co]},co+1:
else 
 for i from 1 to d+1 do
  Q:=Q+op(i*d+1,P)*a[co1]/b[co+i*d]:
  vars:=vars union {a[co1]}:
  co1:=co1+1:
 od:
fi:
Q, vars, co1:
end:


GenNonHomPol:=proc(X,d,a,co) local i, co1, P,G: option remember:
co1:=co:
P:=0:
for i from 0 to d do
 G:=GenHomPol(X,i,a,co1):
 P:=P+G[1]:
 co1:=G[3]:
od:
P:
end: