# Ok to post homework # Lucy Martinez, 03-12-2025, Assignment 14 Pi Help:=proc(): print(`JesusG(N), Chudnovsky(N)`): end: #Problem 1: Write a procedure JesusG(N) # that uses N terms of Jesus Guillera's amazing formula in his homepage # (given in "Ramanujan notation") that computes an approximation to 128/pi^2 # How close is JesusG(10) to 128/pi^2? # IF Digits:=30, then # Answer: They are off at the 28th digit. The difference is 1.*10^(-28) # How close is JesusG(100) to 128/pi^2? # Answer: They are off at the 28th digit. The difference is 1.*10^(-28) JesusG:=proc(N) local n: add((-1)^n*binomial(2*n,n)^5*(820*n^2+180*n+13)/(2^(20*n)), n=0..N): end: # [Pi Hacktaton]: Use any algorithm you can find # (except "cheating ones" like copying the digits of Pi) to compute Pi. # The fastest and/or most elegant entry will win a Pi T-shirt # [Hint: Many series for Pi are of hypergeometric type: # Sum(p(n)*f(n),n=0..infinity) where p(n) is a polynoial and # f(n)/f(n-1) is a rational function. # It is most efficient to just have the current value of f(n) # and keep adding p(n)*f(n)] # We are going to implement the Chudnovsky algorithm binarySplit:=proc(a, b) local Pab, Qab, Rab, m, Pam, Qam, Ram, Pmb, Qmb, Rmb: if b=a+1 then Pab:=-(6*a-5)*(2*a-1)*(6*a-1): Qab:=10939058860032000*a^3: Rab:=Pab*(545140134*a+13591409): else m:=floor((a + b)/2): (Pam, Qam, Ram):=binarySplit(a, m): (Pmb, Qmb, Rmb):=binarySplit(m, b): Pab:=Pam*Pmb: Qab:=Qam*Qmb: Rab:=Qmb*Ram+Pam*Rmb: fi: Pab,Qab,Rab: end: #Chudnovsky(N): Computes the first N digits of Pi using binary splitting # for speeding up numerical evaluation by using the Chudnovsky algorithm Chudnovsky:=proc(N) local P1n, Q1n, R1n, sq: (P1n,Q1n,R1n):=binarySplit(1,N): sq:=evalf(sqrt(10005)): Digits:=N: evalf((426880*sq*Q1n)/(13591409*Q1n+R1n)): end: