import matplotlib.pyplot as plt
from matplotlib import interactive
from numpy import pi,log10
interactive(True)
N = int(input("N terms = "))

#       prepare plot frame and description
plt.figure(figsize=(6,5))
plt.xlabel('log$_{10}$ N',fontsize=14)
plt.ylabel('log$_{10}$ |$x-\pi$|',fontsize=14)
plt.title('Accelerated convergence of Leibniz ser. for $\pi$',fontsize=14)
plt.axis([0.,4.75, -16.,0.])
plt.grid(True,color=(.8,.9,.97))

'''
    Series x = 4/1 -4/3 +4/5 -4/7 +...+(+-1)/(odd N) due to Leibniz
    is accelerated by repeated averaging of neighboring terms.
    This is done with almost zero storage and one series summation
    
    Shanks transformation is also done, turns out equivalent
    to the twice-averaged series. 
    Shanks y_n = x_n - (x_{n} - x_{n-1})**2 / (x_n -2*x_{n-1} +x_{n-2})
'''
plt.figure(figsize=[8,6])
plt.axis(limits=[0,5,-16,0])
x = x_1 = a = b = c = d = const = 4. 
i = 1
while(i < N):    # c is the number of people/obj.
    i += 2
    const = -const
    a_1 = a
    b_1 = b
    c_1 = c
    d_1 = d
    x_2 = x_1
    x_1 = x
    x = x + const/i
    a = (x+x_1)*0.5
    b = (a+a_1)*0.5
    c = (b+b_1)*0.5
    d = (c+c_1)*0.5
    f = x - (x-x_1)**2/(x-2*x_1+x_2)
    # print("i",i,"dx  ",x-pi,"   da,db ",a-pi,b-pi)
    if(i<900 or (i>900 and (i-1)%32==0 and (i-1)%128==0)):
        i0 = log10(i)
        y0 = log10(max(abs(x-pi),1e-18))
        e1 = log10(max(abs(a-pi),1e-18))
        e2 = log10(max(abs(b-pi),1e-18))
        e3 = log10(max(abs(c-pi),1e-18))
        e4 = log10(max(abs(d-pi),1e-18))
        e5 = log10(max(abs(0.5*(d+d_1)-pi),1e-18))
        fl = log10(max(abs(f-pi),1e-18))
    # plot singe point   
        plt.scatter(i0,y0,s=5,color=(.3,0.,.8),linewidth=2)
        plt.scatter(i0,e1,s=5,color=(.7,1.,.2),linewidth=2)
        plt.scatter(i0,e2,s=4,color=(0.,1.,0.),linewidth=2) 
        plt.scatter(i0,e3,s=4,color=(.2,.4,1.),linewidth=2)
        plt.scatter(i0,e4,s=5,color=(.9,.3,.0),linewidth=2)
        plt.scatter(i0,e5,s=5,color=(.0,.0,.0),linewidth=2)
        plt.scatter(i0,fl,s=6,color=( 1,.7,.0),linewidth=2)
        eps = 2.22e-16
        log_eps = log10(eps)
        l05 = log10(1/2)
    if((i-1)%1000==0): print(int(i/N*100),'%')
plt.plot([0,5],[log_eps,log_eps],linewidth=1)
plt.plot([0.5,5],[log_eps,log_eps-l05+5/2],color=(0,.7,0),linewidth=0.6)
plt.grid(color=(.8,.8,.8))
plt.text(i0+.1,y0,"Leibniz"); plt.text(i0+.2,y0 -.6,"($4.3\cdot 10^5$ yr)")
plt.text(i0+.1,e1,"1av  (2.3 day)")
plt.text(i0+.1,e2,"2av  (11 min)")
plt.text(i0+.1,e3,"3av  (42 s)")
plt.text(i0+.1,-15.5,"4av (8 s)")
plt.text(1.66,-13,"5av (4 s)")
plt.text(i0+.1,fl-.25,"Shanks (7 min)")
plt.text(1.75,-15.55,"machine $\epsilon$")
plt.text(0.3,-14.6,"$\log_{10}\epsilon\sqrt{N/2}$",color=(0,.7,0))
plt.text(0.3,-13.4,"roundoff error",color=(0,.7,0))
plt.text(2.5,-2,"$\sim N^{-1}$")
plt.text(2.5,-4.5,"$\sim N^{-2}$")
plt.text(2.5,-7,"$\sim N^{-3}$")
plt.text(1.2,-8.25,"$\sim N^{-6}$")
plt.ylabel("logarithmic error")
plt.xlabel("log$_{10}$ N")
plt.title("Leibnitz series for $\pi$ and its seedup by averaging and Shanks method")
print(" rough prediction for N to reach machine eps")
print(" including the PyPlot graphics, which is a slow routine")
print(" Leibniz: ",1.5*(1/eps)/500/3.15e7,' yr')
print(" 1 average",1.5*(1/eps)**(1/2)/500/(24*3600),' d')
print(" 2 averages (Shanks)",2*(1/eps)**(1/3)/500/60,' min')
print(" 3 averages ",2.5*(1/eps)**(1/4)/500/60,' min')
print(" 4 averages ",3*(1/eps)**(1/5)/500,' s')
print(" 5 averages ",5*(1/eps)**(1/6)/500,' s')

ans = input("ok?")