Solving Quintics - Introduction

In this section, we explain how to exploit the icosahedral symmetry discussed before to solve an arbitrary polynomial of order five (also known as a quintic polynomial).

We begin by returning to a result derived earlier: in the Section on icosahedral symmetry, we showed that the the icosahedral face centres can be grouped into five octohedra:

\begin{aligned} T(u,v) &= \prod_{r=1}^5 t_I(u\,\epsilon^{-r/2},v\,\epsilon^{r/2}) = \prod_{r=1}^5 \epsilon^{15-3r} \,t_I(u,v\,\epsilon^r) \\ &= u^{30} +522\left(u^{25}\,v^5-u^5\,v^{25}\right) - 10005\left(u^{20}\,v^{10} + u^{10}\,v^{20}\right) + v^{30} \,. \end{aligned}

(1)

where

t_I(u,v) := (u^2+v^2)(u^2-2z_n\,uv-v^2)(u^2-2z_m\,uv-v^2) \,.

(2)

# Check of the above:
from utils import e_simp, reduce_multiply
from sympy import symbols, expand, simplify, discriminant


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(u, v):
    return (
        u**30
        + 522 * (u**25 * v**5 - u**5 * v**25)
        - 10005 * (u**20 * v**10 + u**10 * v**20)
        + v**30
    )


T_eval = e_simp(
    expand(reduce_multiply([e ** (15 - 3 * k) * tI(u, v * e**k) for k in range(5)])), e
)
T_eval.equals(T(u, v))

True

Following the work on octahedral invariants, it’s usef

What we intend to do next is to demonstrate the existence of a map from points $z$ on the complex plane (away from the edge midpoints) to a particular set of roots of a Brioschi quintic, whose Brioschi parameter can also be computed from $z$ .

To do this, we begin by considering the polynomial $P(X)$ given by

P(X) = \prod_{r=1}^5 \left(X - \epsilon^{15-3r} \,t_I(u,v\,\epsilon^r)\right) \,.

(3)

Expanding, we find that

P(X) = X^5 - 10\,f(u,v)\,X^3+45\,f(u,v)^2\,X - T(u,v)

(4)

where $f(u,v) = uv\left(u^{10}+11\,u^5v^5-v^{10}\right)$ , as in the previous Section.

# Check of the above:
X = symbols("X")
result = e_simp(
    reduce_multiply([X - e ** (15 - 3 * k) * tI(u, v * e**k) for k in range(5)]), e
).collect(X)
result

Loading...

# Check the whole polynomial
def f(u, v):
    return u * v * (u**10 + 11 * u**5 * v**5 - v**10)


result.equals(X**5 - 10 * f(u, v) * X**3 + 45 * f(u, v) ** 2 * X - T(u, v))

True

Therefore, if we define the homogeneous coordinates

\tilde t_r(u,v) := -\frac{f^2(u,v)}{T(u,v)}\,\epsilon^{15-3r}\,t_I(u,v\,\epsilon^r) \,,

(5)

(which are degree-zero complex functions in $u$ and $v$ , defined everywhere except for at the edge midpoints) then a second polynomial can be readily constructed as follows:

P_B(X,z) = \prod_{r=1}^5 \left(X- \tilde t_r(z,1)\right) \,.

(6)

The polynomial $P_B(X)$ is indeed a Brioschi quintic with Brioschi parameter $B(z) = B(z,1)$ where

B(u,v) := -\frac{f^5(u,v)}{T^2(u,v)}

(7)

such that

P_B(X,z) = X^5 +10\,B(z)\, X^3 + 45\,B^2(z)\,X + B^2(z) \,.

(8)

For completeness, we recall the syzygy derived earlier which relates the three icosahedral polynomial invariants $f(u,v)$ , $T(u,v)$ , and $H(u,v)$ :

T^2(u,v) - H^3(u,v) = 12^3\,f^5(u,v)

