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