#Pi Day contest
HelpJG:=proc(): print(`JesusGah(N), JesusGjb(N), JesusGjt(N), JesusJGjk(N), JesusJGkw(N),  JesusJGlm(N), JesusGme(N), JesusGpb(N),  JesusGsm(N), JesusGtk(N)  `): end:
Help:=proc(): print(` ChudnovskyJB(N),  ChudnovskyJT(N),  ChudnovskyJK(N), ChudnovskyKW(N), ChudnovskyLM(N),  ChudnovskyME(N), NB(), PiCompPB(N), estPiSM(N) `): 
print(`ChudnovskyTK(N), VietesPiTK(N), BBPtk(N) `):
end:

#Start Aurora Hively
JesusGah := proc(N) local n:
print(`By Aurora Hively`):
print(``):
evalf(add( (-1)^n * binomial(2*n,n)^5 * (820*n^2 + 180*n + 13)/(2^(20*n)) , n=0..N)):
end:

# abs(JesusG(10) - 128/Pi^2);
# output: 1E-32
# abs(JesusG(100) - 128/Pi^2);
# output: 0.
#End Aurora Hively

#Start James Betti
JesusGjb := proc(N) local total,i;
print(`By James Betti`):
    total := 0;
    for i from 0 to N do
        total := total+(-1)^i*((2*i)!/(2^i*i!)^2)^5*(20*i*(9+41*i)+13)/2^(10*i);
    od;
    evalf(sqrt(128/total));
end:

# https://en.wikipedia.org/wiki/Chudnovsky_algorithm
# Chudnovsky algorithm using binary splitting as in the above article.
ChudnovskyJB := proc(N) local split,P,Q,R;
print(`By James Betti`):
    split := proc(a,b) local P1,Q1,R1,P2,Q2,R2;
        P,Q,R;
        if b=a+1 then
            P := -(6*a-1)*(2*a-1)*(6*a-5);
            Q := 10939058860032000*a^3;
            R := P*(545140134*a+13591409);
        else
            P1,Q1,R1 := split(a,floor((a+b)/2));
            P2,Q2,R2 := split(floor((a+b)/2),b);
            P := P1*P2;
            Q := Q1*Q2;
            R := Q2*R1+P1*R2;
        fi;
        evalf(P),evalf(Q),evalf(R);
    end;
    P,Q,R := split(1,N);
    evalf(426880*sqrt(10005)*Q/(13591409*Q+R));
end:

#End James Betti

#Start Jeff Tang
JesusGjt := proc(N) local S, n:
print(`By Jeff Tang`):
    S := evalf(sum((-1)^n * binomial(2*n, n)^5 * (820*n^2 + 180*n + 13) / 2^(20*n), n = 0 .. N),500):
end:

ChudnovskyJT := proc(n) local i, numerator, denominator, term, sum_value, pi_approx:
print(`By Jeff Tang`):
    sum_value := 0:
    
    for i from 0 to n do
        numerator := factorial(6 * i) * (545140134 * i + 13591409);
        denominator := factorial(3 * i) * (640320^(3 * i + 3/2));
        term := (-1)^i * numerator / denominator;
        
        sum_value := sum_value + term;
    od;
    
    pi_approx := 12 * evalf(sum_value);
    
    evalf(1 / pi_approx);
end:

#End Jeff Tang

#Start Joseph Koutsoutis
JesusGjk := proc(N) local n:
print(`By Joseph Koutsoutis`):
  add((-1)^n * binomial(2*n, n)^5 * (820*n^2 + 180*n + 13) / 2^(20*n), n=0..N):
end:

#This is a Maple implementation of the Chudnovsky algorithm described here: 
#https://en.wikipedia.org/wiki/Chudnovsky_algorithm

BSjk:=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 :=BSjk(a, m):
    Pmb, Qmb, Rmb :=BSjk(m, b):
        
    Pab := Pam * Pmb:
    Qab := Qam * Qmb:
    Rab := Qmb * Ram + Pam * Rmb:
  fi:
 Pab, Qab, Rab:
end:


ChudnovskyJK := proc(n) local P1n, Q1n, R1n:
print(`By Joseph Koutsoutis`):
    P1n, Q1n, R1n :=BSjk(1, n):
    426880 * sqrt(10005) * Q1n / (13591409*Q1n + R1n):
end:

#end Joseph Koutsoutis


#Start Kaylee Weatherspoon
JesusGkw:=proc(N) local T,n:
print(`By Kaylee Weatherspoon`):
T:=sum(((-1)^n)*binomial(2*n, n)^5*(820*n^2+180*n+13)/(2^(20*n)), n=0..N):   
evalf(T);
end: 


#want roughly 12.969

#Even for N as small as 4, this works really well! I got 12.96911151, and the true value (truncated appropriately) is 12.96911150. 

#2:

#used Chudnovsky algorithm wikipedia article

insideKW:=proc(j)
return(evalf((-(6*j-1)*(2*j-1)*(6*j-5))/(10939058860032000*j^3))):
end:

ChudnovskyKW:=proc(N) local T, Acoeff, Baddend, Cprod, D, k, prodfactor:
print(`By Kaylee Weatherspoon`):
Acoeff:=1/(426880*evalf(sqrt(10005))):
Baddend:=13591409:
prodfactor:=1:
#splitting up as advised_______
for k from 1 to N do
 prodfactor:=prodfactor*evalf(insideKW(k)):
od:

Cprod:=545140134*N+13591409:
D:=Acoeff*(Baddend+prodfactor*Cprod):
T:=1/D:
return(evalf(T)):

end:

#End Kaylee Weatherspoon



# Start Lucy Martinez, 03-12-2025, Assignment 14 Pi

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
binarySplitLM:=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):=binarySplitLM(a, m):
 (Pmb, Qmb, Rmb):=binarySplitLM(m, b):
 Pab:=Pam*Pmb:
 Qab:=Qam*Qmb:
 Rab:=Qmb*Ram+Pam*Rmb:
fi:
Pab,Qab,Rab:
end:

#ChudnovskyLM(N): Computes the first N digits of Pi using binary splitting
# for speeding up numerical evaluation by using the Chudnovsky algorithm
ChudnovskyLM:=proc(N) local P1n, Q1n, R1n, sq:
(P1n,Q1n,R1n):=binarySplitLM(1,N):
sq:=evalf(sqrt(10005)):
Digits:=N:
evalf((426880*sq*Q1n)/(13591409*Q1n+R1n)):
end:
# End Lucy Martinez, 03-12-2025, Assignment 14 Pi

# Start Matthew Esaia

#JesusG(n): uses N terms of Jesus Guillera's formula to compute an
#approximation to 128/pi^2
JesusG:=proc(N) local n:
print(`By Matthew Esaia`):
	  evalf(add((-1)^n * (binomial(2*n, n))^5 * (820*n^2 + 180*n + 13) / (2^(20*n)), n=0..N), 500):
end:
# JesusG(10) has an error of 1.048 x 10^-32
# JesusG(100) has an error of 4.23 x 10^-304

#RealFrac(): true value of 128/pi^2
RealFrac:=proc():
	return evalf(128 / Pi^2, 500):
end:

## PROBLEM 2
ChudnovskyME:=proc(E) local i,err,sum,const:
print(`By Matthew Esaia`):
	  sum := 0:
	  err := evalf(Pi - sum):
	  const := 1 / (426880 * sqrt(10005)):
	  for i from 0 while err > 10^(-E) do
	      sum := sum + ( (-1)^i * (6*i)! * (545140134*i + 13591409) / ((3*i)! * (i!)^3 * (640320)^(3*i))):
	      	  err := evalf(abs(Pi - (1 /(const*sum))), E + 10):
		      #print(i, evalf(1 /(const*sum), E)):
		      od:
		      return evalf(1 /(const*sum), E):
end:
# End Matthew Esaia

#Start Nick Belov

