# Root art: the roots of  P(x) = 16 x^10 + 500i (t2-1) x^7 + 500i (1 + i t1 + t1^2 - t1^3) x^6 + (25/16) x^2
# for many random pairs (t1, t2) on the unit circle, drawn as translucent pastel ink on paper.
# After the mathematical art of Simone Conradi (profconradi.github.io). Needs: pip install numpy matplotlib
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import hsv_to_rgb

N = 200_000          # number of pairs (t1, t2); try 20_000 first
SIZE = 1000          # pixels per side
C_CONST = 25 / 16    # the constant term of P; try 1.5626 or 18.345 (the picture's central detail grows with it)
HALF = 4.45          # the picture shows the square [-HALF, HALF]^2
PAPER = (0.957, 0.941, 0.902)   # background colour; try (0, 0, 0) for black

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

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

# P = x^2 * Q, with Q of degree 8; make Q monic and take the eigenvalues of its companion matrix
c5 = 500j * (t2 - 1) / 16                       # coefficient of x^5
c4 = 500j * (1 + 1j * t1 + t1**2 - t1**3) / 16   # coefficient of x^4
c0 = np.full(N, C_CONST / 16, dtype=complex)   # constant term
C = np.zeros((N, 8, 8), dtype=complex)
C[:, 1:, :-1] = np.eye(7)                        # ones below the diagonal
C[:, 0, 2] = -c5                                 # first row: minus the coefficients of x^7 ... x^0
C[:, 0, 3] = -c4
C[:, 0, 7] = -c0
roots = np.linalg.eigvals(C)                     # shape (N, 8); the other two roots of P are 0

# colour of each root: its hue is the angle of t1 (use a2 or the root index for other looks)
hue = np.repeat(a1 / (2 * np.pi), 8)
x, y = roots.real.ravel(), roots.imag.ravel()

# histograms: how many roots per pixel, and the average hue direction (cos and sin) per pixel
rng_ = [[-HALF, HALF], [-HALF, HALF]]
kw = dict(bins=SIZE, range=rng_)
n, _, _ = 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), 8), **kw)   # mean of the second angle, for the script colours

# tone map: logarithm of the density, scaled by a high percentile, then alpha-composite ink over the paper
dens = np.log1p(n / max(np.percentile(n[n > 0], 6), 1e-9))
v = np.clip(dens / np.percentile(dens[n > 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(n, 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(n, 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])   # [::-1]: imaginary axis pointing up
plt.axis("off")
plt.savefig("root_art.png", dpi=150, bbox_inches="tight", pad_inches=0)
print("saved root_art.png")