Mandelbrot set Coloring Error Around Period-2 Bulb (not colormap related)

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()

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.