JesusG:=proc(N)
sum((-1)^n*binomial(2*n,n)^5*(820*n^2+180*n+13)/(2^(20*n)), n=0..N):
end:
#eval(abs(JesusG(10) - 128/pi^2))=10^-8
#JesusG(eval(abs(JesusG(100) - 128/pi^2))) was also 10^8
#infact Jesus of everything after 5 was 10^-8 so i imagine thats the point
#where the accuracy exceeds the capability of floating point arithmatic


#pi approximation:

NB:=proc():evalf((🥧:=3+sin(3))+sin(🥧)); : end:

#I found the idea for this here https://codegolf.stackexchange.com/questions/22009/pi-calculation-code-golf
#I hard coded two iterations of newtons method in as few characters possible.
#End  Nick Belov

#start  Pablo Blanco
JesusGpb:=proc(N) local i,j:
 13-add((mul(2*j+1,j=0..i-1)/mul(2*j,j=1..i))^5 * (20*i*(50+41*(i-1)) +13)/(2^(10*i))  ,i=1..N)
end:


# Comparisons:
# JesusG(10):  0.00005167 from true value
# JesusG(100): 0.00005167 ... from true value ?


# opting out of this one
PiCompPB:=proc(N):
 evalf(sqrt(128/JesusG(N)),10*N^2):
end:


# Remark: Since I'm sure Jesus's computation should improve... I assume it exceeds Maple rounding limits?
#End  Pablo Blanco

#Start Salman Manzoor 


JesusGsm:=proc(N) local n,i:
print(`By Salman Manzoor`):
if evalb(N=1) then RETURN(13): fi:
evalf(13 + add(    (-1)^n * (product(2*i-1,i=1..n)/product(2*i-1,i=1..n) )^5 *
( 20*n*(9+41*n) +13 ) / 2^(10*n)  , n=1..N-1  )) :
end: 

#Plouffe (2006) formula exactly equal to pi
estPiSM:=proc(n) local i:
print(`By Salman Manzoor`):
evalf(add( (1/i)*(72/(fooSM(i))-96/(fooSM(2*i))+24/(fooSM(4*i))),i=1..n)):
end:

# Andrew Fraker (2014) estimates e^pi around (5^4/3^3)
# we approximate foo_true(i) = 1/(e^(Pi*x)-1)

fooSM:=proc(x)
evalf(1/((5^4/3^3)^x-1)):
end:

#End  Salman Manzoor 

#start Tijil Kiran
JesusGtk := proc(N)
    local n, term, sum;
 print(` By Tijil Kiran`):
    sum := 0;
    term := 1;
    for n from 0 to N-1 do
        term := (binomial(2*n, n)^3 * (42*n + 5)) / (2^(12*n + 4));
        sum := sum + term;
    end do;
    return sum;
end proc:

ChudnovskyTK := proc(N)
    local n, term, sum;
 print(` By Tijil Kiran`):
    sum := 0;
    term := 1;
    for n from 0 to N-1 do
        term := (factorial(6*n) / (factorial(n)^3 * factorial(3*n))) * (13591409 + 545140134*n) / (640320^(3*n + 3/2));
        sum := sum + term;
    end do;
    return 1 / (sum * 12);
end proc:

#VietesPiTK(N)
VietesPiTK := proc(N)
    local k, term, product;
 print(` By Tijil Kiran`):
    product := 1;
    term := sqrt(2);
    for k from 1 to N do
        product := product * (2 / term);
        term := sqrt(2 + term);
    end do;
    return 2*product;
end proc:

#BBPtk(N)
BBPtk := proc(N)
    local k, sum;
 print(` By Tijil Kiran`):
    sum := 0;
    for k from 0 to N-1 do
        sum := sum + (1 / 16^k) * (
            4 / (8*k + 1) -
            2 / (8*k + 4) -
            1 / (8*k + 5) -
            1 / (8*k + 6)
        );
    end do;
    return evalf(sum);
end proc:

#End Tijil Kiran