# feigenbaum.py
# http://keithbriggs.info/documents/how-to-calc.pdf
'''
  How to calculate the Feigenbaum constants on your PC, Australian Math. Soc. Gazette 16 (1989). pp.89-92.
  
  This is suitable for low precision only, owing to 
  its slow convergence.

  Feigenbaum constants calculation;
  d = 4.669 201 609 102 990 671 853 203 820 466 ...
  alpha = 2.502 907 875 095 892 822 283 902 873 218... , 
  http://en.wikipedia.org/wiki/Feigenbaum_constant
'''

import decimal
from decimal import Decimal as D

NUM_ITERS = 20
A_ITERS = 10

def calcAIter(n, a):
  # perform Briggs equations 7, 8, 9 (first 3)
  for _ in range(A_ITERS):
    b = D(0)
    bp = D(0)
    for _ in range(2**n):
      bp = 1 - 2*bp*b   # represents bPrime()
      b = a - b**2      # represents b(); equ 5
    a = a - (b / bp)
  return a, bp


d = D('3.2')    # will become the feigenbaum constant
a2 = D(0)  # the ordering is [ a2, a1, a, ...]
a1 = D(1)
bpOld = D(1)
# perform Briggs equations 10, 11, 12 (last three)
print(" i      delta          alpha")
for i in range(2, NUM_ITERS+1):
  a, bp = calcAIter(i, a1 + ((a1 - a2) / d))
  d = (a1 - a2)/(a - a1)  
  alpha = d*bpOld/bp
  bpOld = bp
  print(f"{i:2d}  {d:.10f}   {alpha:.10f}")
  # prepare for next iteration; move elems backwards
  a2 = a1
  a1 = a