(9)

where

H(u,v) := u^{20} - 228\left(u^{15}v^5-u^5v^{15}\right)+494\,u^{10}v^{10}+v^{20} \,.

(10)

def H(u, v):
    return u**20 - 228 * (u**15 * v**5 - u**5 * v**15) + 494 * u**10 * v**10 + v**20


(T(u, v) ** 2 - H(u, v) ** 3).equals(12**3 * f(u, v) ** 5)

True

We also define a second degree-zero ratio of polynomial invariants, namely

I(u,v) := -\frac{H^3(u,v)}{12^3f^5(u,v)} = 1 + \frac{1}{12^3B(u,v)} \,,

(11)

for future use, which conveniently maps the icosahedral vertices, edge midpoints, and face centres to infinity, one, and zero, respectively.

Tschirnhaus Transformation and Principal Quintics¶

In this section, we aim to demonstrate that solving a single quadratic equation allows us to trade any general quintic of the form

P(x) = x^5+ p_1 \,x^4 + p_2 \, x^3 + p_3 \, x^2+p_4 \, x + p_5 = 0

(12)

for a principal quintic with has the $x^4$ and $x^3$ coefficients equal to zero.

To this end, let’s define a transformation

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

(13)

and show how to solve for $(A,B)$ such that the polynomial equation $Q(y)=0$ with coefficients $q_i$ with $q_1 = q_2 = 0$ .

The polynomial $Q(y)$ after this transformation is given by the Resolvent, $\mathrm{Res}[P(x), y-x^2-Ax-B](y)$ :

\begin{aligned} Q(y) &= y^{5} + y^{4} \left(A p_{1} - 5 B - p_{1}^{2} + 2 p_{2}\right) \\ &+ y^{3} \left(A^{2} p_{2} - 4 A B p_{1} - A p_{1} p_{2} + 3 A p_{3} + 10 B^{2} + 4 B p_{1}^{2} - 8 B p_{2} - 2 p_{1} p_{3} + p_{2}^{2} + 2 p_{4}\right) + \ldots \end{aligned}

(14)

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

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 = symbols("y A B")

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

Loading...

Solving for the coefficients $A$ and $B$ , we find

\begin{aligned} A &= \frac{5 B + p_{1}^{2} - 2 p_{2}}{p_{1}} \\ B &=\frac{p_{1}^{2} p_{2} + 3 p_{1} p_{3} - 4 p_{2}^{2} \pm p_{1} \sqrt{\frac{8 p_{1}^{3} p_{3} - 3 p_{1}^{2} p_{2}^{2} + 16 p_{1}^{2} p_{4} - 38 p_{1} p_{2} p_{3} + 12 p_{2}^{3} - 40 p_{2} p_{4} + 45 p_{3}^{2}}5} }{2\left(2 p_{1}^{2} - 5 p_{2}\right)} \,. \end{aligned}

(15)

# Check of the above:
y4_coeff = expr2.coeff(y**4)
A_val = solve(y4_coeff, A)[0]
Eq(A, A_val)

Loading...

exprA2 = expr2.coeff(y**4).subs(p1, 0)
exprA2

Loading...

exprB2 = expr2.coeff(y**3).subs(p1, 0).subs(p3, 0).subs(B, 2 * p2 / 5)
exprB2

Loading...

A_vals = solve(exprB2, A)
A_vals[1]

Loading...

exprB = (expr2.coeff(y**3).subs(A, A_val) * p1**2).expand().collect(B)
exprB

Loading...

d = discriminant(exprB, B)
d

Loading...

q4 = symbols("q4")
d.subs(p1, 1).subs(p2, 1).subs(p3, 1).subs(p4, 1 - q4 / 120)

Loading...

B_vals = solve(expr2.coeff(y**3).subs(A, A_val), B)
B_vals[0]

Loading...

p1_part = 3
p2_part = 6
p3_part = 7
p4_part = -51
p5_part = -54

B_part = (
    B_vals[0]
    .subs(p1, p1_part)
    .subs(p2, p2_part)
    .subs(p3, p3_part)
    .subs(p4, p4_part)
    .subs(p5, p5_part)
)
B_part

Loading...

A_part = A_val.subs(B, B_part).subs(p1, p1_part).subs(p2, p2_part)
A_part

Loading...

expr2_part = (
    expr2.subs(p1, p1_part)
    .subs(p2, p2_part)
    .subs(p3, p3_part)
    .subs(p4, p4_part)
    .subs(p5, p5_part)
    .subs(A, A_part)
    .subs(B, B_part)
    .expand()
    .collect(y)
)
expr2_part

Loading...

expr2_part = (
    expr2.subs(p1, p1_part)
    .subs(p2, p2_part)
    .subs(p3, p3_part)
    .subs(p4, p4_part)
    .subs(p5, p5_part)
    .subs(A, 4)
    .subs(B, -8)
    .expand()
    .collect(y)
)
expr2_part

Loading...

alpha_part = expr2_part.coeff(y, 2) / 5
alpha_part

Loading...

beta_part = expr2_part.coeff(y, 1) / 5
beta_part

Loading...

gamma_part = expr2_part.coeff(y, 0)
gamma_part

Loading...

B_vals[1].subs(p1, 0).subs(p3, 0)  # .subs(p4, 1 - q4*15/2)

Loading...

A_val.subs(B, 2 * p2 / 5)

Loading...

B_vals[1].subs(p1, 1).subs(p2, 1).subs(p3, 1).subs(p4, 1 - q4 * 15 / 2)

Loading...

expr3 = (
    expr2.subs(p1, 0)
    .subs(p3, 0)
    .subs(B, 2 * p2 / 5)
    # .subs(A**5, (3 * p2 / 5 - 2 * p4 / p2) ** 2 * A)
    # .subs(A**3, (3 * p2 / 5 - 2 * p4 / p2) * A)
    # .subs(A**2, 3 * p2 / 5 - 2 * p4 / p2)
    .expand()
    .collect(y)
)
expr3

Loading...

from sympy import Rational


alpha = expr3.coeff(y**2) * Rational(1, 5)
beta = expr3.coeff(y) * Rational(1, 5)
gamma = expr3.coeff(y, 0)

(alpha**4 - beta**3 + alpha * beta * gamma).expand()

Loading...

Therefore, at the expense of a single square root, we can trade any general quintic for a principal quintic of the form

Q(y) = 0 \quad{\rm where}\quad Q(y) = y^5 + 5\alpha\, y^2 + 5\beta\, y + \gamma

(16)

for some $\alpha$ , $\beta$ , $\gamma$ computable in terms of the original $p_i$ . (The factors of 5 are included for later convenience.)

From Brioschi and Principal Quintics to the Icosahedral Equation¶

Returning to the Brioschi quintic, we wish to show that every Brioschi quintic can be associated with a corresponding principal quintic by transforming the roots as follows:

y_r = W_r \left(\sigma + \tau\,t_r\right)

(17)

where $W_r(u,v)$ is the polynomial invariant for the faces of the $r$ -th inscribed octahedron:

\begin{aligned} W_r(u,v) &= \epsilon^{20-4r}\,u^8 - \epsilon^{15-3r}\,u^7v+ 7\,\epsilon^{10-2r}\,u^6v^2+7\,\epsilon^{5-r}\,u^5v^3 \\ &\quad -7\,\epsilon^{r}\,u^3v^5+7\,\epsilon^{2r}\,u^2v^6 + \epsilon^{3r}\,uv^7 +\epsilon^{4r}\,v^8\,. \end{aligned}

(18)

(Recall that the roots of the original, Brioschi quintic were given by $\tilde t_r := -f^2\,t_r/T$ and so

y_r = W_r\left(\sigma - \frac{T\,\tau}{f^2}\,\tilde t_r\right)

(19)

in terms of the $\tilde t_r$ ’s.)

from utils import hessian, z_simp


e, r = symbols("e r")
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))

Summing over $r$ , we find that both $\sum_r Y_r$ and $\sum_r Y_r^2$ vanish:

# Check of the above:
def sum_over_r(expr, r):
    return e_simp(reduce_multiply([expr.subs(r, x) for x in range(5)]), e)


sigma, tau = symbols("sigma tau")
Y_r = sigma * W_r + tau * t_r * W_r

# sum_over_r(Y_r, r).equals(0) and sum_over_r(Y_r**2, r).equals(0)

Expanding the polynomial $P_y:= \prod_r(y-y_r)$ , we find a principal quintic with

\begin{aligned} \alpha &= -8\,f^2\,\sigma^3 - T\, \sigma^2\tau -72\,f^3\, \sigma\tau^2-fT\,\tau^3 \\ \beta &= H \left[f\,\sigma^4 - 18\,f^2\,\sigma^2\tau^2- T\,\sigma\tau^3-27\,f^3\,\tau^4\right]\\ \gamma &= - H^2 \left[ \sigma^5-10\,f\,\sigma^3\tau^2+45\,f^2\,\sigma\tau^4\,+T\,\tau^5\right] \end{aligned}

(20)

# Check of the above:
y = symbols("y")
P_y = e_simp(
    reduce_multiply([y - e_simp(Y_r.subs(r, x), e) for x in range(5)]).expand(),
    e,
)
# P_y.coeff(y**5), P_y.coeff(y**4), e_simp(P_y.coeff(y**3), e)

P_y_y2 = e_simp(P_y.coeff(y**2), e) / 5
P_y_y2.collect(sigma)

Loading...

P_y_y2.equals(
    -8 * f(u, v) ** 2 * sigma**3
    - T(u, v) * sigma**2 * tau
    - 72 * f(u, v) ** 3 * sigma * tau**2
    - f(u, v) * T(u, v) * tau**3
)

True

P_y_y1 = e_simp(P_y.coeff(y), e) / 5
P_y_y1.collect(sigma)

Loading...

P_y_y1.equals(
    H(u, v) * f(u, v) * sigma**4
    - 18 * H(u, v) * f(u, v) ** 2 * sigma**2 * tau**2
    - H(u, v) * T(u, v) * sigma * tau**3
    - 27 * f(u, v) ** 3 * H(u, v) * tau**4
)

True

P_y_y0 = e_simp(P_y - y**5 - 5 * y**2 * P_y_y2 - 5 * y * P_y_y1, e)
P_y_y0.collect(sigma)

Loading...

P_y_y0.equals(
    H(u, v) ** 2
    * (
        -(sigma**5)
        + 10 * sigma**3 * tau**2 * f(u, v)
        - 45 * sigma * tau**4 * f(u, v) ** 2
        - tau**5 * T(u, v)
    )
)

True

Since we will want the coefficients to have weight zero in $u$ , $v$ (so that they depend only on the ratio $z = u/v$ ), we will find it convenient to instead use

y_r = \tilde W_r \left(\lambda - \mu\,\tilde t_r\right)\qquad{\rm where}\quad \tilde W_r := \frac{W_r\,f}H \,,\quad \tilde t_r = -\frac{t_r\,f^2}T

(21)

since $\tilde W_r$ , $x_r$ are explicitly degree-zero in $u$ , $v$ . From here, we wish to show that, for every choice of $(\alpha,\beta,\gamma)$ defining a principal quintic, there’s a corresponding choice of $(\lambda,\mu,z)$ reproducing the same polynomial. As it turns out, solving for $\lambda$ , $\mu$ can be done using radicals (namely, a single square-root), whereas selecting the appropriate $z$ will require solving the Icosahedral Equation (to be defined shortly).

