#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+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+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: