Demonstrates how to compute a Julia fractal efficiently in pure Python + Numpy

vectorized-julia.py

import numpy
from matplotlib import pyplot as plot

N = 256
C = -0.4 + 0.6j
max_epoch_count = 1024
    
# Initialization
T = numpy.linspace(-2., 2., N)
X = numpy.tile(T, (N, 1))
Y = X.T
Z = X + 1j * Y
Z_count = numpy.zeros((N, N))
    
# Iterations
for epoch in range(max_epoch_count):
  mask = numpy.abs(Z ) < 2.
  Z[mask] = Z[mask] ** 2 + C
  Z_count[mask] += 1
        
# Plot
plot.imshow(Z_count, cmap = 'Spectral')
plot.show()