#Nathan Fox #Homework 314159 #I give permission for this file to be posted online ##Read old files read(`C314159.txt`): #Help procedure Help:=proc(): print(` ascfact(k,n) , JG1(N) , JG2(N) , JG3(N) , JG4(N) `): print(` JG5(N) , JG6(N) , g(n) , JGfast(N) , BBP(N) `): print(` BBP2(n) , AndreFF(n) , AndrePi(n) `): end: ##Problem 1 #Ascending factorial, (k)_n in JG's notation ascfact:=proc(k,n) local val, ret, i: ret:=1: val:=k: for i from 1 to n do ret:=ret * val: val:=val + 1: od: ret: end: #Serie 1 on page 22 JG1:=proc(N) local n: 2/add((-1)^n*(ascfact(1/2,n)^3/((n!)^3))*(4*n+1),n=0..N): end: #This converges really slowly #Serie 2 on page 25 JG2:=proc(N) local n: 4/add(1/(2^(2*n))*(ascfact(1/2,n)^3/((n!)^3))*(6*n+1),n=0..N): end: #This gives about half a digit per turn #Serie 3 on page 25 JG3:=proc(N) local n: 2*sqrt(2)/add((-1)^n/(2^(3*n))* (ascfact(1/2,n)^3/((n!)^3))*(6*n+1),n=0..N): end: #This gives about one digit per turn #Serie 4 on page 26 JG4:=proc(N) local n: 8/add((-1)^n/(2^(2*n))* (ascfact(1/2,n)*ascfact(1/4,n)*ascfact(3/4,n)/((n!)^3))* (20*n+3),n=0..N): end: #This gives about one digit per turn #Serie 5 on page 27 JG5:=proc(N) local n: 16/add(1/(2^(6*n))*(ascfact(1/2,n)^3/((n!)^3))*(42*n+5),n=0..N): end: #This gives about two digits per turn #Serie 6 on page 27 JG6:=proc(N) local n: 32*sqrt(2)/add((-1)^n*(3^(3*n))/(2^(9*n))* (ascfact(1/2,n)*ascfact(1/6,n)*ascfact(5/6,n)/((n!)^3))* (154*n+15),n=0..N): end: #This gives about one digit per turn ##Problem 2 #f(n)/f(n-1) in JG formula g:=proc(n): -(6*n-1)*(6*n-2)*(6*n-3)*(6*n-4)*(6*n-5)/(3981312000*n^5): end: #Optimized JG formula JGfast:=proc(N) local ret, tot, n: ret:=29: tot:=1: for n from 1 to N do tot:=tot * g(n): ret:=ret + tot * (5418*n^2+693*n+29): od: sqrt(128*sqrt(5)/ret): end: #JG(1000) took 3.521 seconds #JGfast(1000) took 1.688 seconds ##Problem 3 #BBP formula BBP:=proc(N) local n: add(1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)),n=0..N): end: ##Problem 4 ##BBP digit extraction BBP2:=proc(n) local s1,s2,s3,s4,k,m,ret: if n=0 then return 2: fi: m:=n+10: s1:=add((16&^(n-k) mod (8*k+1))/(8*k+1),k=0..n)+ add(16^(n-k)/(8*k+1),k=n+1..m): s2:=add((16&^(n-k) mod (8*k+4))/(8*k+4),k=0..n)+ add(16^(n-k)/(8*k+4),k=n+1..m): s3:=add((16&^(n-k) mod (8*k+5))/(8*k+5),k=0..n)+ add(16^(n-k)/(8*k+5),k=n+1..m): s4:=add((16&^(n-k) mod (8*k+6))/(8*k+6),k=0..n)+ add(16^(n-k)/(8*k+6),k=n+1..m): ret:=4*s1-2*s2-s3-s4: if ret < 0 then ret:=ret - trunc(ret) + 1: fi: ret:=ret-trunc(ret): ret:=ret * 16: trunc(ret): end: ##Problem 5 #I could not get the formula to work properly on this problem. #Sum((1/16)^i/(8*i+j+1),i=0..100) is larger than 1 and #int(1/(1-x^8)),x=0..1/16) is much smaller than 1 ##Problem 6 #Andre using recurrence AndreFF:=proc(n) local k: option remember: if n=0 then return 1: elif n=1 then return 1: fi: add(binomial(n-1,2*k-1)*AndreFF(2*k-1)*AndreFF(n-2*k), k=1..trunc(n/2)): end: #AndreF(1000) took 31.122 seconds #AndreFF(1000) took 16.716 seconds #I'm not running these for n=4000. It will take a long time. #AndrePi(n): computes 2*n*A[n-1]/A[n]->Pi AndrePi:=proc(n) 2*n*AndreFF(n-1)/AndreFF(n): end: #AndrePi(10)-Pi=7.1209e-5 #AndrePi(100)-Pi=0 (according to Maple) #AndrePi(1000)-Pi=0 (according to Maple)