import matplotlib as mpl
from matplotlib import pyplot as plt
import numpy as np

mpl.rcParams.update({
    'text.usetex': True,
    'text.latex.preamble': r'\usepackage{amsfonts}'
})

N_ts = 3_000_000
ts = np.exp(np.random.uniform(0, 2*np.pi, 2 * N_ts) * 1j)
ts = ts.reshape((N_ts, 2))
K = ts[:,1]**4 - 1.0j * ts[:,1]**2 - 1.0
J = ts[:,0]**4 + ts[:,0]**2 - 1.0j * ts[:,0] - 1.0
z1 = (-K + np.sqrt(K**2 - 32*J)) / 16
z2 = (-K - np.sqrt(K**2 - 32*J)) / 16
x1 = np.sqrt(z1)
x2 = np.sqrt(z2)
x3 = -np.sqrt(z1)
x4 = -np.sqrt(z2)
pts = np.column_stack([x1, x2, x3, x4]).flatten()

fig,ax = plt.subplots(figsize=(30,30))
fig.set_facecolor("#f4f0e7")
fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None)
for spine in ['top', 'right','left','bottom']:
    ax.spines[spine].set_visible(False)
ax.scatter(
    x=pts.real,
    y=pts.imag,
    c="#262626",
    s=0.035,
    linewidths=1e-6)
ax.set_axis_off()
ax.set_aspect('equal')
ax.set_title('$8.0x^4+(1.0t_2^4-1.0it_2^2-1.0)x^2+1.0t_1^4+1.0t_1^2-1.0it_1-1.0$ \n' +  r'$t_1,t_2 \in \mathbb{C}, \ \ \ |t_1|=|t_2|=1$', fontsize=50)
plt.show()