Quintic Example - Introduction

In this section, we’ll work through the procedure for computing the roots of a quintic in the case of a specific example with known roots, three real and two complex, of the polynomial given by

\begin{aligned} P(x) &:= (x+3)(x+1)(x-2)(x^2+x+9) \\ &= x^5 + 3\,x^4+6\,x^3+7\,x^2-51\,x-54 \,. \end{aligned}

(1)

Since the roots are known, we can run the polynomial coefficients through our algorithm and verify that the appropriate roots are found.

The first step is to perform a Tschirnhaus transformation, in order to get the coefficients of the corresponding principal quintic. Taking

y = x^2+A\,x+B

(2)

with $A = -6$ , $B=-3$ , we find

Q(y) = y^5 -11065\,y^2-88260\,y+529056

(3)

and so $\alpha = -11065$ , $\beta = -17652$ , and $\gamma = 529056$ .

# Check of the above:
from sympy import sqrt, symbols, resultant, simplify


p1, p2, p3, p4, p5 = symbols("p_1 p_2 p_3 p_4 p_5")
x, y = symbols("x y")
P = x**5 + p1 * x**4 + p2 * x**3 + p3 * x**2 + p4 * x + p5
p_vals = {p1: 3, p2: 6, p3: 7, p4: -51, p5: -54}
B = (
    p1**2 * p2 / 2
    + 3 * p1 * p3 / 2
    + p1
    * sqrt(
        40 * p1**3 * p3
        - 15 * p1**2 * p2**2
        + 80 * p1**2 * p4
        - 190 * p1 * p2 * p3
        + 60 * p2**3
        - 200 * p2 * p4
        + 225 * p3**2
    )
    / 10
    - 2 * p2**2
) / (2 * p1**2 - 5 * p2)
B_val = B.subs(p_vals)
B_val

Loading...

A = (5 * B + p1**2 - 2 * p2) / p1
A_val = A.subs(p_vals)
A_val

Loading...

Q = resultant(P, y - x**2 - A_val * x - B_val, x).subs(p_vals)
Q

Loading...

alpha, beta, gamma = symbols("alpha beta gamma")
q_vals = {alpha: Q.coeff(y, 2) / 5, beta: Q.coeff(y, 1) / 5, gamma: Q.coeff(y, 0)}
q_vals

{alpha: -2213, beta: -17652, gamma: 529056}

Now that the coefficients of the principal quintic have been identified, we can proceed by solving for the values appearing in the expression

y = y_r := \frac{W_r(z_0,1)f(z_0,1)}{H(z_0,1)}\left(\lambda_\pm + \mu_\pm\,\frac{t_r(z_0,1)f^2(z_0,1)}{T(z_0,1)}\right) \,.

(4)

First of all, $\lambda_\pm$ is found by solving the quadratic equation

\left(\alpha^4-\beta^3+\alpha\beta\gamma\right)\lambda^2 + \left(\alpha\gamma^2-2\,\beta^2\gamma-11\,\alpha^3\beta\right)\lambda + 64\,\alpha^2\beta^2-27\,\alpha^3\gamma-\beta\gamma^2 =0 \,.

(5)

Solving for $\lambda$ , given the values of $(\alpha,\beta,\gamma)$ above, we find

\lambda = \lambda_\pm = \frac{5551876412706}{182368882091} \pm \frac{4470309903150 \sqrt{7} i}{182368882091} \,.

(6)

# Check of the above:
l, mu = symbols("lambda mu")
l_eqn = (
    simplify(
        -(beta * l + gamma)
        * alpha
        * (
            27 * alpha**2
            + (l * alpha + 8 * beta) ** 2 / mu**2 / B
            - l * gamma
            - mu**2 * B * beta
        ).subs(mu**2, -1 / B * (l * beta + gamma) / alpha)
    )
    .expand()
    .collect(l)
)
l_eqn

Loading...

from sympy import solve


l_val = solve(l_eqn.subs(q_vals), l)[0]
l_val

Loading...

From here, the corresponding values of $I$ are