Switching to homogenous variables $(\lambda,\mu)$ we find

\begin{aligned} 12^3\,I(z) \, \alpha &= 8\,\lambda^3 + \lambda^2\mu - 72\,B(z)\,\lambda\mu^2 - B(z)\,\mu^3 \\ 12^3\,I(z) \, \beta &= - \lambda^4 -18\,B(z)\,\lambda^2\mu^2 - B(z)\,\lambda\mu^3 + 27\,B^2(z)\,\mu^4 \\ 12^3\,I(z) \, \gamma &= \lambda^5 + 10\,B(z)\,\lambda^3\mu^2 + 45\,B^2(z)\,\lambda\mu^4 + B^2(z)\,\mu^5 \end{aligned}

(22)

where, as before, $B(z) := -f^5(z,1)/T^2(z,1)$ .

Combining these equations, we notice that they satisfy

$$

\lambda,\beta + \gamma = -B(z),\mu^2,\alpha ,.

No input for math node

# Check of the above:
K, l, mu, B, alpha, beta, gamma = symbols("K lambda mu B alpha beta gamma")
# where K = 12^3*I(z)
alpha_expr = 1 / K * (8 * l**3 + l**2 * mu - 72 * B * l * mu**2 - B * mu**3)
beta_expr = (
    1 / K * (-(l**4) - 18 * B * l**2 * mu**2 - B * l * mu**3 + 27 * B**2 * mu**4)
)
gamma_expr = (
    1 / K * (l**5 + 10 * B * l**3 * mu**2 + 45 * B**2 * l * mu**4 + B**2 * mu**5)
)
(l * beta + gamma).subs(beta, beta_expr).subs(gamma, gamma_expr).equals(
    (-B * mu**2 * alpha).subs(alpha, alpha_expr)
)

True

Similarly, we can the compute a second combination and find

\lambda\,\gamma+\mu^2B(z)\beta = \frac{\left(\lambda^2+3\,B(z)\,\mu^2\right)^3}{12^3\,I(z)}

(24)

# Check of the above:
(l * gamma + mu**2 * B * beta).subs(gamma, gamma_expr).subs(beta, beta_expr).factor()

Loading...

We can obtain the same combination using the $\alpha$ and $\beta$ equations as follows:

27\,\alpha^2 + \frac{(\lambda\,\alpha + 8\,\beta)^2}{\mu^2 B(z)} = \frac{\left(\lambda^2+3\,B(z)\,\mu^2\right)^3}{12^3\,I(z)}

(25)

# Check of the above:
(K * (27 * alpha**2 + (l * alpha + 8 * beta) ** 2 / B / mu**2)).subs(
    alpha, alpha_expr
).subs(beta, beta_expr).subs(K, 1728 + 1 / B).factor()

Loading...

Equating these cubic terms, we find a quadratic expression in $\lambda$ , from which we can eliminate $\mu^2 B(z)$ using eq. mu^2 B and find

\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 \,.

(26)

# Check of the above:
alpha, beta, gamma = symbols("alpha beta gamma")
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...

At this point, we see that a second square root must be taken to determine $\lambda$ in terms of the values $(\alpha,\beta,\gamma)$ in the principal polynomial. Interestingly, the discriminant of the quadratic above,

\Delta := \alpha^2\left(\gamma^4 + 256\,\beta^5-320\,\alpha\beta^3\gamma+90\,\alpha^2\beta\gamma^2-135\,\alpha^4\beta^2+108\,\alpha^5\gamma\right) \,,

(27)

is proportional to the discriminant $D$ of the principal quintic:

D := (-1)^{n(n-1)/2} \, \mathrm{Res}[Q(y),Q'(y)] = \frac{5^5}{\alpha^2}\,\Delta \,.

(28)

# Check of the above:
l_disc = discriminant(l_eqn, l).factor()
l_disc

Loading...

discriminant(y**5 + 5 * alpha * y**2 + 5 * beta * y + gamma, y).equals(
    5**5 / alpha**2 * l_disc
)

True

Denoting the solutions to the quadratic equation for $\lambda$ as $\lambda_\pm$ , we can now write

$$ \mu^2B(z) = -\frac1\alpha\left(\beta,\lambda_\pm+\gamma\right) ,,

No input for math node

which allows us to solve for $I(z)$ in terms of the coefficients of the principal quintic: $$

I(z) = I_0 := \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)} ,.

$$

(Note that $I(z)$ has scaling dimension of zero in $y$ , since $[\alpha]=[y]^3$ , $[\beta]=[y]^4$ , $[\gamma] = [y]^5$ , and $[\lambda_\pm]=[y]$ .)

# Check of the above:
mu2B_eqn = l * gamma + mu**2 * B * beta - (l**2 + 3 * B * mu**2) ** 3 / K
K_soln = solve(mu2B_eqn, K)[0].subs(mu**2 * B, -1 / alpha * (beta * l + gamma)).factor()
K_soln

Loading...

Once $I(z)$ is known, the syzygy ensures that $B(z) = 1/12^3(I(z)-1)$ (away from the edge midpoints, where $I(z) = 1$ ) which allows us to solve for $\mu$ as well (at the expense of one last square root in Eq. (29)):

\mu = \mu_\pm := \pm\sqrt{\frac{-\beta\lambda_\pm-\gamma}{\alpha B(z_0)}} \,,

(30)

once we have some $z_0$ satisfying $I(z_0) = I_0$ . All in all, the roots $y_r$ of the principal quintic polynomial are given by

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) \,.

(31)

Inverting the map $I(z_0) = I_0$ to determine $z_0$ (given some $I_0$ computed in terms of $\alpha$ , $\beta$ , and $\gamma$ ) is the topic to which we now turn our attention.

Solving the Icosahedral Equation¶

The last step is to show how the equation $I(z=z_0) = I_0$ can be inverted to determine $z_0 := \phi(I=I_0)$ where $\phi := I^{-1}$ . As we shall see, this can only be achieved up to a Mobiüs transformation, but that is sufficient, since once a single value for $z$ is determined, then all five roots are identified using

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) \,.

(32)

J, K, L = symbols("J K L")
l = symbols("lambda")
alpha, beta, gamma = symbols("alpha beta gamma")
eqn = J * l**2 + K * l + L
lp, lm = solve(eqn, l)

(lp * lm).expand()

Loading...

(lp + lm).expand()

Loading...

l1 = -K * Rational(1, 2) / J
l1

Loading...

Jval = alpha**4 - beta**3 + alpha * beta * gamma
Kval = alpha * gamma**2 - 2 * beta**2 * gamma - 11 * alpha**3 * beta
Lval = 64 * alpha**2 * beta**2 - beta * gamma**2 - 27 * alpha**3 * gamma

denom = (
    12**3
    * alpha**2
    * (2 * J) ** 7
    * (
        (alpha * gamma - beta**2) ** 2 * L / J
        - 2 * (alpha * gamma - beta**2) * l1 * beta * gamma
        + beta**2 * gamma**2
    )
)
denom = denom.expand().factor()
denom

Loading...

denom2 = denom.subs(J, Jval).subs(K, Kval).subs(L, Lval).factor()
denom2

Loading...

s2 = (K**2 - 4 * J * L).subs(J, Jval).subs(K, Kval).subs(L, Lval).factor()
s2

Loading...

s = symbols("s")
l2 = Rational(1, 2) / J
lpm = l1 + s * l2
num = (
    (alpha * lpm**2 - 3 * beta * lpm - 3 * gamma) ** 3
    * (
        (alpha * gamma - beta**2) * l1
        - beta * gamma
        - s * l2 * (alpha * gamma - beta**2)
    )
    * (2 * J) ** 7
)
num.expand()

