# Degree-11 root art, after the 2025 image by Simone Conradi (profconradi.github.io):
# -i x^11 + x^10 + (t2 (-0.537-0.114i) - 0.571+0.414i) x^9 + x^8 + q(t1) x^7 + i x^6 - i x^5 + x^4 + i x^3 - x^2 + i x - i,
# q(t) = t^5 (-18.8-3.25i) + t^4 (-4.85+14i) + t^3 (5.81+12.6i) + t^2 (-8.74+4.29i) + t (5.81+0.687i) - 2.44+3.23i,   |t1| = |t2| = 1.
# Needs: pip install numpy matplotlib
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import hsv_to_rgb

N = 100_000          # number of pairs (t1, t2); try 10_000 first
SIZE = 1000
HALF = 3.0
PAPER = (0.957, 0.941, 0.902)

# Colour scheme: "paper" (translucent pastel ink on a light background), "script" (hue = angle of the first
# parameter, saturation 0.8, brightness 0.2 ... 1 from the second angle, on black) or "gradient" (your own colours below)
SCHEME = "paper"
GRADIENT = [(0.0, "#05030f"), (0.25, "#3a0f6b"), (0.5, "#c2185b"), (0.75, "#ff8a3d"), (1.0, "#fff3c4")]   # (position 0..1, colour)
GRADIENT_FOR = "density"   # "density" (dark = few roots) or "source" (the gradient is looped over the colour source; make the last colour equal the first)
COLOR_BY = "index"   # "index" (position of the root in the eigenvalue list), "t1" or "t2"

rng = np.random.default_rng(2)
a1, a2 = rng.random(N) * 2 * np.pi, rng.random(N) * 2 * np.pi
t1, t2 = np.exp(1j * a1), np.exp(1j * a2)

q = (t1**5 * (-18.8 - 3.25j) + t1**4 * (-4.85 + 14j) + t1**3 * (5.81 + 12.6j)
     + t1**2 * (-8.74 + 4.29j) + t1 * (5.81 + 0.687j) - 2.44 + 3.23j)
c9 = t2 * (-0.537 - 0.114j) - 0.571 + 0.414j
fixed = {11: -1j, 10: 1, 8: 1, 6: 1j, 5: -1j, 4: 1, 3: 1j, 2: -1, 1: 1j, 0: -1j}   # coefficients that do not move
coef = np.zeros((N, 12), dtype=complex)           # column k holds the coefficient of x^(11-k)
for deg, c in fixed.items():
    coef[:, 11 - deg] = c
coef[:, 11 - 9] = c9
coef[:, 11 - 7] = q
coef = coef / coef[:, :1]                         # divide by the leading coefficient -i: the polynomial becomes monic

C = np.zeros((N, 11, 11), dtype=complex)          # companion matrices
C[:, 1:, :-1] = np.eye(10)
C[:, 0, :] = -coef[:, 1:]
roots = np.linalg.eigvals(C)                      # shape (N, 11)

if COLOR_BY == "t1":
    hue = np.repeat(a1 / (2 * np.pi), 11)
elif COLOR_BY == "t2":
    hue = np.repeat(a2 / (2 * np.pi), 11)
else:
    hue = np.tile(np.arange(11) / 11, N)
x, y = roots.real.ravel(), roots.imag.ravel()

kw = dict(bins=SIZE, range=[[-HALF, HALF], [-HALF, HALF]])
cnt, _, _ = np.histogram2d(y, x, **kw)
cs, _, _ = np.histogram2d(y, x, weights=np.cos(2 * np.pi * hue), **kw)
sn, _, _ = np.histogram2d(y, x, weights=np.sin(2 * np.pi * hue), **kw)
un, _, _ = np.histogram2d(y, x, weights=np.repeat(a2 / (2 * np.pi), 11), **kw)   # mean of the second angle, for the script colours

dens = np.log1p(cnt / max(np.percentile(cnt[cnt > 0], 6), 1e-9))
v = np.clip(dens / np.percentile(dens[cnt > 0], 99.7), 0, 1)
alpha = np.clip(0.9 * v**0.55, 0, 1)[..., None]
h = (np.arctan2(sn, cs) / (2 * np.pi)) % 1
conc = np.clip(np.hypot(cs, sn) / np.maximum(cnt, 1) * 1.5 + 0.3, 0, 1)
ink = hsv_to_rgb(np.dstack([h, 0.78 * conc, np.full_like(h, 0.93)]))
if SCHEME == "paper":
    img = np.array(PAPER) * (1 - alpha) + ink * alpha
elif SCHEME == "script":
    m = np.maximum(cnt, 1)
    bright = 0.2 + 0.8 * np.clip(un / m, 0, 1)
    col = hsv_to_rgb(np.dstack([h, np.full_like(h, 0.8), bright]))
    a2_ = np.clip(v**0.6, 0, 1)[..., None]
    img = col * a2_                                  # on black
else:
    from matplotlib.colors import to_rgb
    pos = np.array([p for p, _ in GRADIENT]); cols = np.array([to_rgb(c) for _, c in GRADIENT])
    lut = lambda t: np.stack([np.interp(t, pos, cols[:, k]) for k in range(3)], axis=-1)
    if GRADIENT_FOR == "density":
        img = lut(v)
    else:
        img = lut(h) * alpha + np.array(PAPER) * (1 - alpha)

plt.figure(figsize=(8, 8))
plt.imshow(img[::-1], extent=[-HALF, HALF, -HALF, HALF])
plt.axis("off")
plt.savefig("poly11_art.png", dpi=150, bbox_inches="tight", pad_inches=0)
print("saved poly11_art.png")