\begin{aligned} I &= I_\pm = \frac{(\alpha\lambda_\pm^2-3\,\beta\lambda_\pm-3\gamma)^3}{12^3\alpha^2(\alpha\gamma\lambda_\pm-\beta^2\lambda_\pm-\beta\gamma)} \\ &= \frac{15433393006320276128482185197792403047552759957120800000}{18338381126299237208601401250567585364768902047827121041} \pm \frac{2477631834917432833053465047627874055898615726894880000 \sqrt{7} i}{201722192389291609294615413756243439012457922526098331451} \,. \end{aligned}

(7)

# Check of the above:
I0 = (alpha * l**2 - 3 * beta * l - 3 * gamma) ** 3 / (
    1728 * alpha**2 * (alpha * gamma * l - beta**2 * l - beta * gamma)
)
I0

Loading...

from sympy import radsimp


I_val = simplify(radsimp(I0.subs(q_vals).subs(l, l_val)))
I_val

Loading...

We also need a value for $\mu$ , which we take to be

\begin{align*} \mu &= \mu_\pm := \pm\sqrt{\frac{-\beta\lambda_\pm-\gamma}{\alpha B(z_0)}}\\ &=\sqrt{\frac{1106519134731851526970487630634061104470044531690575633549858636228880896}{36787850718980739066465358072503539144355288353260839478659255944041} + \frac{1965995425815591592119531912388367125043216243943615692796437685016972800 \sqrt{7} i}{36787850718980739066465358072503539144355288353260839478659255944041}} \,. \end{align*}

(8)

(Recall that $B(z_0) = 1+1/(1728\,I(z_0))$ .)

# Check of the above:
from sympy import Rational


mu_val = (1728 * I0 * alpha**2 - 8 * alpha * l**3 - 72 * l * (beta * l + gamma)) / (
    l**2 * alpha + beta * l + gamma
)
mu_val

Loading...

mu_val = radsimp(mu_val.subs(q_vals).subs(l, l_val).expand())
mu_val

Loading...

Now, we do have to compute $z_0$ using the icosahedral equation below. (This is the only part of the computation that can’t be done using radicals.) In doing so, we recall that the inverse is given by a ratio of hypergeometric functions:

s(x) = - (1728\,x)^{1/5}\,\frac{{}_2F_1(b,1+b-c;1+b-a;1/x)}{{}_2F_1(a, 1+a-c;1+a-b;1/x)} \,.

(9)

from sympy import hyper


def g_ratio_inf(sg_0, sg_1, sg_inf, x):
    nu_0, nu_1, nu_inf = 3, 2, 5
    a = Rational(1, 2) * (
        1
        - sg_0 * Rational(1, nu_0)
        - sg_1 * Rational(1, nu_1)
        + sg_inf * Rational(1, nu_inf)
    )
    b = Rational(1, 2) * (
        1
        - sg_0 * Rational(1, nu_0)
        - sg_1 * Rational(1, nu_1)
        - sg_inf * Rational(1, nu_inf)
    )
    c = 1 - sg_0 * Rational(1, nu_0)
    num = lambda z: z ** (-b) * hyper([b, 1 + b - c], [1 + b - a], 1 / z)
    denom = lambda z: z ** (-a) * hyper([a, 1 + a - c], [1 + a - b], 1 / z)
    return -((1728) ** Rational(1, 5)) * num(x) / denom(x)


z0_val = g_ratio_inf(1, 1, 1, I_val).evalf()
z0_val

Loading...

With $z_0 := s(I_0)$ in tow, we can now compute the roots $y_r$ of the principal quintic:

y_r = \frac{W_r(z_0,1)f(z_0,1)}{H(z_0,1)} \left(\lambda_\pm + \mu_\pm\,\frac{t_r(z_0,1)f^2(z_0,1)}{T(z_0,1)}\right) \,.

(10)

from utils import hessian, z_simp


e, r = symbols("e r")
z_n, u, v = symbols("z_n u v")
t_0 = (
    (u**2 + v**2)
    * (u**2 - 2 * z_n * u * v - v**2)
    * (u**2 - 2 * (-1 - z_n) * u * v - v**2)
)
t_r = simplify(z_simp((e ** (15 - 3 * r) * t_0.subs(v, v * e**r)).expand(), z_n))
W_r = simplify(hessian(t_r, u, v).expand() / (-400))

from sympy import I, pi, exp
from icosahedron import f, T, H