Loading...

num2 = num.expand().subs(J, Jval).subs(K, Kval).subs(L, Lval).expand().collect(s)
num2

Loading...

def s_simp(expr, s):
    poly = Poly(expr, s)
    coeffs = poly.all_coeffs()
    while len(coeffs) > 2:
        top_factor = coeffs[0]
        coeffs[2] += (K**2 - 4 * J * L) * top_factor
        coeffs.pop(0)
    return Poly.from_list(coeffs, s).as_expr()


num3 = s_simp(num2, s).subs(J, Jval).subs(K, Kval).subs(L, Lval).expand().collect(s)
num3.factor()

Loading...

denom2.factor()

Loading...

num4 = (
    num3.factor() / (3 * alpha * gamma - 4 * beta**2) ** 3 / Jval / (-64) / alpha**3
).collect(s)
num4

Loading...

num4.coeff(s, 0).factor()

Loading...

num4.coeff(s, 1).factor()

Loading...

denom3 = (
    denom2.factor() / (3 * alpha * gamma - 4 * beta**2) ** 3 / Jval / (-64) / alpha**3
).factor()
denom3

Loading...

num3.coeff(s, 0).factor()

Loading...

num3.coeff(s, 1).factor()

Loading...

A, B = symbols("A B")
alpha_val = A * B**2 - 56 * B**3
beta_val = -5 * A * B**3 + 408 * B**4
gamma_val = 44 * A * B**4 - 4704 * B**5 - B**4


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


A_simp(
    num3.coeff(s, 0)
    .subs(alpha, alpha_val)
    .subs(beta, beta_val)
    .subs(gamma, gamma_val)
    .expand(),
    A,
    B,
).factor()

Loading...

from sympy import sqrt

num_eval_s0 = _.subs(A, sqrt(15)).subs(B, -5).evalf()
num_eval_s0

Loading...

A_simp(
    num3.coeff(s, 1)
    .subs(alpha, alpha_val)
    .subs(beta, beta_val)
    .subs(gamma, gamma_val)
    .expand(),
    A,
    B,
).factor()

Loading...

num_eval_s1 = _.subs(A, sqrt(15)).subs(B, -5).evalf()
num_eval_s1

Loading...

s_eval = (
    s2.subs(alpha, alpha_val)
    .subs(beta, beta_val)
    .subs(gamma, gamma_val)
    .subs(A, sqrt(15))
    .subs(B, -5)
    .evalf()
)
s_eval

Loading...

num_eval_s0 + sqrt(s_eval) * num_eval_s1

Loading...

A_simp(
    s2.subs(alpha, alpha_val).subs(beta, beta_val).subs(gamma, gamma_val).expand(),
    A,
    B,
).factor()

Loading...

A_simp(
    num3.coeff(s, 1)
    .subs(alpha, alpha_val)
    .subs(beta, beta_val)
    .subs(gamma, gamma_val)
    .expand(),
    A,
    B,
).coeff(A, 1).factor()

Loading...

A_simp(
    (Jval**5)
    .subs(alpha, alpha_val)
    .subs(beta, beta_val)
    .subs(gamma, gamma_val)
    .expand(),
    A,
    B,
).factor()

Loading...

566155008 / 7776

72808.0

2381104350228864331230155414357818864626761245655040 / 532283307395841360870185084583173907477283001597952

4.473377836848294

denom4 = A_simp(
    denom3.subs(alpha, alpha_val).subs(beta, beta_val).subs(gamma, gamma_val).expand(),
    A,
    B,
)
denom4.factor()

Loading...

_.subs(A, sqrt(15)).subs(B, -5).evalf()

Loading...

3136 / 3

1045.3333333333333

d0 = denom4.coeff(A, 0)
d1 = denom4.coeff(A, 1)
(d0**2 - (-3 * B) * d1**2).factor()

Loading...