#OK to post homework
#Yuxuan Yang, March 13th, Assignment 15
with(combinat):

#1
DecDigits:=proc(a,K) local i,dig,num:
num:=a:
dig:=[]:
for i from 1 to K do
dig:=[op(dig),floor(num)]:
num:=10*(num-floor(num)):
od:
dig:
end:

#2
WrongDigits:=proc(a,b,K) local i,A,B,dig:
A:=DecDigits(a,K):
B:=DecDigits(b,K):
dig:=[]:
for i from 1 to K do
 if A[i]<>B[i] then
  dig:=[op(dig),i]:
 fi:
od:
dig:
end:


#Digits:=1000:
#WrongDigits(Pi, NorthM(5000000), 100)
#[8, 22, 35, 48, 49, 50, 61, 62, 63, 73, 75, 77, 78, 85, 86, 87, 89, 90, 91, 98, 100]
#seems not in OEIS

#4
Archimedes:=proc(k,N) local n:
n:=k*2^N:
evalf([n*sin(Pi/n),n*tan(Pi/n)]):
end:
#Archimedes(3, 5)
#[3.141031951, 3.142714599]

#5
#BestMachin:=proc(a) local best:
#we need to run over all integer k? No idea.

#6
#time(PiApc(2000))   is 4.328
#time(PiAp(200))     is 2.375 and 2000 is hopeless.
#Ac(n) is much faster than A(n) for n=2000


#7
#have done this last year, the previous proc is here
PiJG:=proc(dig) local i,j:
Digits:=3*dig+10:
i:=13:
for j from 1 to dig do
i:=evalf(i+(-1)^j*binomial(2*j,j)^5*(820*j^2+180*j+13)/2^(20*j)): 
od:
Digits:=dig:
evalf(sqrt(128/i)):
end:
#PiJG(10000) gives the same output as Digits:=10000:evalf(Pi): on my computer.


#the proc for this year
Jesus:=proc(n) local j,i:
i:=13:
for j from 1 to n-1 do
i:=evalf(i+(-1)^j*binomial(2*j,j)^5*(820*j^2+180*j+13)/2^(20*j)): 
od:
i:
end:
#evalf(Jesus(20) - 128/Pi^2)
#-6.38555069218951279826720567836398509897*10^(-60)



#Maple code for Lecture 15
Help15:=proc(): print(` A(n), PiAp(n), NorthM(N), ArcTanT(x,k), ATAN(L) , FindMcahin(a,k), Ac(n), PiApc(n) `):end:


#A(n) The integral of (x*(1-x))^n/(1+x^2) from x=0 to x=1
A:=proc(n) local x:    int(x^n*(1 - x)^n/(x^2 + 1), x = 0 .. 1): end :

#PiAp(n): The approximation of Pi by a rational number obtained by pretending that A(4*n) above is 0.
PiAp:=proc(n)    local a:    a := A(4*n):    -coeff(a, Pi, 0)/coeff(a, Pi, 1):end:


#ArcTanT(x,k): The truncation of the Taylor series of arctan(x) after k terms, followed by the rigorous error bound
#x^(2*k+1)/(2*k+1)
ArcTanT:=proc(x,k) local i:
[4*add((-1)^i*x^(2*i+1)/(2*i+1),i=0..k), 4*x^(2*k+3)/(2*k+3)]:
end:

#ATAN(L): Inputs a list of numbers (or symbols) L finds the quantity A such that
#arctan(A)= arctan(L[1])+...+arctan(L[k])
#
#arctan(a)+arctan(b)=arctan((a+b)/(1-a*b)):
ATAN:=proc(L) local a,b:
if nops(L)=1 then
 RETURN(L[1]):
fi:
a:=ATAN([op(1..nops(L)-1,L)]):
b:=L[nops(L)]:
normal((a+b)/(1-a*b)):
end:


#FindMachin(a,k): finds the unique b such that ATAN([a$k,b])=1
FindMachin:=proc(a,k) local b:
solve(ATAN([a$k,b])=1,b):
end:


#NorthM(N): The truncated Gregory-Leibnitz formula at N, giving most correct digits of Pi when N=10^k/2, e.g. NorthM(500000); It was discovered emprically, in 1988, by
#at the time undergraudate R.D. North
NorthM:=proc(N) local k:
evalf(4*add((-1)^(k-1)/(2*k-1),k=1..N)):
end:



##ADDED AFTER CLASS
#Ac(n): Same as A(n) but using the third-order recurrence gotten from the amazing Almkvist-Zeilberger algorithm
#1/2*(45*n^2-88*n+32)/(2*n-1)/(5*n-7)/N-(25*n^2-50*n+18)/(2*n-1)/(5*n-7)/N^2+(n-2)*(5*n-2)/(2*n-1)/(5*n-7)/N^3
Ac:=proc(n) local L,n1,newguy:option remember:
if n=0 then
Pi/4:
elif n=1 then
-1+1/2*ln(2)+1/4*Pi:
elif n=2 then
 -2/3+ln(2):
else
L:=[Pi/4,-1+1/2*ln(2)+1/4*Pi, -2/3+ln(2)]:

for n1 from 3 to n do

newguy:=1/2*(45*n1^2-88*n1+32)/(2*n1-1)/(5*n1-7)*L[3]-(25*n1^2-50*n1+18)/(2*n1-1)/(5*n1-7)*L[2]+(n1-2)*(5*n1-2)/(2*n1-1)/(5*n1-7)*L[1]:
L:=[L[2],L[3],newguy]:
od:
RETURN(L[3]):
fi:

end:





#PiApc(n): The approximation of Pi by a rational number obtained by pretending that Ac(4*n) above is 0.
PiApc:=proc(n)    local a:    a := Ac(4*n):    -coeff(a, Pi, 0)/coeff(a, Pi, 1):end: