# Iterate, print and plot any discrete map 
# accelerated by Shanks recursion of Aitken's process
N = 322;	M = 9	# length of sequence, number of Aitken processes + 1
import numpy as np
import matplotlib.pyplot as plt
plt.figure()
z = np.zeros((N)).astype(complex)   	# sequence of z_n
A = np.zeros((N+1,M)).astype(complex)  # accelerated series
z_inf = 0.43828293672703215+0.3605924718713855j

z0 = 1j     # imaginary unit is the starting value 
z[0] = z0; A[0,0] = z0

for i in range(1,N):  	# mapping is given by:
    z[i] = z0**z[i-1]  # complex exponentiation
    A[i,0] = z[i]

# compute A_n, B_n, C_n,...  (k = 1,2,3,..)
for k in range(1,M):	
    p = k-1   # 0,1,2,..
    max_n = 55+ 120/k
    for i in range(2*k,N):  # start filling from 2,4,6,..
        denom = (A[i,p]-A[i-1,p])-(A[i-1,p]-A[i-2,p])
        if(abs(denom) < 1e-15): print("denom warning ",k,i)
        # + 1e-30 safeguards against division by zero
        A[i,k] = A[i,p] - (A[i,p]-A[i-1,p])**2/(1e-30+denom)
        if(i > max_n):  A[i,k] = z_inf + 0.8e-16

# plot & print
for k in range(M):
    r = range(2*k,N)
    zz = A[2*k:N,k]    # k-times Aitken-processed z-sequence
    dist = abs(zz-z_inf)
    plt.plot(r,np.log10(dist))
    print(" k = ",k,"   # elem ",N-2*k)
    print(len(zz),"\n",zz)
    for i in range(N-2*k):
        print(i," ",zz[i])

#plt.scatter(zz.real,zz.imag) 
#plt.plot(zz.real,zz.imag)
# finish
plt.xlabel("iteration number")
plt.ylabel("log10(error)   =   log10 | A(k,n) - z* |")
plt.title("...i^i^i  sequence Aitken, Shanks-accelerated") 
plt.show()