Question: I wrote some code to render the Mandelbrot set with continuous coloring, but the bulbs on the period-2 blub are not colored correctly. regions that should not escape are colored as though the escape after 100 iterations, which happens to be my iteration limit.  What I tried 1. Increasing resolution and max iteration count 2. Increasing precision (np.float32 -> np.float64). Unfortunately, I could not finish execution because of the time complexity. import numpy as np import matplotlib.pyplot as plt def mandelbrot(cmax, width, height, maxiter): real = np.linspace(-2, 0.5, width, dtype=np.float32) imag = np.linspace(0, cmax.imag, height // 2 + 1, dtype=np.float32) c_real, c_imag = np.meshgrid(real, imag) output = np.zeros((height // 2 + 1, width), dtype=np.uint16) z_real = np.zeros_like(c_real) z_imag = np.zeros_like(c_imag) mask = np.ones_like(c_real, dtype=bool) p = np.sqrt((c_real - 0.25)**2 + c_imag**2) cardioid = c_real < p - 2 * p**2 + 0.25 period2_bulb = (c_real + 1)**2 + c_imag**2 < 0.0625 mask[cardioid | period2_bulb] = False output[~mask] = maxiter - 1 epsilon = 1e-10 for i in range(maxiter): zr2 = z_real[mask] ** 2 zi2 = z_imag[mask] ** 2 z_imag_new = 2 * z_real[mask] * z_imag[mask] + c_imag[mask] z_real[mask] = zr2 - zi2 + c_real[mask] z_imag[mask] = z_imag_new diverged = zr2 + zi2 >= 4.0 abs_z = zr2 + zi2 + epsilon output[mask] = i + 1 - np.log(np.log(abs_z)) / np.log(2) mask[mask] = ~diverged output[output == maxiter - 1] = 0 full_output = np.vstack([np.flipud(output), np.flipud(output[-2::-1, :])]) return full_output width, height, maxiter = 12000, 9000, 100 cmax = complex(0, 1) mandelbrot_set = mandelbrot(cmax, width, height, maxiter) plt.figure(figsize=(16, 12)) plt.imshow(mandelbrot_set, extent=[-2, 0.5, -1, 1], cmap="inferno", interpolation="bilinear") plt.colorbar(label="Iterations to Divergence") plt.title("Mandelbrot Set") plt.xlabel("Re(c)") plt.ylabel("Im(c)") plt.show() 

Answer: The updating of the output and the mask inside the loop (iterations) is incorrect. We need to focus only on the points inside mask and escaped in current iteration. def mandelbrot(cmax, width, height, maxiter): real = np.linspace(-2, 0.5, width, dtype=np.float32) imag = np.linspace(0, cmax.imag, height // 2 + 1, dtype=np.float32) c_real, c_imag = np.meshgrid(real, imag) output = np.zeros((height // 2 + 1, width), dtype=np.uint16) z_real = np.zeros_like(c_real) z_imag = np.zeros_like(c_imag) mask = np.ones_like(c_real, dtype=bool) #p = np.sqrt((c_real - 0.25)**2 + c_imag**2) #cardioid = c_real < p - 2 * p**2 + 0.25 q = (c_real - 0.25)**2 + c_imag**2 cardioid = q*(q + c_real - 1/4) <= 1/4*c_imag**2 period2_bulb = (c_real + 1)**2 + c_imag**2 < 1/16 #0.0625 mask[cardioid | period2_bulb] = False # these points never escape epsilon = 1e-10 for i in range(maxiter): z_imag_new = 2 * z_real[mask] * z_imag[mask] + c_imag[mask] z_real[mask] = z_real[mask]**2 - z_imag[mask]**2 + c_real[mask] z_imag[mask] = z_imag_new diverged = z_real**2 + z_imag**2 >= 4.0 # escaped in this iteration abs_z = z_real[mask & diverged]**2 + z_imag[mask & diverged]**2 + epsilon output[mask & diverged] = i + 1 - np.log(np.log(abs_z)) / np.log(2) # only the points inside mask that escaped in the current iteration are need to be updated in the output mask[mask & diverged] = False # only the escaped points inside the mask are no longer needed to be considered #output[output == maxiter - 1] = 0 full_output = np.vstack([np.flipud(output), np.flipud(output[-2::-1, :])]) return full_output   Although this may not be a very efficient implementation and for more efficient implementations one could use numba jit , as suggested here.