#Homework 15 by Ben Miles for Experimental math OK to post
#Number 1
An:=proc(n) local i,c: option remember: with(combinat):
c:=0:
if n=0 then
 1:
elif n=1 then
 1:
elif n=2 then
 1:
else
 for i from 1 to (n)/2 do:
  c:=c+binomial(n-1,2*i-1)*An(2*i-1)*An(n-2*i):
 od:
 RETURN(c):
fi:
end:
#this gets nowhere near pi, but very close to 2/pi. In fact evalf(2/evalf((An(200)/(200*An(199))))-Pi);
# returns something on the order of 10^-94




#Number 2
MaxDev:=proc(a,d,w,N) local P,c,i,j,w1:
Digits:=N*d^(nops(w))*20:
P:=RtoS(evalf(a),d,N*d^(nops(w))):
c:=0:
for i from 1 to nops(P)-nops(w)+1 do
w1:=[seq(P[i+j],j=0..nops(w)-1)]:
if w1=w then
 c:=c+1:
fi:
od:
c:

end:




#Number 3
J:=proc(n) local t:
int(t^(4*n)*(1-t)^(4*n),t=0..1);
end:

First40J:=proc() local i,L:
L:=[seq(J(i),i=1..40)];
GuessH(L,n,N);
end:

Jfast:=proc(M) local H:
H:=First40J():
evalf(HtoSeq(H,n,N,[1/630],M)[M]);
end:

timecompare:=proc(M):
time(Jfast(M))/time(J(M)):
end:








#old stuff
Help15:=proc(): print(`IsUD(pi),AnStupid(n), rf(a,k) ,SR1(M) `): 
print(` RtoS(a,d,k), RtoS(a,d,k) `):  end:


 with(combinat):

  #IsUD(pi): Is the permutation pi Up-Down
  #For example IsUD([1,3,2]); is true 
   #IsUD([1,2,3]); is false
  IsUD:=proc(pi) local i:
  evalb(
  {seq( evalb(pi[2*i-1]<pi[2*i]), i=1..nops(pi)/2),
   seq( evalb(pi[2*i]>pi[2*i+1]), i=1..(nops(pi)-1)/2)}=
     {true}):
   end:

   #AnStupid(n): the number of up-down permutations of {1, ...,n} 
   AnStupid:=proc(n) local co, S, pi:
   S:=convert(permute(n),set):

   co:=0:

   for pi in S do
     if IsUD(pi) then
         co:=co+1:
     fi:
  od:

  co:
  end:



 #rf(a,k):=a(a+1)(a+2)...*(a+k-1)
 rf:=proc(a,k) local i:
  mul(a+i,i=0..k-1):
 end:

 #using M terms in Ramanujan's slow formula for 2/Pi
 SR1:=proc(M) local k:
 2/evalf(add((-1)^k*(4*k+1)*rf(1/2,k)^3/k!^3,k=0..M)):
 end:


 #RtoS(a,d,k): the list of the first k digits in the base-d
 #rep. of the real number a between 0 and 1
 RtoS:=proc(a,d,k) local a1,L,kirsten,i:

  a1:=evalf(a):
  if a1>=1 or a1<=0 then
    RETURN(FAIL):
  fi:

 L:=[]:
  
  
 for i from 1 to k do
  a1:=a1*d:
   kirsten:=trunc(a1):
  L:=[op(L),kirsten]:
   a1:=a1-kirsten:
 od:


  if add(L[i]/d^i,i=1..k)-evalf(a)>1/d^(k-1) then

    RETURN(FAIL):

   else
    RETURN(L):

   fi:

end:






   









   


 with(combinat):

  #IsUD(pi): Is the permutation pi Up-Down
  #For example IsUD([1,3,2]); is true 
   #IsUD([1,2,3]); is false
  IsUD:=proc(pi) local i:
  evalb(
  {seq( evalb(pi[2*i-1]<pi[2*i]), i=1..nops(pi)/2),
   seq( evalb(pi[2*i]>pi[2*i+1]), i=1..(nops(pi)-1)/2)}=
     {true}):
   end:

   #AnStupid(n): the number of up-down permutations of {1, ...,n} 
   AnStupid:=proc(n) local co, S, pi:
   S:=convert(permute(n),set):

   co:=0:

   for pi in S do
     if IsUD(pi) then
         co:=co+1:
     fi:
  od:

  co:
  end:
#Number 1
EZlist:=proc(F,n,k,M) local n1,k1:
eval([seq(sum(subs({n=n1,k=k1},F),k1=0..n1),n1=1..M)]):
end:



EmpiricalZeilberger:=proc(F,n,k,N,M) local L:
L:=EZlist(F,n,k,M):
GuessH(L,n,N):
end:

#Number 2
#EmpiricalZeilberger(binomial(n,k)^4,n,k,N,55);
#EmpiricalZeilberger(binomial(n,k)^2*binomial(n+k,k)^2,n,k,N,55);
#EmpiricalZeilberger(binomial(n,k)^3*binomial(n+k,k)^2,n,k,N,200);

#Number 3
#The first two are the same. The third one is different,

#Number 4
FirstMterms:=proc(F,n,k,M1,M2) local N,P,a,L,L1,i;
P:=EmpiricalZeilberger(F,n,k,N,M1):
L:=EZlist(F,n,k,M1):
a:=degree(P,N):
L1:=[seq(L[i],i=1..a)]:
HtoSeq(P,n,N,L1,M2):
end:

#Number 5
qGuessH:=proc(L,n,N,q):

end:
Helpold:=proc(): print( ` GuessH(L,n,N) , SpellOut(ope,n,N,f)`): 

print(`HtoSeq(ope,n,N,Ini,M)`):
end:
#old stuff from C13.txt
Help13:=proc(): print(` LIF(P,u,N), GuessH1(L,ORD,DEG,n,N)`): end:
 #Lagrange Inversion applied to tree-enumeration
   #of ordered
 #Starting the holonomic ansatz (a.k.a P-recursive)
 #LIF: Inputs an expression P of the variable 
 #u that possesses
 #a Maclaurin expansion in u, and outputs the
  #first N terms of the series expansion of the
  #unique f.p.s u(t) satisfying the functional eq.
  #u(t)=t*P(u(t)). For example,
  #LIF(1+u^2,u,10); should return something with
 #Catalan numbers
 LIF:=proc(P,u,N) local n:
 [seq(coeff(taylor(P^n, u=0,n),u,n-1)/n,n=1..N)]:
  end:
  #GuessH1(L,ORD,DEG,n,N): Inputs a sequence of numbers L
  #and pos. integer ORD and non-neg. integer DEG and letters
  #n and N and outputs an linear recurrence operator with
  #poly. coeffs. (conj.) annihilating the sequence L
  #nops(L)>=(1+ORD)*(1+DEG)+ORD+6  
  GuessH1:=proc(L,ORD,DEG,n,N) local ope,a,i,j,var,eq,n1,var1:
  if nops(L)<(1+ORD)*(1+DEG)+ORD+6 then
   print(`We need at least`, (1+ORD)*(1+DEG)+ORD+6, `terms in L`):
   RETURN(FAIL): 
  fi:
  ope:=add(add(a[i,j]*n^i*N^j,i=0..DEG),j=0..ORD):
  var:={seq(seq(a[i,j],i=0..DEG),j=0..ORD)}:
   eq:={seq(add(subs(n=n1, coeff(ope,N,j))*L[n1+j],j=0..ORD),
    n1=1..(1+DEG)*(1+ORD)+6)}:
  var1:=solve(eq,var):
  ope:=subs(var1,ope):
   if ope=0 then
     RETURN(FAIL):
    fi:
  ope:=subs({seq(v=1, v in var)}, ope):
  ope:=add(factor(coeff(ope,N,j))*N^j, j=0..degree(ope,N)):
eq:=expand({seq(add(subs(n=n1, coeff(ope,N,j))*L[n1+j],j=0..ORD),
    n1=1..nops(L)-ORD)}):
  if eq<>{0} then
   print(ope, `didn't work out`):
    RETURN(FAIL):
   else
     RETURN(ope):
   fi:
   end:
#end old stuff
 #start new stuff
 #GuessH(L,n,N): Inputs a sequence of numbers L
 # and symbols n and N
  #and finds a minimal (in the sense of ord+deg)
  # linear recurrence operator with
  #poly. coeffs. (conj.) annihilating the sequence L with
  #nops(L)>=(1+ORD)*(1+DEG)+ORD+6  
  GuessH:=proc(L,n,N) local K,ope,DEG,ORD:
  ORD:=1:
   DEG:=0:
  for K from 1 while nops(L)>=(1+ORD)*(1+DEG)+ORD+6 do
    for  ORD from 1 to K do
         DEG:=K-ORD:
         ope:=GuessH1(L,ORD,DEG,n,N):
         if ope<>FAIL then
            RETURN(ope):
          fi:
 od:
 od:
  FAIL:
   end:   
#SpellOut(ope,n,N,f)
#writes the fact ope(n,N)f(n)=0 in common languae 
SpellOut:=proc(ope,n,N,f) local i:
add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0:
end:
#HtoSeq(ope,n,N,Ini,M): Inputs a linear recurrence operator
#ope phrased in terms of the discrete variable n and
#its corresponding shift operator N, and ORD=degree(ope,N)
#values for the sequence outputs the first M terms of
#the sequence. For example,
#HtoSeq(n+1-N,n,N,[1],4); should yield
#[1,2,6,24]
#N^(-ORD)*ope(N,n)
HtoSeq:=proc(ope,n,N,Ini,M)  local L,n1,ORD,BenMiles,kirsten,i:
ORD:=degree(ope,N):
if nops(Ini)<>ORD then
  RETURN(FAIL):
fi:
BenMiles:=expand(subs(n=n-ORD,ope)/N^ORD):
L:=Ini:
for n1 from nops(Ini)+1 to M do
kirsten:=subs(n=n1,coeff(BenMiles,N,0)):
if kirsten=0 then
print(`blows up at`, n1):
print(`so far we have`, L):
  RETURN(FAIL):
fi:
L:=[op(L),
-add(subs(n=n1,coeff(BenMiles,N,-i))*L[n1-i],i=1..ORD)/kirsten]:
od:
L:
end: