# Hofstadter Butterfly Fractal
# http://en.wikipedia.org/wiki/Hofstadter%27s_butterfly
# Wolfgang Kinzel/Georg Reents,"Physics by Computer" Springer Press (1998)
# FB36 - 20130922
import math
from PIL import Image
imgSize = 200
image = Image.new("RGB", (imgSize, imgSize))
pixels = image.load() 

def gcd(a, b): # Greatest Common Divisor
    if b == 0: return a
    return gcd(b, a % b)

pi2 = math.pi * 2.0
MAXX = imgSize + 1
MAXY = imgSize + 1
qmax = imgSize
for q in range(4, qmax, 2):
    print str(100 * q / qmax).zfill(2) + "%"
    for p in range(1, q, 2):
        if gcd(p, q) <= 1:
            sigma = pi2 * p / q
            nold = 0
            ie = 0
            for ie in range(0, MAXY + 2):
                e = 8.0 * ie / MAXY - 4.0 - 4.0 / MAXY
                n = 0
                polyold = 1.0
                poly = 2.0 * math.cos(sigma) - e
                if polyold * poly < 0.0: n += 1

                for m in range(2, q / 2):
                    polynew = (2.0 * math.cos(sigma * m) - e) * poly - polyold
                    if poly * polynew < 0.0: n += 1
                    polyold = poly
                    poly = polynew

                polyold = 1.0
                poly = 2.0 - e
                if polyold * poly < 0.0: n += 1
                polynew = (2.0 * math.cos(sigma) - e) * poly - 2.0 * polyold
                if poly * polynew < 0.0: n += 1
                polyold = poly
                poly = polynew

                for m in range(2, q / 2):
                    polynew = (2.0 * math.cos(sigma * m) - e) * poly - polyold
                    if poly * polynew < 0.0: n += 1
                    polyold = poly
                    poly = polynew

                polynew = (2.0 * math.cos(sigma * q / 2.0) - e) * poly - 2.0 * polyold
                if poly * polynew < 0.0: n += 1

                polyold = 1.0
                poly = 2.0 - e
                if polyold * poly < 0.0: n += 1
                polynew = (2.0 * math.cos(sigma) - e) * poly - 2.0 * polyold
                if poly * polynew < 0.0: n += 1
                polyold = poly
                poly = polynew

                for m in range(2, q / 2):
                    polynew = (2.0 * math.cos(sigma * m) - e) * poly - polyold
                    if poly * polynew < 0.0: n += 1
                    polyold = poly
                    poly = polynew

                polynew = (2.0 * math.cos(sigma * q / 2.0) - e) * poly - 2.0 * polyold
                if poly * polynew < 0.0: n += 1
                if n > nold:
                    pixels[int(MAXY - ie), int(MAXX * p / q)] = (255, 255, 255)
                    pixels[int(MAXX * p / q), int(MAXY - ie)] = (255, 255, 255)
                nold = n

image.save("HofstadterButterflyFractal.png", "PNG")