def y_r(z_0, r_val):
    e_val = (exp(2 * I * pi / 5)).evalf()
    W_val = W_r.subs(u, z_0).subs(v, 1).subs(e, e_val).subs(r, r_val)
    t_val = t_r.subs(u, z_0).subs(v, 1).subs(e, e_val).subs(r, r_val)
    f_val = f.subs(u, z_0).subs(v, 1)
    H_val = H.subs(u, z_0).subs(v, 1)
    T_val = T.subs(u, z_0).subs(v, 1)
    return (W_val * f_val / H_val * (l_val + mu_val * t_val * f_val**2 / T_val)).evalf()

from sympy import conjugate


def ys(z_0, conj=False):
    if conj:
        return [conjugate(y_r(z_0, i)) for i in range(5)]
    return [y_r(z_0, i) for i in range(5)]


def xs(z_0, conj=False):
    yss = ys(z_0, conj=conj)
    return [solve(yss[i] - x**2 - A_val * x - B_val, x)[0] for i in range(5)]


def icosahedral_orbit(z_seed, tol=1e-10):
    """Return all 60 preimages from a regular, finite seed."""
    eps = exp(2 * I * pi / 5).evalf()
    a = eps - eps**4
    b = eps**3 - eps**2

    points = [z_seed]
    i = 0

    while i < len(points):
        z = points[i]
        candidates = (
            (eps * z).evalf(),
            ((a * z + b) / (b * z - a)).evalf(),
        )
        for candidate in candidates:
            already_seen = any(
                abs(candidate - old) <= tol * max(1, abs(candidate), abs(old))
                for old in points
            )
            if not already_seen:
                points.append(candidate)
                if len(points) > 60:
                    raise RuntimeError("Orbit did not close; try a larger tolerance.")
        i += 1

    return points

z_vals = icosahedral_orbit(z0_val)
len(z_vals)

60

for val in z_vals[:5]:
    print(ys(val))

[-8.49999999999997 + 20.7062792408486*I, -8.50000000000002 - 20.7062792408486*I, -11.0 + 5.5169898597468e-15*I, 24.0 - 2.28072769004051e-14*I, 4.0 + 5.97253095504374e-16*I]
[4.00000000000001 + 3.2816089255705e-15*I, -8.49999999999997 + 20.7062792408486*I, -8.50000000000002 - 20.7062792408486*I, -11.0 + 2.99115116356163e-15*I, 24.0 - 2.01198649754653e-14*I]
[-8.49999999999997 + 20.7062792408486*I, -11.0 + 1.62564325065576e-14*I, -8.5 - 20.7062792408486*I, 3.99999999999999 - 1.54896305200672e-14*I, 24.0 - 1.91632733986813e-16*I]
[24.0 - 2.95001653313753e-14*I, 4.00000000000001 - 2.21564845463706e-15*I, -8.49999999999998 + 20.7062792408487*I, -8.50000000000001 - 20.7062792408486*I, -11.0 + 2.10476167744611e-15*I]
[4.00000000000001 - 9.60273798412264e-15*I, -8.50000000000001 - 20.7062792408486*I, -8.49999999999997 + 20.7062792408486*I, 24.0 - 2.03325854417069e-14*I, -11.0 + 2.07939134660995e-15*I]

from sympy import Abs, N


def root_permutation(reference, roots, tol=1e-8):
    """Return labels p such that roots[i] ≈ reference[p[i] - 1]."""
    if len(reference) != len(roots):
        raise ValueError("Root lists must have equal lengths.")

    permutation = []
    for root in roots:
        matches = [
            j + 1 for j, ref in enumerate(reference) if Abs(N(root - ref, 25)) < tol
        ]
        if len(matches) != 1:
            raise ValueError("Root match is missing or ambiguous.")
        permutation.append(matches[0])

    if len(set(permutation)) != len(reference):
        raise ValueError("Root matching is not bijective.")

    return tuple(permutation)

all_permutations = []
for val in z_vals:
    all_permutations.append(root_permutation(ys(val), ys(z_vals[0])))
    all_permutations.append(root_permutation(ys(val, conj=True), ys(z_vals[0])))

