# hw15JikeLiu.txt
# Math 640 — Homework 15
# Jike Liu
#

# ------------------------------------------------------------
# 1) Another sequence converging to Pi
# ------------------------------------------------------------

# a(n) = 16^n/(n*binomial(2*n,n)^2)

a := proc(n)
  16^n/(n*binomial(2*n,n)^2):
end:

# SeqPi(N): first N floating-point terms of a(n)
SeqPi := proc(N) local n;
  [seq(evalf(a(n)), n=1..N)]:
end:

# TestPiConv(N): first N absolute errors |a(n)-Pi|
TestPiConv := proc(N) local n;
  [seq(evalf(abs(a(n)-Pi)), n=1..N)]:
end:

# TestPiConv2(N): triples [n, a(n), |a(n)-Pi|]
TestPiConv2 := proc(N) local n;
  [seq([n, evalf(a(n)), evalf(abs(a(n)-Pi))], n=1..N)]:
end:


# A_pi(N): first N exact terms
A_pi := proc(N) local n;
  [seq(a(n), n=1..N)]:
end:

# DigitsPi(x): number of correct digits of Pi in x
DigitsPi := proc(x) local OLD,res;
  OLD := Digits:
  Digits := 600:
  res := floor(-log10(abs(evalf(x-Pi)))):
  Digits := OLD:
  res:
end:

# Experimental proof in Maple that a(n)->Pi:
# TestPiConv(20);


# 1000-th term and number of correct digits:
# evalf(a(1000),30);
# DigitsPi(a(1000));

# This gives about 3 correct digits of Pi.

# Note:
# PiAppxVC(1000) has length 999, since it is
# [seq(2*i*L[i-1]/L[i],i=2..n)].
# Therefore the last term is PiAppxVC(1000)[999], not [1000].

# DigitsPi(PiAppxVC(1000)[999]);

# Output obtained:
# 476

# Conclusion:
# a(1000) gives about 3 correct digits of Pi,
# while PiAppxVC(1000)[999] gives 476 correct digits of Pi.
# So PiAppxVC converges much faster.


# ------------------------------------------------------------
# 2) Implement K(n) directly and Kc(n) by recurrence
# ------------------------------------------------------------

# K(n): direct integral definition
K := proc(n)
  int(x^(4*n)*(1-x)^(4*n)/(1+x^2), x=0..1):
end:

# Kc(n): same sequence using the second-order recurrence
# K(0)=Pi/4, K(1)=22/7-Pi
Kc := proc(n) local m,A,B,C;
  option remember:

  if n=0 then
    RETURN(Pi/4):
  elif n=1 then
    RETURN(22/7-Pi):
  else
    m := n-2:

    A := 4*(4*m+1)*(4*m+3)*(820*m^3+3993*m^2+6428*m+3420)*(m+1)*(2*m+1):

    B := 26843520*m^7+211258528*m^6+688917624*m^5+1202271190*m^4
         +1207256235*m^3+693651577*m^2+209656416*m+25472340:

    C := 2*(8*m+9)*(8*m+11)*(8*m+13)*(8*m+15)*(820*m^3+1533*m^2+902*m+165):

    RETURN(simplify((A*Kc(n-2)-B*Kc(n-1))/C)):
  fi:

end:

# Compare speeds

t1 := time():
L1 := [seq(K(n), n=1..200)]:
time()-t1;

t2 := time():
L2 := [seq(Kc(n), n=1..200)]:
time()-t2;

# Output obtained:
# [seq(K(n),n=1..200)]   -> 170.016
# [seq(Kc(n),n=1..200)]  -> 0.016


# Therefore [seq(Kc(n),n=1..200)] is much faster.



# End