# 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")
Diff to Previous Revision
--- revision 1 2013-09-22 22:25:07
+++ revision 2 2013-09-23 03:18:02
@@ -12,11 +12,9 @@
if b == 0: return a
return gcd(b, a % b)
-pix = [[0 for x in range(imgSize)] for y in range(imgSize)]
+pi2 = math.pi * 2.0
MAXX = imgSize + 1
MAXY = imgSize + 1
-pi2 = math.pi * 2.0
-nold = 0
qmax = imgSize
for q in range(4, qmax, 2):
print str(100 * q / qmax).zfill(2) + "%"
@@ -72,12 +70,8 @@
polynew = (2.0 * math.cos(sigma * q / 2.0) - e) * poly - 2.0 * polyold
if poly * polynew < 0.0: n += 1
if n > nold:
- pix[int(MAXX * p / q)][int(MAXY - ie)] = 1
- pix[int(MAXY - ie)][int(MAXX * p / q)] = 1
+ pixels[int(MAXY - ie), int(MAXX * p / q)] = (255, 255, 255)
+ pixels[int(MAXX * p / q), int(MAXY - ie)] = (255, 255, 255)
nold = n
-# paint
-for ky in range(imgSize):
- for kx in range(imgSize):
- pixels[kx, ky] = (pix[ky][kx] * 255, pix[ky][kx] * 255, pix[ky][kx] * 255)
image.save("HofstadterButterflyFractal.png", "PNG")
