# OK To Post
# Homework 7 — Jike Liu

# ------------------------------------------------------------
# (0) Read C7.txt
# ------------------------------------------------------------
# read `C7.txt`:


# ------------------------------------------------------------
# Recall
# ------------------------------------------------------------

# RF(a,n): rising factorial (Pochhammer symbol) (a)_n = a(a+1)...(a+n-1)
RF := proc(a, n) local i;
    mul(a + i, i = 0 .. n - 1);
end proc:

# HowManyDigits(f,N): number of decimal digits of agreement with Pi for f(N)
# NOTE: Digits is Maple's global working-precision variable.
HowManyDigits := proc(f, N)
    Digits := 10000;
    -trunc(log[10](abs(evalf(f(N) - Pi))));
end:


# ------------------------------------------------------------
# (1) Q1: How many terms in JG1 are needed for >=100 digits?
# ------------------------------------------------------------

# I found that truncating JG1 at n=33 gives >= 100 digits:
# So 34 terms are needed

# result:
# HowManyDigits(JG1,33)=102.
# HowManyDigits(JG1,32)=99.

# ------------------------------------------------------------
# (2) Q2: Five chosen series
# ------------------------------------------------------------
# I chose the five formulas.
# Each JGjlX(N) returns an approximation to Pi using terms n=0..N.

# ------------------------------------------------------------
# (2.1) JGjl1
# Sum_{n>=0} (1/4)_n (1/2)_n (3/4)_n / (n!)^3 * (40n+3) * 1/7^(4n) = 49*sqrt(3)/(9*pi)
# Therefore pi = (49*sqrt(3)) / (9 * sum)
JGjl1 := proc(N) local n, a;
    a := add(evalf(RF(1/4, n)*RF(1/2, n)*RF(3/4, n)*(40*n + 3)
                   /(RF(1, n)^3*7^(4*n))),
             n = 0 .. N);
    evalf(49/9*sqrt(3)/a);
end:

# ------------------------------------------------------------
# (2.2) JGjl2
# Sum_{n>=0} (-1)^n (6n)!/(n!)^6 * (5418 n^2 + 693 n + 29) / 2880^(3n) = 128*sqrt(5)/pi^2
# Therefore pi = sqrt( 128*sqrt(5) / sum )
JGjl2 := proc(N) local n, a;
    a := add(evalf((-1)^n*RF(1, 6*n)*(5418*n^2 + 693*n + 29)
                   /(RF(1, n)^6*2880^(3*n))),
             n = 0 .. N);
    evalf(sqrt(128*sqrt(5)/a));
end:

# ------------------------------------------------------------
# (2.3) JGjl3
# Sum_{n>=0} (6n)!/(n!)^6 * (532 n^2 + 126 n + 9) / 10^(6n) = 375/(4*pi^2)
# Therefore pi = sqrt( (375/4) / sum )
JGjl3 := proc(N) local n, a;
    a := add(evalf(RF(1, 6*n)*(532*n^2 + 126*n + 9)
                   /(RF(1, n)^6*10^(6*n))),
             n = 0 .. N);
    evalf(sqrt(375/4*1/a));
end:

# ------------------------------------------------------------
# (2.4) JGjl4
# Sum_{n>=0} (-1)^n (3n)!(2n)!/(n!)^5 * (14151 n + 827) / 300^(3n) = 1500*sqrt(3)/pi
# Therefore pi = (1500*sqrt(3)) / sum
JGjl4 := proc(N) local n, a;
    a := add(evalf((-1)^n*RF(1, 3*n)*RF(1, 2*n)*(14151*n + 827)
                   /(RF(1, n)^5*300^(3*n))),
             n = 0 .. N);
    evalf(1500*sqrt(3)/a);
end:

# ------------------------------------------------------------
# (2.5) JGjl5
# Sum_{n>=0} ( (1/2)_n^3 (1/4)_n (3/4)_n / (n!)^5 )*(1/16)^n*(120 n^2 + 34 n + 3) = 32/pi^2
# Therefore pi = sqrt( 32 / sum )
JGjl5 := proc(N) local n, a;
    a := add(evalf(RF(1/2, n)^3*RF(1/4, n)*RF(3/4, n)*(1/16)^n*(120*n^2 + 34*n + 3)
                   /RF(1, n)^5),
             n = 0 .. N);
    evalf(sqrt(32/a));
end:


# ------------------------------------------------------------
# (2.6) Compute digit agreement at N=100 (i.e., 101 terms)
# ------------------------------------------------------------

# For my run, the outputs were:
# HowManyDigits(JGjl1, 100);  -> 341
# HowManyDigits(JGjl2, 100);  -> 576
# HowManyDigits(JGjl3, 100);  -> 135
# HowManyDigits(JGjl4, 100);  -> 545
# HowManyDigits(JGjl5, 100);  -> 122


# End of file.