-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsuperguesss
147 lines (135 loc) · 3.81 KB
/
superguesss
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
# Harm Derksen
# last change: 7/15/1999
# prim-ed Oct 30 2007
kernelopts(printbytes=false):
`guesss/cfraction`:=proc(list,acc,x)
local n,i,j,a,guess,val,lco,newguess,newval,newlco,result,newresult;
n:=nops(list);
a:=linalg[matrix](n,acc+n,0);
for i from 1 to n do
for j from 1 to acc do
a[i,j]:=coeff(list[i],x,j-1);
od;
a[i,acc+i]:=1;
od;
a:=linalg[submatrix](linalg[gausselim](a),1..n,(acc+1)..(acc+n));
for i from 1 to n do guess[i]:=linalg[row](a,i); od;
result:=linalg[multiply](a,list);
for i from 1 to n do
val[i]:=ldegree(result[i],x);
lco[i]:=tcoeff(result[i],x);
od;
while true do
newguess:=map(expand,linalg[scalarmul](guess[1],x^(val[2]-val[1])));
newresult:=expand(x^(val[2]-val[1])*result[1]);newval:=val[2];
newlco:=lco[1];val[n+1]:=val[n]+1;
for i from 2 to n do
while (newval<val[i+1]) do
newresult:=rem(expand(newresult-newlco/lco[i]*result[i]*
x^(newval-val[i])),x^acc,x);
newguess:=linalg[matadd](newguess,guess[i],1,-newlco/lco[i]*x^(newval-val[i]));
if newresult=0 then RETURN(op(newguess)) fi;
newval:=ldegree(newresult);
newlco:=tcoeff(newresult);
od;
od;
for i from 1 to n-1 do
result[i]:=result[i+1];
guess[i]:=eval(guess[i+1]);
val[i]:=val[i+1];
lco[i]:=lco[i+1];
od;
result[n]:=newresult;
guess[n]:=eval(newguess);
val[n]:=newval;
lco[n]:=newlco;
od;
end:
`guesss/monomials`:=proc(list,d,acc,x)
local i,monsmall,mon,jj;
if list=[] then RETURN([1]) fi;
mon:=[];
for i from 0 to d do
monsmall:=`guesss/monomials`(list[2..nops(list)],d-i,acc,x);
mon:=[op(mon),seq(rem(expand(list[1]^i*monsmall[jj]),x^acc,x),jj=1..nops(monsmall))];
od;
mon;
end:
`guesss/searchbinomial`:=proc(n)
local i,j,c;
c:=0; while(n>=binomial(2*c,c)) do c:=c+1 od;c:=c-1;
for j from c by -1 to 1 do
i:=1;while (n>binomial(j+i,j)) do i:=i+1 od;
if n=binomial(j+i,j) then RETURN([i,j]) fi;
od;
end:
`guesss/makelist`:=proc(f,n,d,acc,x)
local g,list,i;
g:=f;list:=[f];
for i from 1 to n-1 do
g:=diff(g,x);list:=[op(list),expand(x^i*g)]
od;
`guesss/monomials`(list,d,acc,x);
end:
`guesss/formalmonomials`:=proc(list,d,x)
local i,monsmall,mon,jj;
if list=[] then RETURN([1]) fi;
mon:=[];
for i from 0 to d do
monsmall:=`guesss/formalmonomials`(list[2..nops(list)],d-i,x);
mon:=[op(mon),seq(expand(list[1]^i*monsmall[jj]),jj=1..nops(monsmall))];
od;
mon;
end:
`guesss/makeformallist`:=proc(f,n,d,x)
local g,list,i;
g:=f;list:=[f];
for i from 1 to n-1 do
g:=diff(g,x);list:=[op(list),expand(x^i*g)]
od;
`guesss/formalmonomials`(list,d,acc,x);
end:
# GUESS Sequence
guesss:=proc()
local sequ,i,F,j,jj,n,f,g,polyrel,guess,level,lisst,result,eqns,solutions,var,pq,
numextension,x;
sequ:=args[1];
printf("%s\n",`CUT`);
if nargs>1 then
numextension:=args[2]
else
numextension:=7;
fi;
n:=nops(sequ);
for level from 2 to trunc(n/2) do
pq:=`guesss/searchbinomial`(level);
f:=sum(sequ[j]*x^(j-1),j=1..(n-1));
lisst:=`guesss/makelist`(f,pq[1],pq[2],n-1,x);
polyrel:=`guesss/cfraction`(lisst,n-1,x);
g:=f+sum(guess[j]*x^(n-2+j),j=1..numextension);
lisst:=`guesss/makelist`(g,pq[1],pq[2],n+numextension-1,x);
result:=expand(sum(lisst[j]*polyrel[j],j=1..level));
eqns:={seq(coeff(result,x,jj),jj=(n-1)..(n+numextension-2))};
var:={seq(guess[jj],jj=1..numextension)};
solutions:=solve(eqns,var);
solutions:=subs(solutions,[seq(guess[jj],jj=1..numextension)]);
if solutions[1]=sequ[n] then
lisst:=`guesss/makeformallist`(F(x),pq[1],pq[2],x);
printf("%s\n",`SUCCESS`);
printf("%s\n",`Guesss suggests that the generating function F(x) `);
printf("%s\n",`may satisfy the following algebraic or differential equation:`);
printf("%s\n",` `);
lprint(sum(polyrel[i]*lisst[i],i=1..nops(lisst))=0);
printf("%s\n",` `);
printf("%s\n",`If this is correct the next 6 numbers in the sequence are:`);
printf("%s\n",` `);
lprint([op(2..numextension,solutions)]);
printf("%s\n",` `);
RETURN(); fi;
od;
printf("%s\n",`FAILED`);
RETURN();
end:
guesss(lis);
quit;