all_permutations = []
for val in z_vals:
    all_permutations.append(root_permutation(xs(val), xs(z_vals[0])))
    all_permutations.append(root_permutation(xs(val, conj=True), xs(z_vals[0])))

from itertools import permutations

exp_perm = []
for x in permutations(range(1, 6)):
    exp_perm.append(x)

len(exp_perm)

120

exp_perm[:6]

[(1, 2, 3, 4, 5),
 (1, 2, 3, 5, 4),
 (1, 2, 4, 3, 5),
 (1, 2, 4, 5, 3),
 (1, 2, 5, 3, 4),
 (1, 2, 5, 4, 3)]

for perm in all_permutations:
    if perm in exp_perm:
        exp_perm.remove(perm)
    else:
        print("fail")

len(exp_perm)

0

y_vals = [y_r(z0_val, i) for i in range(5)]
y_vals

[-8.49999999999997 + 20.7062792408486*I,
 -8.50000000000002 - 20.7062792408486*I,
 -11.0 + 5.5169898597468e-15*I,
 24.0 - 2.28072769004051e-14*I,
 4.0 + 5.97253095504374e-16*I]

x_vals = [solve(y_vals[i] - x**2 - A_val * x - B_val, x) for i in range(5)]
x_vals

[[-0.499999999999996 - 2.9580398915498*I, 6.5 + 2.9580398915498*I],
 [-0.499999999999996 + 2.9580398915498*I, 6.5 - 2.9580398915498*I],
 [2.0 - 2.7584949298734e-15*I, 4.0 + 2.7584949298734e-15*I],
 [-3.0 + 1.90060640836709e-15*I, 9.0 - 1.90060640836709e-15*I],
 [-1.0 - 7.46566369380468e-17*I, 7.0 + 7.46566369380468e-17*I]]

from utils import reduce_multiply

reduce_multiply([(x - x_vals[i][1]) for i in range(5)]).expand()

Loading...

reduce_multiply([(x - x_vals[i][0]) for i in range(5)]).expand()

Loading...

reduce_multiply([(y - y_vals[i]) for i in range(2)]).expand()

Loading...

reduce_multiply([(x - x_vals[i][1]) for i in range(2)]).expand()

Loading...

((x - x_vals[0][0]) * (x - x_vals[-1][0])).expand()

Loading...

from utils import e_simp, hessian

u, v, e = symbols("u v e")
z_n = e**2 + e**3
z_m = e + e**4


def tI(u, v):
    return (
        (u**2 + v**2)
        * (u**2 - 2 * z_n * u * v - v**2)
        * (u**2 - 2 * z_m * u * v - v**2)
    )


def t_r(r, z):
    return e_simp(e ** (15 - 3 * r) * tI(z, v * e**r), e)


def W_r(r, z):
    return (
        simplify(e_simp(hessian(t_r(r, u), u, v).expand() / (-400), e))
        .expand()
        .subs(u, z)
        .subs(v, 1)
    )


t_r(2, u).subs(u, z0_val).evalf()

Loading...

W_r(2, u).subs(u, z0_val).evalf()

Loading...

from icosahedron import H, f, T

u, v = symbols("u v")
T.subs(u, z0_val).subs(v, 1).expand()

Loading...

# Check of the above:
from sympy import Add, Eq, Poly, discriminant, resultant, solve, symbols

x, p1, p2, p3, p4, p5 = symbols("x p1 p2 p3 p4 p5")
expr = x**5 + p1 * x**4 + p2 * x**3 + p3 * x**2 + p4 * x + p5
Eq(expr, 0)

Loading...

y, A, B, C = symbols("y A B C")

expr2 = Poly(resultant(expr, y - x**2 - A * x - C, x), y).terms()
expr2 = Add(*[b * y ** a[0] for a, b in expr2], evaluate=False)
expr2

Loading...

from sympy import Rational

expr3 = expr2.subs(p1, 0).subs(p3, 0).subs(C, 2 * p2 * Rational(1, 5))
expr3

Loading...

A2val = solve(expr3.coeff(y**3), A**2)[0]
A2val

Loading...

q4, q5 = symbols("q_4 q_5")
A2val.subs(p2, 10 * B).subs(p4, 45 * B**2 - 5 * B * q4).expand()

Loading...

expr4 = (
    expr3
    # .subs(A**5, A2val**2 * A)
    # .subs(A**4, A2val**2)
    # .subs(A**3, A2val * A)
    # .subs(A**2, A2val)
    .subs(A**2 * p2, (A2val * p2).expand())
    # .subs(A**2*p2**2, (A2val*p2**2).expand())
    # .subs(p4,Rational(45,100)*p2**2+q4)
    # .subs(p5,p2**2*Rational(1,100)+q5)
    # .subs(p2,10*B)
    .expand()
    .collect(y)
)
expr4

Loading...

alpha = Rational(1, 5) * expr4.coeff(y, 2)
alpha

Loading...

alpha.subs(A**2, A2val).subs(p2, 10 * B).subs(p4, 45 * B**2).subs(p5, B**2)

Loading...

beta = Rational(1, 5) * expr4.coeff(y, 1)
beta

Loading...

beta.subs(A**2, A2val).subs(A**4, A2val**2).subs(A**3, A2val * A).subs(p2, 10 * B).subs(
    p4, 45 * B**2
).subs(p5, B**2)

Loading...

gamma = expr4.coeff(y, 0)
gamma

Loading...

gamma.subs(A**5, A2val**2 * A).subs(A**2, A2val).subs(A**4, A2val**2).subs(
    A**3, A2val * A
).subs(p2, 10 * B).subs(p4, 45 * B**2).subs(p5, B**2)

Loading...

gamma2 = (
    gamma.subs(A**5, (-3 * B + q4) ** 2 * A)
    .subs(A**2, -3 * B + q4)
    .subs(p5, B**2 + q5)
    .subs(p4, 45 * B**2 - 5 * B * q4)
    .subs(p2, 10 * B)
    .expand()
)
gamma2

Loading...

gamma_B = gamma2.subs(q4, 0).subs(q5, 0)
gamma_B

Loading...

beta2 = (
    beta.subs(A**4, (-3 * B + q4) ** 2)
    .subs(A**3, (-3 * B + q4) * A)
    .subs(p4, 45 * B**2 - 5 * B * q4)
    .subs(p5, B**2 + q5)
    .subs(p2, 10 * B)
    .expand()
)
beta2

Loading...

beta_B = beta2.subs(q4, 0).subs(q5, 0)
beta_B

Loading...

alpha2 = (
    alpha.subs(A**2, -3 * B + q4)
    .subs(p4, 45 * B**2 - 5 * B * q4)
    .subs(p2, 10 * B)
    .subs(p5, B**2 + q5)
    .expand()
)
alpha2

Loading...

alpha_B = alpha2.subs(q4, 0).subs(q5, 0)
alpha_B

Loading...

l = symbols("lambda")
alpha, beta, gamma = symbols("alpha beta gamma")
l_eqn = (
    (
        (alpha**4 - beta**3 + alpha * beta * gamma) * l**2
        + (alpha * gamma**2 - 2 * beta**2 * gamma - 11 * alpha**3 * beta) * l
        + 64 * alpha**2 * beta**2
        - 27 * alpha**3 * gamma
        - beta * gamma**2
    )
    .expand()
    # .subs(A**2, A2val)
    # .subs(A**4, A2val**2)
    # .subs(A**3, A2val * A)
    # .expand()
    .collect(l)
)
l_eqn

Loading...

l_eqn_B = (
    l_eqn.subs(alpha, alpha_B)
    .subs(beta, beta_B)
    .subs(gamma, gamma_B)
    .expand()
    .collect(A)
)
l_eqn_B

Loading...

def A_simp(expr, A, B):
    poly = Poly(expr, A)
    coeffs = poly.all_coeffs()
    while len(coeffs) > 2:
        top_factor = coeffs[0]
        coeffs[2] += (-3 * B) * top_factor
        coeffs.pop(0)
    return Poly.from_list(coeffs, A).as_expr()


l_eqn_B2 = A_simp(l_eqn_B, A, B).expand().collect(l)
l_eqn_B2

Loading...

discriminant(l_eqn_B, l).factor()

Loading...

A_simp(discriminant(l_eqn_B2, l), A, B).factor()

Loading...

l_vals = solve(l_eqn_B2, l)
l_vals[0].factor()

Loading...

l_vals[1]

Loading...