with(combinat):
Help:=proc()

if (args=NULL) then
print(`Functions included are:`):
print(`guessRat(fns,L,maxdeg,x)`):
print(`TryTel(f,fns,numTerm,maxdeg,x)`):
print(`approxLim(summand,n)`):
print(`IPowWScore(number,irrats,maxpow,threshold)`):
print(`findSumRat(summand,n,maxdegree)`):
print(`findSumGenRat(summand,atoms,n,maxdegree)`):
print(`findSumHyper(summand,n,maxdegree)`):
print(`findSumMultiBasic(summand,atoms,n,maxdegree)`):

elif(args=guessRat) then
print(`x is left as a varaible occurring in`):
print(`fns which is a set of base expressions`):
print(`then, this procedure will attemt to return`):
print(`an expression that is rational in those`):
print(`and has maximum total degree of maxdeg`):
print(`for both top and bottom, that agrees with L`):

elif(args=TryTel) then
print(`f is the formula for the summand`):
print(`fns are the basic building blocks to look for in the solution`):
print(`numTerm is how long of a partial sum to consider`):
print(`maxdeg is the highest degree of a building block to use in the solution`):
print(`x(left as a symbol) is the variable appearing in the functions`):
print(`This function attempts to find an expression for the partial sums of summation f`):
print(` that starts at x=1 going to infinity`):
print(`for example to compute partial sums of csch(2^n), you would set e:=exp(1) and run`):
#print(`TryTel(2/(e^(2^n) -e^(-2^n)),{n,e^(2^n),e^(-2^n)},10,1,n);`):
print(`TryTel(-(n^2+13*n-3)/((n^2+n+3)*(n^2-n+3)),{n},20,2,n);`);
print(`or:`):
print(`TryTel((4*(2^(n-1)-2^n-1))/((3+n+2^n)*(2+n+2^(n-1))),{n,2^n},10,1,n);`):

elif(args=approxLim) then
print(`summand is an expression in n(left as a symbol)`):
print(`This function tries to estimate a rational that this sum from 1 to infinity evaluates to`):
print(`for example, approxLim(-(n^4-2*n^3-13*n^2+8*n+3)/((n^3-n+3)*(n^3-3*n^2+2*n+3)),n) = 5/3`):
print(`if it runs for a long time, it will return FAIL, this is what will happen if`):
print(`it is a REALLY messy rational, or, more likely, not rational`):
print(`The summand doesn't need to be rational, for example, try running:`):
print(`approxLim((3*(2^n-2^(n-1)+1))/((3+n+2^n)*(2+n+2^(n-1))),n);`):

elif(args=IPowWScore) then
print(`number is the number to be approximated`):
print(`irrats is the set of irrational numbers whoose powers should be considered`):
print(`maxpow is the highest that the sum of the absolute value of the powers can be`):
print(`threshold determines how much of a jump in the P.F. decomp is needed to stop`):
print(`higher values of threshold generally allow more complicated fractions to appear`):
print(`for example:`):
print(`IPowWScore(add(1/n^2,n=1..200),{Pi},3,50);`):
print(`or `):
print(`IPowWScore(add(1/n^4,n=1..400),{Pi},4,100);`):

elif(args= findSumGenRat) then
print(`tries to find a rational expression in n with degree at most maxdegree`):
print(`from 1 to n of summand`):
print(`for example:`):
print(`findSumRat(1/(n^2+n),n,1);`):

elif(args= findSumRat) then
print(`tries to find a rational expression in atoms, each of which are expressions in n`):
print(`with degree at most maxdegreefrom 1 to n of summand`):
print(`for example:`):
print(`expr:= (2^n-2+6^n)/((2^n)^2-3*2^n+6);`):
print(`summand:= normal(expr - subs(n=n-1,expr));`):
print(`findSumGenRat(summand,{2^n,3^n},n,2);`):

elif(args = findSumHyper) then
print(`tries to find a formula for the summation of summand from 1 to n`):
print(`which is the sum of a constant and a hypergeometric which has ratio with`):
print(`degree at most maxdegree`):
print(`for example:`):
print(``):

fi:
end:

TryTel:=proc(f,fns,numTerm,maxdeg,x) local i,maxI,guess:
guess :=normal(guessRat(fns,[seq(sum(subs(x=i,f),i=1..maxI),maxI=1..numTerm)],maxdeg,x)-approxLim(f,x)):
#guessRat(fns,[seq(sum(subs(x=i,f),i=1..maxI),maxI=1..numTerm)],maxdeg,x):
if(0<>normal(guess - subs(x=x-1,guess)-f)) then
print(`failed check`):
fi:
guess;
end:


#tries to find a polynomial of max degree at most deg with variables taken from fns that results in L
#pretty generous in terms of allowing false positives, only requires one more data point than degree of freedom
guessRat:=proc(fns,L,maxdeg,x) local deg:
for deg from 1 to maxdeg do
if(guessRat2(fns,L,deg,x)<>FAIL) then
return(simplify(guessRat2(fns,L,deg,x))):
fi:
od:
FAIL:
end:

guessRat2:=proc(fns,L,deg,x) local i,j,k,num,den,vars,sol,eqs,topCon,botCon,Var: 
if(nops(L)<=2+2*add(nops(composition(nops(fns)+i,nops(fns))),i=1..deg)) then
print(`need more data to guess polynomials degree `,deg):
return(FAIL):
fi:
vars:='vars';
botCon:='botCon';
topCon:='topCon';
num := add(add(vars[i][j]*mul(fns[k]^(composition(nops(fns)+i,nops(fns))[j][k]-1),k=1..nops(fns)),j=1..nops(composition(nops(fns)+i,nops(fns)))),i=1..deg)+topCon:
den := add(add(vars[i][-j]*mul(fns[k]^(composition(nops(fns)+i,nops(fns))[j][k]-1),k=1..nops(fns)),j=1..nops(composition(nops(fns)+i,nops(fns)))),i=1..deg)+botCon:
#firstly, we need that it evaluates to L
#secondly, we need that we are not dividing by zero
eqs:={seq(subs(x=i,1 = botCon),i=1..nops(L)),seq(subs(x=i,num/den) = L[i],i=1..nops(L))}:
Var:= {topCon,botCon,seq(seq(op({vars[i][j],vars[i][-j]}),j=1..nops(composition(nops(fns)+i,nops(fns)))),i=1..deg)}:
sol:=solve(eqs,Var):
#print(sol):
if(sol=NULL) then
return(FAIL):
fi:
simplify(subs(sol,num/den)):
end:

verifyRat:=proc (summand,guess,x)
return(evalb(normal(guess - subs(x=x-1,guess))=normal(summand)))
end:

approxLim:=proc(summand,n) local i,top,guess:
Digits:=10000;
top:=10:
guess:=identifyRat(add(subs(n=i,summand),i=1..top),100);
top:=top*10;
while(guess<>identifyRat(add(subs(n=i,summand),i=1..top),100)) do
guess:=identifyRat(add(subs(n=i,summand),i=1..top),100);
top:=top*10;
if (top>100000) then
print(`this is gonna take a long time, limit is probably not rational number`):
return(FAIL):
fi:
od:
guess;
end:

identifyRat:=proc(number,threshold) local q,fr,i,intpart,mynum:
Digits:=10000;
intpart:=trunc(number):
mynum:=number-intpart: 
fr :=convert(mynum,'confrac','q'):
i:=0:
while((i<nops(q)) and (fr[i+1]<threshold)) do
i:=i+1:
od:
intpart+q[i]:
end:

#idenifies rational with score. Higher score means more likely rational
IRatWScore:=proc(number,threshold) local intpart,q,mynum,i,fr:
Digits:=10000;
intpart:=trunc(number):
mynum:=number-intpart: 
fr :=convert(mynum,'confrac','q'):
i:=0:
while((i<nops(q)) and (fr[i+1]<threshold)) do
i:=i+1:
od:
if i< nops(fr) then
[intpart+q[i],evalf(fr[i+1]/(3+intpart+i^2*Statistics[Mean](fr[1..i])),10),fr]:
else
[intpart+q[i],0,fr]:
fi
end:

#Tries to identify the number as a rational times some powers of the numbers in irrats
#It will output the number and a score that if it's bigger, we are more sure, if smaller, less sure
IPowWScore:=proc(number,irrats,maxpow,threshold) local i, bestscore,num,score,bestnum,foo,bestfra,fra:
#print(number,irrats,maxpow);
if nops(irrats)=0 then
RETURN(IRatWScore(number,threshold)):
fi:
Digits:=10000:
bestscore :=-1:
bestnum:=1:
bestfra:=[]:
for i from -maxpow to maxpow do
#foo:=IRatWScore(number/mul(irrats[k]^(i-comps[j][k]),k=1..nops(irrats)),threshold):
foo :=IPowWScore(number/irrats[1]^i,irrats[2..nops(irrats)],maxpow-abs(i),threshold):
num:=foo[1]:
score:=foo[2]:
fra:=foo[3]:
if score > bestscore then
bestscore := score:
bestnum := num*irrats[1]^i:
bestfra:=fra:
fi:
od:
[bestnum,bestscore,bestfra]:
end:

brookeSum:= proc(summandPoly,summandHyper,n,N,irrats,threshold) local estimate,newSummandPoly,i:
Digits:=10000:
estimate := IPowWScore(evalf(add(1+subs(n=i,summandPoly*summandHyper),i=1..N)),irrats,1,threshold);
print(estimate);
newSummandPoly := summandPoly - estimate[1];
sum(newSummandPoly*summandHyper,n=1..infinity);
end:



polyForm:= proc(degree,n,a) local i:
[add(a[i]*n^i,i=0..degree),{seq(a[i],i=0..degree)}]
end:

rationalForm:=proc(degree,n,a,b) local out1,out2:
out1 := polyForm(degree,n,a):
out2 := polyForm(degree,n,b):
[out1[1]/out2[1],out1[2] union out2[2]]
end:

multiForm:=proc(maxdeg,ns,a) local deg,comps,comp,ret,vars:
ret:=0;
vars:={}:
for deg from 0 to maxdeg do
comps := composition(deg+nops(ns),nops(ns)):
for comp in comps do
ret:=ret+a[comp]*convert(ns^~(comp-~1),`*`):
vars:=vars union{a[comp]}:
od:
od:
[ret,vars]:
end:

multiRatForm:=proc(degree,ns,a,b) local out1,out2:
out1 := multiForm(degree,ns,a):
out2 := multiForm(degree,ns,b):
[out1[1]/out2[1],out1[2] union out2[2]]
end:

findSumRat:= proc(summand,n,maxdegree) local degree,out,rat,a,b:
for degree from 1 to maxdegree do:
a:='a':
b:='b':
out := rationalForm(degree,n,a,b):
rat :=out[1]:
out :=guessSum(summand,n,k->subs(n=k,out[1]),out[2]):
if(out<>FAIL) then
#print(rat):
rat:=simplify(subs(out,rat)):
if(simplify(normal(rat-subs(n=n-1,rat)) / normal(summand))<>1) then
print(`this only works for first many terms, might not be really true`):
fi:
RETURN(rat):
fi:
od:
end:

findSumGenRat:=proc(summand,atoms,n,maxdegree) local rat,degree,out,a,b:
for degree from 1 to maxdegree do:
out := multiRatForm(degree,atoms,a,b):
rat :=out[1]:
out :=guessSum(summand,n,k-> subs(n=k,rat),out[2]):
if(out<>FAIL) then
#print(rat):
rat:=simplify(subs(out,rat)):
if(simplify(normal(rat-subs(n=n-1,rat)) / normal(summand))<>1) then
print(`this only works for first many terms, might not be really true`):
fi:
RETURN(rat):
fi:
od:
end:

#tries to get that the sum is a constant plus a hypergeometric
#This is exactly the scope that the real Gosper's algorithm would handle
findSumHyper:=proc(summand,n,maxdegree)local degree,out,oute,rat,a,b,c,d,start,solution:
for degree from 1 to maxdegree do:
a:='a':
b:='b':
out := rationalForm(degree,n,a,b):
#print(out[1]);
rat:=out[1]:
c:='c':
start:=polyForm(degree,n,c):
d:='d':
oute :=guessSum(summand,n,k->d+hyperGeoForm(rat,start[1],n)(k),{d,op(start[2]),op(out[2])}):
#print(oute):
if(oute<>FAIL) then
print(`it has constant term `):
print(subs(oute,d)):
print(`and hypergeometric ratio `):
print(simplify(subs(oute,rat))):
print(`starting with`):
print(simplify(subs(oute,start[1]))):
solution :=subs(oute,d)+simplify(summand * (subs(oute,rat))/(subs(oute,rat) - 1)):
RETURN(solution):
fi:
od:
print(`failed`):
end:

#tries to get that the sum is a constant plus a hypergeometric
#This is exactly the scope that the real Gosper's algorithm would handle
findMultiBasic:=proc(summand,atoms,n,maxdegree)local degree,out,oute,rat,a,b,c,d,start,solution:
for degree from 1 to maxdegree do:
a:='a':
b:='b':
out := multiRatForm(degree,atoms,a,b):
c:='c':
start:= multiForm(degree,atoms,c):
#print(out[1]);
rat:=out[1]:
d:='d':
oute :=guessSum(summand,n,k->d+multiBasicForm(rat,start[1],n)(k),{d,op(start[2]),op(out[2])}):
#print(oute):
if(oute<>FAIL) then
print(`it has constant term `):
print(subs(oute,d)):
print(`and ratio of consecutive terms `):
print(simplify(subs(oute,rat))):
print(`starting with`):
print(simplify(subs(oute,start[1]))):
solution :=subs(oute,d)+simplify(summand * (1)/(subs(oute,rat) - 1)):
RETURN(solution):
fi:
od:
print(`failed`):
end:

guessSum:=proc(summand,n,formOfSolution,vars) local i,eqs,sols,N,j:
N := nops(vars)+3;
eqs:={}:
for i from 1 to N do
	eqs :=eqs union {add(subs(n=j,summand),j=1..i) = formOfSolution(i)};
od:
#print(eqs,vars):
sols := solve(eqs,vars);
if((subs(sols,vars)={0}) or (sols=NULL)) then
return(FAIL):
fi:
#print(formOfSolution):
#print(subs(sols,formOfSolution));
#return(subs(sols,formOfSolution)):
return(sols):
end:

hyperGeoStep:=proc(rationalForm,c,n,k)
if(k=0) then
return(c);
fi:
return(subs(n=k,rationalForm)*hyperGeoStep(rationalForm,c,n,k-1)):
end:

hyperGeoForm:=(rationalForm,c,n)->k->hyperGeoStep(rationalForm,c,n,k):

multiBasicForm:=(multiRatForm,c,n)->k->multiBasicStep(multiRatForm,c,n,k):

multiBasicStep:=proc(multiRatForm,c,n,k)
if(k=0) then
return(c);
fi:
return(subs(n=k,multiRatForm)*multiBasicStep(multiRatForm,c,n,k-1)):
end: