Issue specifying datatype for the Mandelbrot set

Question

I have optimized a bit on calculating the Mandelbrot set, & I now wish to be able to specify whether my arrays should be float64 or float32 instead of the easier implementation with type complex128 or complex64. I use the fact that for a complex number (a+jb)^2 = a^2-b^2 + (2ab)j, but this seems to give me a slightly different wrong mandelbrot set. The code is seen below:

from timeit import default_timer as timer
import numpy as np
from numexpr import evaluate
import matplotlib.pyplot as plt

#%% Inputs
N = 5000
I = 20
T = 2 #Thresholdenter code here

#%% Functions
def mandel_brot_vector(I,C,T,datatype):
    Cre = np.array(C.real,dtype=datatype)
    Cim = np.array(C.imag,dtype=datatype)

    M = np.zeros(Cre.shape,dtype=datatype)
    zreal=0
    zimag=0
    for i in range(I):
        M[zreal*zreal+zimag*zimag<T**2] = i/I
        zreal = evaluate("zreal*zreal-zimag*zimag+Cre") #complex multiplication rule
        zimag = evaluate("2*zreal*zimag+Cim") #complex multiplication rule

    N = len(M[0])
    M = np.reshape(np.array(M),(N//2,N)).astype(datatype)
    M = np.concatenate((M,M[::-1]),axis=0)
    return M

def create_C(N,split):
    C_re = np.linspace(np.full((1,N),-2)[0],np.full((1,N),1)[0],N).T
    C_im = np.linspace(np.full((1,N),1.5*1j)[0],np.full((1,N),-1.5*1j)[0],N)
    C = C_re+C_im
    C = C[:N//2,:]

    if split != 0:
        C_split = np.array_split(C,split)
    else:
        C_split = C

    return np.array(C_split)

C = create_C(N, 0)

t0_32 = timer()
M32 = mandel_brot_vector(I,C,T,np.float32)
t_32 = timer() - t0_32

t0_64 = timer()
M64 = mandel_brot_vector(I,C,T,np.float64)
t_64 = timer() - t0_64



plt.matshow(M64,cmap="hot")

print(" "*10,f"N={N}")
print(f"{'Float 32':<20}{t_32:<40}",
      f"\n{'Float 64':<20}{t_64:<40}"
      )

Currently the image I get: wrong mandelbrot. For reference, the following function will produce the correct mandelbrot set but with complex128:

def mandel_brot(I,C,T):
    M = np.zeros(C.shape)
    z=0
    for i in range(I):
        M[np.abs(z)<T] = i/I
        z = evaluate("z*z+C")

    N = len(M[0])
    M = np.reshape(np.array(M),(N//2,N)).astype(datatype)
    M = np.concatenate((M,M[::-1]),axis=0)
    return M

Hope someone can help solve this issue, thanks in advance. Btw do not bother with the split of the C array, it is set up to run with multiprocessing which I am not using in the code attached.