This is a minor change, with a 50% improvement to @gareth's answer. Changing `escaped = abs(z_live) > threshold` to `escaped = z_live.real**2 + z_live.imag**2>threshold**2` while less nice looking is about 50% faster because it saves a square-root of all the elements. To time it, I used t1 = time() image = julia_set(mapping=mapping, min_coordinate=-1.5 - 1j, max_coordinate=1.5 + 1j, iterations_count = 255, width=1920, height=1080) print(time() - t1) Before: 3.16s After 3.21s