# Feigenbaum constant calculation
# FB - 201011151
# http://en.wikipedia.org/wiki/Feigenbaum_constant
# Source:
# http://keithbriggs.info/documents/how-to-calc.pdf
import math

def b(a, k):
    if k == 0:
        return 0
    else:
        return a - math.pow(b(a, k - 1), 2)

def b_prime(a, k):
    if k == 0:
        return 0
    else:
        return 1 - 2 * b_prime(a, k - 1) * b(a, k - 1)

# main
print 'Feigenbaum constant calculation:'
maxIt = 10
maxItJ = 10
ai_1 = 1
ai_2 = 0
di_1 = 3.2
print 'i', 'di'
for i in range(2, maxIt):
    ai = ai_1 + (ai_1 - ai_2) / di_1
    for j in range(maxItJ):
        ai = ai - b(ai, math.pow(2, i)) / b_prime(ai, math.pow(2, i))
    #
    di = (ai_1 - ai_2) / (ai - ai_1)
    print i, di
    # prepare for the next iteration
    di_1 = di
    ai_2 = ai_1
    ai_1 = ai