# OK to post homework
# Vikrant, Mar 13 2022, Assignment 15

# ================================================================================
# 0. Code that has been given.
# ================================================================================

read("C15.txt"):

# ================================================================================
# 1. First k digits in a list.
# ================================================================================

DecDigits:=
    proc(a,K)
        if not type(a,float) then
            return DecDigits(evalf(a,K+1),K):
        fi:
        if a > 1 then
            return DecDigits(a/10,K):
        fi:
        if K = 0 then
            return []:
        fi:
        [trunc(a*10), op(DecDigits(frac(a*10),K-1))]:
    end:

# ================================================================================
# 2. Disagreeing digits up to kth place.
# ================================================================================

WrongDigits:=
    proc(a,b,K)
        local i:
        local A:= DecDigits(a,K):
        local B:= DecDigits(b,K):
        [seq(`if`(A[i]<>B[i],i,NULL),i=1..K)]:
    end:

(*

Digits:=1000:

# Positions in which digits differ:
# [7, 18, 19, 30, 41, 42, 43, 52, 53, 54, ...]
WrongDigits(Pi,NorthM(500000),100);

# Positions in which digits differ:
# [8, 22, 35, 48, 49, 50, 61, 62, 63, 73, ...]
WrongDigits(Pi, NorthM(5000000), 200);

# The positions themselves do not appear on OEIS.
# The actual differences in the runs of differing digits in both executions are [2, -2, 1, -122, 277, ...], which doesn't seem to appear on OEIS either (even unsigned).

*)

GroupDiffPos:=
    proc(L)
        if nops(L)=0 then
            return []:
        fi:
        local i:= 2:
        local j:
        local s:= [0$nops(L)]:
        s[1]:= L[1]:
        while i<=nops(L) and L[i] = L[i-1]+1 do
            s[i]:= L[i]:
            i:= i+1:
        od:
        [[seq(s[j],j=1..i-1)],op(GroupDiffPos([seq(L[j],j=i..nops(L))]))]:
    end:

NumberAtDiff:=
    proc(A,L)
        local l:
        local x:= 0:
        for l in L do
            x:= x*10:
            x:= x+A[l]:
        od:
        x:
    end:

DifferingDigitDiffs:=
    proc(a,b,K)
        local i:
        local A:= DecDigits(a,K):
        local B:= DecDigits(b,K):
        local GPDs:= GroupDiffPos(WrongDigits(a,b,K)):
        [seq(NumberAtDiff(A,i),i in GPDs)] - [seq(NumberAtDiff(B,i),i in GPDs)]:
    end:

# ================================================================================
# 3. <Skipped>
# ================================================================================

# ================================================================================
# 4. Archimedes' Method.
# ================================================================================

Archimedes:=
    proc(k,N)
        if N=0 then
            return [k*sin(Pi/k),k*tan(Pi/k)]:
        fi:
        local pP:= Archimedes(k,N-1):
        local P2:= 2*pP[1]*pP[2]/(pP[1]+pP[2]):
        [sqrt(P2*pP[1]),P2]:
    end:

(*

# [3.14103, 3.14275]
evalf(Archimedes(3,5),6);

*)

# ================================================================================
# 5. Find the best Machin.
# ================================================================================

BestMachin:=
    proc(a)
        local i:
        min[index]([seq(abs(FindMachin(a,i)),i=1..floor(abs(Pi/arctan(a))))]):
    end:

(*

[seq(BestMachin(a),a=2..10)].

# A000012 is the sequence in OEIS.
# I think the sequence rightfully belongs in OEIS, but I don't think using BestMachin as a basis for justifying it being in OEIS is right;
# BestMachin(a >= 2) does not seem a natural way to define 1,1,1,...

*)

# ================================================================================
# 6. PiApc vs PiAp.
# ================================================================================

(*

l:= 1:
# PiApc(2222) takes >5 seconds.
h:= 2222:
while h-l>3 do
    m:= floor((l+h)/2):
    t:= time(PiApc(m)):
    if t>5 then
        h:= m-1:
    else
        l:= m:
    fi:
od:
# PiApc(2078) is roughly the threshold after which PiApc takes >5 seconds.
print(h);

# PiAp(222) already takes >10 seconds. Feels hopeless compared to PiApc(222) which takes roughly 0.05 seconds.

*)

# ================================================================================
# 7. Guillera's formula for 128/Pi^2.
# ================================================================================

Jesus:=
    proc(n)
        local i:
        add((-1)^i*binomial(2*i,i)^5*(820*i^2+180*i+13)/(2^(20*i)),i=0..n):
    end:

(*

# About 6*10^(-63).
evalf(Jesus(20)-128/Pi^2);

*)

# ================================================================================
# 8. <Skipped>
# ================================================================================

# ================================================================================
# 9. <Skipped>
# ================================================================================