Icosahedral Symmetry

Icosahedral Symmetry#

The icosahedral group is defined as the group of rotations that map the icosahedron into itself. Placing two antipodal vertices along the \(\zeta\) axis at \((0,0,\pm 1)\), and one next-nearest vertex at \({\bf\hat n}_n = (-\sin\,2\alpha,0,\cos\,2\alpha)\), we find that there is a two-fold symmetry about the axis \({\bf \hat n}_2:=(-\sin\alpha,0,\cos\alpha)\), in addition to a five-fold symmetry about \({\bf \hat n}_5:=(0,0,1)\):

\[\begin{split}\begin{align*} IV:&\quad {\mathcal O}({\bf\hat n}_2, \pi)\\ V:&\quad {\mathcal O}({\bf\hat n}_5, 2\pi/5) \end{align*}\end{split}\]

It can be shown that the above symmetry generators can be combined to reproduce the 60-fold symmetries of the icosahedral group. In particular, the three-fold symmetry of one of the icosahedron faces nearest to the top vertex can be written as

\[{\mathcal O}_C({\bf\hat n}_3, 2\pi/3) = {\mathcal O}_B({\bf\hat n}_5, -2\pi/5)\,{\mathcal O}_A({\bf\hat n}_2, \pi) \,,\]

given an appropriate choice of \({\bf\hat n}_2\) (or \(\alpha\)) and \({\bf\hat n_3}\).

To determine \(\alpha\), note that the \(\{{\bf a}, d\}\) parametrizations of the above transformations are

\[\begin{split}\begin{align*} a_{1A} = -\sin\alpha\,,\quad a_{3A} = \cos\alpha \,,\quad a_{2A} = d_A =0 \\ a_{1B} = a_{2B} = 0 \,,\quad a_{3B} = -\sin(\pi/5) \,,\quad d_B = \cos(\pi/5) \,. \end{align*}\end{split}\]

Recalling the composition property of the \(d\) parameter, \(d_C = -a_{1A} a_{1B} - a_{2A} a_{2B} - a_{3A} a_{3B} + d_A d_B\), we find that

\[d_C = \cos\alpha\,\sin(\pi/5) = \cos(\pi/3) = \frac12\]

for the composed transformation. Therfore,

\[\cos\alpha = \frac1{2\,\sin(\pi/5)} \,.\]

Similarly,

\[\begin{split}\begin{align*} a_{1C} &= a_{1A} d_B + a_{1B} d_A - a_{2A} a_{3B} + a_{2B} a_{3A} = -\sin\alpha\,\cos(\pi/5)\\ a_{2C} &= a_{2A} d_B + a_{2B} d_A - a_{3A} a_{1B} + a_{3B} a_{1A} = \sin\alpha\,\sin(\pi/5) \\ a_{3C} &= a_{3A} d_B + a_{3B} d_A - a_{1A} a_{2B} + a_{1B} a_{2A} = \cos\alpha\,\cos(\pi/5) \end{align*}\end{split}\]

which together imply that

\[{\bf\hat n}_3 = \frac2{\sqrt{3}}(-\sin\alpha\,\cos(\pi/5), \sin\alpha\,\sin(\pi/5), \cos\alpha\,\cos(\pi/5)) \,.\]

(\({\bf\hat n}_3\) is guaranteed to be unit-normalized whenever \(\alpha\) satisfies the condition above.)

Moving to stereographically-projected coordinates, the vertex at \({\bf\hat n}_n\) is given by

\[z_n := -\frac{\sin 2\alpha}{(1-\cos 2\alpha)} = -\cot\alpha \,.\]

However, we will find that it is more convenient to instead write

\[z_n = \epsilon^2 + \epsilon^3 \,,\quad{\rm where}\quad \epsilon:=e^{2\pi i /5} \,,\]

since the other four vertices related to \(z_n\) under the five-fold symmetry transformation \(V\) are given by

\[z_k := \epsilon^k\,z_n \,,\quad k \in \{1, 2, 3, 4, 5\} \,.\]
from sympy import I, Matrix, N, Poly, Rational, conjugate, cos, diff, oo, pi, simplify, sin, sqrt, symbols, tan, Eq, solve, evalf, acos, csc  # fmt: skip

u, v = symbols("u v")
(2 * cos(pi / 5)).evalf()
\[\displaystyle 1.61803398874989\]
((1 + sqrt(5)) / 2).evalf()
\[\displaystyle 1.61803398874989\]

The equivalence of the two expressions above follows from the defining relation

\[\epsilon^5 = 1 \quad\implies\quad \epsilon + \epsilon^2 + \epsilon^3 + \epsilon^4 + \epsilon^5 = 0 \]

which together imply that

\[(\epsilon^4 - \epsilon)(\epsilon^2 - \epsilon^3) = \epsilon + \epsilon^4 - \epsilon^2 - \epsilon^3 = \sqrt{5} \]

(as can be checked by squaring the second equality above). Noting that \(\sin (\pi/5) = (\epsilon^2-\epsilon^3)/2i\), the previous condition on \(\alpha\) becomes

\[\cos\alpha = \frac i{\epsilon^2-\epsilon^3} =\frac i{\sqrt{5}}(\epsilon^4 - \epsilon) \quad\implies\quad \sin\alpha = \frac i{\sqrt{5}}(\epsilon^3-\epsilon^2)\]

and so

\[z_n = - \cot\alpha = \frac{(\epsilon - \epsilon^4)}{(\epsilon^3-\epsilon^2)} \times\frac{(\epsilon^2+\epsilon^3)}{(\epsilon^2+\epsilon^3)} = \epsilon^2+\epsilon^3 \,.\]

Across from \(z_n\) is an antipodal vertex \(z_m:= -1/z_n\) given by

\[z_m = -\frac1{(\epsilon^2+\epsilon^3)} = \frac{-1\times(\epsilon+\epsilon^4)}{(\epsilon^2+\epsilon^3)(\epsilon+\epsilon^4)} = \epsilon+\epsilon^4 \,.\]

In addition to the vertices \(z_0 = 0\), \(z_\infty = \infty\) and the five \(z_k\) (where \(z_5 = z_n\)), the remaining five are related to \(z_m\) under the five-fold symmetry transformation \(V\) as follows:

\[z_k := \epsilon^k\,z_m \,,\quad k \in \{6, 7, 8, 9, 10\} \]

(where \(z_{10} = z_m\), similarly to \(z_n\)).

Relating \(z_n\), \(z_m\) back to to their explicit values, we see that

\[z_n = - 2\,\cos\pi/5 := -\phi \,,\quad z_m = 2\,\cos 2\pi/5 = \phi^2- 2\]

(where the last equality follows from the cosine half-angle formula). Recalling that \(z_n\) and \(z_m\) are antipodal, we see that \(\phi\) must satisfy

\[z_n\,z_m = -1 \quad\leftrightarrow\quad (\phi+1)(\phi^2 - \phi -1) =0 \,.\]
# Check of the above
z_n, z_m, phi = symbols("z_n z_m phi")

num = z_n * z_m + 1
Eq(num, 0)
\[\displaystyle z_{m} z_{n} + 1 = 0\]
solve(num.subs(z_n, -phi).subs(z_m, phi**2 - 2), phi)

[-1, 1/2 - sqrt(5)/2, 1/2 + sqrt(5)/2]

[x.evalf() for x in _]

[-1.00000000000000, -0.618033988749895, 1.61803398874989]

Therefore, \(\phi\) is the golden ratio: \(\phi = (1+\sqrt{5})/2\) (since \(\phi\), by definition, is positive). Additionally, since the \(\epsilon\)’s satisfy \(\sum_{k=1}^5\epsilon^k = 0\), we find

\[z_n + z_m + 1 =0 \,,\quad z_n^2+z_n-1=0\,,\quad z_m^2+z_m-1=0\,,\]

which is consistent with the above.

Lastly, the homogeneous transformation corresponding to the symmetry transformation \(IV\) is given by

\[\begin{split}\begin{align*} u' = i\,\cos\alpha \, u - i\,\sin\alpha\,v &= \frac1{\sqrt{5}}\left(\epsilon - \epsilon^4\right) u +\frac1{\sqrt{5}}\left(\epsilon^3 - \epsilon^2\right) v \\ v' = -i\,\sin\alpha\,u - i\,\cos\alpha\, v &= \frac1{\sqrt{5}}\left(\epsilon^3 - \epsilon^2\right) u -\frac1{\sqrt{5}}\left(\epsilon - \epsilon^4\right) v \,, \end{align*}\end{split}\]

whereas transformation \(V\) corresponds to \(u' = \epsilon^{1/2}\, u\), \(v' = \epsilon^{-1/2}\, v\).

Icosahedral Polynomial Invariants#

Multiplying together the homogeneous polynomials \((u-z_k\,v)\) for \(k\in \{0,1,\ldots,10,\infty\}\), we find

\[f(u,v) = \prod_k\left(u - z_k\,v\right) = uv\left(u^{10}+11\,u^5v^5-v^{10}\right) \,.\]
# Check of the above:
from functools import reduce

# e = exp(2 * pi * I / 5)
e = symbols("e")
n = e**2 + e**3
m = e + e**4

iso_verts = []
for i in range(5):
    iso_verts.append(n * e**i)
    iso_verts.append(m * e**i)

f = u * v * reduce(lambda x, y: x * y, [u - x * v for x in iso_verts])
f
\[\displaystyle u v \left(u - v \left(e^{3} + e^{2}\right)\right) \left(u - v \left(e^{4} + e\right)\right) \left(- e v \left(e^{3} + e^{2}\right) + u\right) \left(- e v \left(e^{4} + e\right) + u\right) \left(- e^{2} v \left(e^{3} + e^{2}\right) + u\right) \left(- e^{2} v \left(e^{4} + e\right) + u\right) \left(- e^{3} v \left(e^{3} + e^{2}\right) + u\right) \left(- e^{3} v \left(e^{4} + e\right) + u\right) \left(- e^{4} v \left(e^{3} + e^{2}\right) + u\right) \left(- e^{4} v \left(e^{4} + e\right) + u\right)\]

Given the complexity of this expression, it makes sense to introduce a function which will efficiently expand out any polynomial and \(\epsilon\) and enforce the conditions

\[\epsilon^5 = 1 \,,\qquad \epsilon^4 = - \epsilon^3 - \epsilon^2 - \epsilon - 1\]

such that the maximum power of \(\epsilon\) appearing is 3. These considerations motivate the definition of the function e_simp() below:

def e_simp(expr, e, full=True, cleanup=True):
    result = Poly(expr, e)
    e_coeffs = reversed(result.all_coeffs())
    simp_coeffs = [0] * 5
    # enforce e**5 = 1
    for i, coeff in enumerate(e_coeffs):
        simp_coeffs[i % 5] += coeff
    # enforce e**4 = - e**3 - e**2 - e - 1
    if full:
        simp_coeffs = [simp_coeffs[i] - simp_coeffs[4] for i in range(4)]
    if cleanup:
        return simplify(Poly.from_list(reversed(simp_coeffs), e).as_expr())
    return Poly.from_list(reversed(simp_coeffs), e).as_expr()


f = e_simp(f, e)
f
\[\displaystyle u v \left(u^{10} + 11 u^{5} v^{5} - v^{10}\right)\]

As a second check, it’s worth pointing out that, for any real antipodal pair \((z_n, z_m)\) with \(z_m=-1/z_n\), the following result holds:

(1)#\[\begin{split}\begin{aligned} P(u,v,z_n,z_m) :=\prod_{k=1}^5\left(u - \epsilon^k\,z_n\,v\right)\left(u - \epsilon^k\,z_m\,v\right) &= \prod_{k=1}^5\left[u^2-\epsilon^k\left(z_n + z_m \right)\,uv-\epsilon^{2k}\,v^2\right] \\ &= u^{10} - \left(z_n^5 + z_m^5\right)u^5v^5 - v^{10} \,. \end{aligned}\end{split}\]
from utils import reduce_multiply


e_simp(
    reduce_multiply([(u - e**k * z_n * v) * (u - z_m * e**k * v) for k in range(5)]),
    e,
    cleanup=False,
)
\[\displaystyle u^{10} - u^{5} v^{5} z_{m}^{5} - u^{5} v^{5} z_{n}^{5} + v^{10} z_{m}^{5} z_{n}^{5}\]
# Second check:
z_n = symbols("z_n")


def antipodal_pair(m):
    return u**2 - e**m * (z_n - 1 / z_n) * u * v - e ** (2 * m) * v**2


f2 = reduce(lambda x, y: x * y, [antipodal_pair(m) for m in range(5)])
e_simp(f2, e)
\[\displaystyle u^{10} - u^{5} v^{5} z_{n}^{5} + \frac{u^{5} v^{5}}{z_{n}^{5}} - v^{10}\]

This, combined with the product \((u-z_0\,v)(u-z_\infty\,v) = uv\), demonstrates the expected formula for the polynomial invariant \(f(u,v)\) associated with the icosahedral vertices, since

\[\begin{split}\begin{aligned} z_n^5 &= z_n\left(1-z_n\right)\left(1-z_n\right) \\ &= z_n\left[1-2\,z_n+\left(1-z_n\right)\right] = 2\,z_n -3\,\left(1-z_n\right) \\ &= 5\,z_n - 3 \end{aligned}\end{split}\]

and so

\[\begin{split}\begin{aligned} \frac1{z_n^5}-z_n^5 &= \frac1{\left(5\,z_n - 3\right)}\times\left(\frac{5\,z_m-3}{5\,z_m-3}\right) - \left(5\,z_n - 3\right) \\ &=\frac{5\,z_m-3}{-25 -15\left(z_n+z_m\right)+9} - \left(5\,z_n - 3\right) \\ &=-5\left(z_n+z_m\right)+6 \\ &= 11 \,. \end{aligned} \end{split}\]

(The above follows from repeated use of the defining relations \(z_n^2 = 1-z_n\), \(z_n\,z_m = -1\), and \(z_n+z_m=-1\).)

# Check of the above:
e_simp(n**5, e).equals(5 * n - 3)

True

e_simp(-(m**5) - n**5, e)
\[\displaystyle 11\]

It should be straightforward to check that \(f(u,v)\) is indeed a polynomial invariant under transformations \(IV\), \(V\):

\[f(u',v') = f(u,v) \,.\]
# Check of the above, under IV:
orth1 = 1 / sqrt(5) * Matrix([[e - e**4, e**3 - e**2], [e**3 - e**2, e**4 - e]])
u, v = symbols("u v")
U, V = symbols("U V")
uv = Matrix([u, v])

uv_prime = orth1 * uv

f_prime = f.subs([(u, U), (v, V)]).subs([(U, uv_prime[0]), (V, uv_prime[1])])
e_simp(f_prime, e).equals(f)

True

# Check of the above, under V:
orth2 = Matrix([[sqrt(e), 0], [0, 1 / sqrt(e)]])
uv = Matrix([u, v])

uv_prime = orth2 * uv
f_prime = f.subs([(u, U), (v, V)]).subs([(U, uv_prime[0]), (V, uv_prime[1])])
e_simp(e**5 * f_prime, e).equals(f)

True

Computing the Hessian of \(f\), we find \({\rm Hess}[f](u,v) = -121\,H(u,v)\) where

\[H(u,v) := u^{20} + v^{20} - 228\left(u^{15}\,v^5-u^5\,v^{15}\right) + 494\,u^{10}\,v^{10} \,.\]
# Check of the above:
from utils import hessian

H = simplify(hessian(f, u, v) / (-121))
H
\[\displaystyle u^{20} - 228 u^{15} v^{5} + 494 u^{10} v^{10} + 228 u^{5} v^{15} + v^{20}\]

We shall see later that this Hessian is precisely the polynomial invariant for the face centres of the icosahedron.

The Rotated Octahedron#

Earlier, we demonstrated the that function \(t(u,v) = u^5\,v-u\,v^5\) was an absolute invariant under tetrahedral symmetry, and that its square was an absolute invariant under octahedral symmetry. We now wish to consider how this particular product of roots changes under a rotation of the octahedron about its centre.

To this end, consider the rotated points \(z'_k\), related to \(z_k \in \{0, \infty, \pm 1,\pm i\}\) by

\[z'_k(\alpha, \theta) = e^{i\,\theta}\left[\frac{\cos(\alpha/2)\,z_k+\sin(\alpha/2)}{-\sin(\alpha/2)\,z_k + \cos(\alpha/2)}\right]\,.\]

(The rotation above corresponds to a rotation about the \(\eta\) axis by an angle \(\alpha\), followed by a rotation about the \(\zeta\) axis by an angle \(\theta\).) Computing the product \(t'(u,v) := \prod_{k}(u-z_k'\,v)\) we find

\[\begin{split}\begin{align*} t'(u,v, \tan(2\alpha), e^{i\theta}) &= u^6 + \frac{4\,u^5\,v\,e^{i\theta}}{\tan(2\alpha)}-5\,u^4\,v^2\,e^{2\,i\theta} \\ &- 5\,u^2\,v^4\,e^{4\,i\theta} - \frac{4\,u\,v^5\,e^{5\,i\theta}}{\tan(2\alpha)} + v^6\,e^{6\,i\theta} \,. \end{align*}\end{split}\]

Note that the original expression for \(t(u,v)\) arises only in the singular limit \(\alpha\to 0\).

# Check of the above:
def reduce_multiply(any_list):
    return reduce(lambda x, y: x * y, any_list)


zs = [0, oo, 1, -1, I, -I]


def rotate(z, alpha, prephase=1):
    if z is oo:
        return -cos(alpha / 2) / sin(alpha / 2)
    return (cos(alpha / 2) * prephase * z + sin(alpha / 2)) / (
        -sin(alpha / 2) * prephase * z + cos(alpha / 2)
    )


alpha, p, w = symbols("alpha p w")


zps = [rotate(z, alpha) for z in zs]

# p <--> exp(i*theta)
tp = reduce_multiply([u - p * z * v for z in zps])
tp = simplify(tp.expand())
tp
\[\displaystyle p^{6} v^{6} - \frac{4 p^{5} u v^{5}}{\tan{\left(2 \alpha \right)}} - 5 p^{4} u^{2} v^{4} - 5 p^{2} u^{4} v^{2} + \frac{4 p u^{5} v}{\tan{\left(2 \alpha \right)}} + u^{6}\]

Icosahedral Edge Midpoints#

If we now consider the product of 5 such functions, representing the 30 edge midpoints of the icosahedron, with each inscribed octahedron rotated by \(e^{2\pi i/5}\) relative to one another, we find (introducing \(w\) as a stand-in for \(\tan(2\alpha)\))

\[\begin{split}\begin{aligned} T(u,v,w) &:= \prod_{k=1}^5 t'(u,v,w, e^{2\pi i k/5}) \\ &= u^{30} + \left(\frac{ 580\, w^4 + 1600 \, w^2 + 1024 }{w^5}\right)\left(u^{25}\,v^5-u^5\,v^{25}\right) \\ &- \left(\frac{ 6725 \,w^4 +11840 \,w^2 + 5120 }{w^4}\right)\left(u^{20}\,v^{10} + u^{10}\,v^{20}\right) + v^{30} \,. \end{aligned}\end{split}\]
# Check of the above:
factors = [tp.subs(p, e**k) for k in range(5)]
T_uvw = reduce_multiply(factors)

T_uvw = Poly(T_uvw, u)


new_coeffs = []


for coeff in T_uvw.all_coeffs():
    new_coeffs.append(e_simp(coeff, e))


T_uvw = Poly.from_list(new_coeffs, u).as_expr()
T_uvw
\[\displaystyle u^{30} + \frac{u^{25} \left(580 v^{5} \tan^{4}{\left(2 \alpha \right)} + 1600 v^{5} \tan^{2}{\left(2 \alpha \right)} + 1024 v^{5}\right)}{\tan^{5}{\left(2 \alpha \right)}} + \frac{u^{20} \left(- 6725 v^{10} \tan^{4}{\left(2 \alpha \right)} - 11840 v^{10} \tan^{2}{\left(2 \alpha \right)} - 5120 v^{10}\right)}{\tan^{4}{\left(2 \alpha \right)}} + \frac{u^{10} \left(- 6725 v^{20} \tan^{4}{\left(2 \alpha \right)} - 11840 v^{20} \tan^{2}{\left(2 \alpha \right)} - 5120 v^{20}\right)}{\tan^{4}{\left(2 \alpha \right)}} + \frac{u^{5} \left(- 580 v^{25} \tan^{4}{\left(2 \alpha \right)} - 1600 v^{25} \tan^{2}{\left(2 \alpha \right)} - 1024 v^{25}\right)}{\tan^{5}{\left(2 \alpha \right)}} + v^{30}\]

Lastly, recognizing that the icosahedral edge midpoints are located at

\[\cos\alpha = \frac1{2\,\sin(\pi/5)} \quad\leftrightarrow\quad \cot\alpha = \phi \,,\]

we find that \(w:=\tan(2\alpha) = 2\phi/(\phi^2-1) = 2\) and so

\[\begin{split}\begin{align*} T(u,v) := T(u,v, 2) &= 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{align*}\end{split}\]
# Check of the above:
phi = symbols("phi")
simplify((2 * phi / (phi**2 - 1)).subs(phi, (1 + sqrt(5)) / 2))
\[\displaystyle 2\]
T = simplify(T_uvw.subs(tan(2 * alpha), 2))
T
\[\displaystyle u^{30} + 522 u^{25} v^{5} - 10005 u^{20} v^{10} - 10005 u^{10} v^{20} - 522 u^{5} v^{25} + v^{30}\]

We also note in passing that \(T(u,v)\) can also be expressed as

\[T(u,v) = \prod_{k=1}^5 t_I(u\,\epsilon^{-k/2},v\,\epsilon^{k/2}) = \prod_{k=1}^5 \epsilon^{15-3k} \,t_I(u,v\,\epsilon^k)\]

where

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

The factor \(\epsilon^{15-3k}\) is engineered such that only non-negative integer powers of \(\epsilon\) appear in the resulting expression, which allows more straightforward simplification. (Recall that \(\epsilon^5=\epsilon^{15} = 1\).)

# Check of the above:
tI = (u**2 + v**2) * (u**2 - 2 * n * u * v - v**2) * (u**2 - 2 * m * u * v - v**2)
e_simp(
    reduce(
        lambda x, y: x * y,
        [e ** (15 - 3 * k) * tI.subs(v, v * e**k).expand() for k in range(5)],
    ),
    e,
)
\[\displaystyle u^{30} + 522 u^{25} v^{5} - 10005 u^{20} v^{10} - 10005 u^{10} v^{20} - 522 u^{5} v^{25} + v^{30}\]

Icosahedral Face Centres#

The icosahedral face centres nearest to the vertex \({\bf n}_5 = (0,0,+1)\) are given by

\[w_i = \epsilon^{i-1}\,w_1 \,,\qquad i\in\{1,2,3,4,5\}\]

where we choose \(w_1\) to be the (\(6\pi/5\))-rotated stereographic projection of \({\bf n}_3\):

\[w_1 = \epsilon^3\left(\frac{-\sin\alpha\,e^{-i\pi/5}}{\sqrt{3}/2 - \cos\alpha\,\cos\pi/5}\right) = \frac{2\,\sin\alpha}{\sqrt{3} + z_n\,\cos\alpha}\,.\]

(The choice of rotating by \(6\pi/5\) is made such that \(w_1\) is positive and real.)

# Check of the above:
alpha = acos(1 / 2 / sin(pi / 5))
w_1 = 2 * sin(alpha) / (sqrt(3) + z_n * cos(alpha))
w1_val = w_1.subs(z_n, -(1 + sqrt(5)) / 2).evalf()
w1_val
\[\displaystyle 2.95629520146761\]

Defining \(k:=\sqrt{3} / \cos\alpha\), and noting that

\[\begin{split}\begin{aligned} \frac1{\cos^2\alpha} = \frac2{1+\cos2\alpha} &= \frac{2\sqrt{5}}{1+\sqrt{5}} = \frac1{\phi}(2\,\phi-1) = 3+z_n \\ \implies k^2 &= 3(3+z_n) \,, \end{aligned}\end{split}\]

we find that we can simplify the expression for \(w_1\) as follows:

\[\begin{split}\begin{aligned} w_1 &= \frac{2\,\tan\alpha}{(k+z_n)} \times \frac{(k-z_n)}{(k-z_n)} = \frac{2\left(-1/z_n\right)(k-z_n)}{k^2-z_n^2}\\ &=\frac{-2(k-z_n)}{z_n\left[3(3+z_n)-z_n^2\right]} = \frac{z_n-k}{2(1+z_n)} \end{aligned}\end{split}\]

where, in deriving the last equality, we make liberal use of the identity \(z_n^2 = 1-z_n\) when rewriting the denominator.

Since, for arbitrary \(A\) and \(B\)

\[\frac1{Az_n+B} = \frac{A(1+z_n)-B}{A^2+AB-B^2} \,,\]

we can also write \(1/(1+z_n) = z_n\) (choosing \(A=B=1\) in the above), which gives

\[w_1 = \frac{1-(1+k)z_n}2 \,.\]

This expression, with only positive powers of \(k\) and \(z_n\), guarantee that successive powers of \(w_1\) can be easily simplified using the relations \(k^2=3(3+z_n)\), \(z_n^2 = 1-z_n\).

# Check of the above:
k = sqrt(3) / cos(alpha)
w_1 = (1 - (1 + k) * z_n) / 2
w_1.subs(z_n, -(1 + sqrt(5)) / 2).evalf().equals(w1_val)

True

w1h = (1 + (k - 1) * z_n) / 2
w1h.subs(z_n, -(1 + sqrt(5)) / 2).evalf()
\[\displaystyle -0.338261212717716\]
h1 = (-1 + (1 - k) * z_n) / 2
h1.subs(z_n, -(1 + sqrt(5)) / 2).evalf()
\[\displaystyle 0.338261212717716\]
((3 + sqrt(5) + sqrt(30 + 6 * sqrt(5))) / 4).evalf()
\[\displaystyle 2.95629520146761\]
((18 - 3 * (1 + sqrt(5))) * ((1 + sqrt(5)) ** 2)).expand()
\[\displaystyle 12 \sqrt{5} + 60\]

Using similar methods, one can easily show that the antipodal point \({\hat w}_1:= -1/w_1\) can be written

\[{\hat w}_1 = \frac{1+(k-1)z_n}2 \,.\]

Therefore, the product of polynomial factors with antipodal pairs, as in (1), are found to be

\[\begin{split}\begin{aligned} P(u,v,w_1,{\hat w}_1) &= u^{10} - (w_1^5+{\hat w}_1^5)u^5v^5-v^{10} \\ &= u^{10} - (64-100\,z_n)u^5v^5-v^{10} \,. \end{aligned}\end{split}\]

However, the above expression only accounts for half of the face centres of the icosahedron.

The other 10 face centres can be found by using the twofold symmetry \((IV)\) of the icosahedron rotating \(\hat w_1\) about \({\bf n}_2\) by \(\pi\):

\[w_6 := \frac{\hat w_1\,z_n+1}{\hat w_1-z_n} = \frac{2+z_n+k}2 \,.\]

From here, the other four face centres (related by symmetry \(V\)) can be easily determined by applying more powers of \(\epsilon\), and their antipodal pairs are given by

\[{\hat w}_6 := -\frac1{w_6} = \frac{2+z_n-k}2 \,.\]
w6 = (2 + z_n + k) / 2
w6.subs(z_n, -(1 + sqrt(5)) / 2).evalf()
\[\displaystyle 1.20905692653531\]
w6h = -1 / w6
w6h = (2 + z_n - k) / 2
w6h.subs(z_n, -(1 + sqrt(5)) / 2).evalf()
\[\displaystyle -0.827090915285202\]

Computing \(P(u,v,w_6, {\hat w}_6)\) as before, we find

\[\begin{split} \begin{aligned} P(u,v,w_6,{\hat w}_6) &= u^{10} - (w_6^5+{\hat w}_6^5)u^5v^5-v^{10} \\ &= u^{10} - [64+100(1+z_n)]u^5v^5-v^{10} \,. \end{aligned}\end{split}\]

Putting this all together, we find that the polynomial invariant resulting from the product of all 20 icosahedral face centres is given by

\[\begin{split}\begin{aligned} H(u,v) &= P(u,v,w_1,{\hat w}_1) \,P(u,v,w_6,{\hat w}_6) \\ &=u^{20} - 228\left(u^{15}v^5 - u^5v^{15}\right)+494\,u^{10}v^{10}+v^{20} \,. \end{aligned}\end{split}\]

(The above results from the general product

\[\begin{split}\begin{aligned} &\left(u^{10} - X\,u^5v^5-v^{10}\right)\left(u^{10} - Y\,u^5v^5-v^{10}\right) \\ &= u^{20}-(X+Y)\left(u^{15}v^5 - u^5v^{15}\right)+\left(XY-2\right)\,u^{10}v^{10}+v^{20} \end{aligned} \end{split}\]

where we take \(X=64-100\,z_n\), \(Y=64+100(1+z_n)\).

Lastly, we can verify that the polynomial invariants discussed above are interrelated as follows:

\[T^2-H^3 = 1728\,f^5 \,.\]
# Check of the above
(T**2 - H**3).equals(1728 * f**5)

True

(T**2).expand()
\[\displaystyle u^{60} + 1044 u^{55} v^{5} + 252474 u^{50} v^{10} - 10445220 u^{45} v^{15} + 100080015 u^{40} v^{20} - 10446264 u^{35} v^{25} + 199655084 u^{30} v^{30} + 10446264 u^{25} v^{35} + 100080015 u^{20} v^{40} + 10445220 u^{15} v^{45} + 252474 u^{10} v^{50} - 1044 u^{5} v^{55} + v^{60}\]

Edge Midpoints Redux#

An entirely orthogonal approach to the computation above is to consider directly the complex coordinates of each midpoint. We introduce here a calculus which can be used to simplify expressions involving these quantities.

We begin by recalling that the icosahedral edge midpoint on the negative real axis is given by

\[\begin{split}\begin{aligned} v_1 &:= -\frac{\sin\alpha}{1-\cos\alpha} = -\cot\left(\alpha / 2\right) = -\cot\alpha-\csc\alpha \\ &= z_n +j \end{aligned}\end{split}\]

where, in the last line, we define \(j:=-\csc\alpha\) (recall that \(z_n = -\cot\alpha\), as found earlier).

By directly computing the square of \(j\), we can verify that \(j^2=2-z_n\):

\[\begin{split}\begin{aligned} j^2 &= \frac1{\sin^2\alpha} = \frac2{1-\cos2\alpha} \\ &=\frac2{1-(1/\sqrt{5})} =\frac{2\sqrt{5}(1+\sqrt{5})}4 \\ &=(2\phi-1)\phi = 2+\phi \\ &= 2-z_n \,. \end{aligned}\end{split}\]
# Check of the above:
j = -csc(alpha)
(j**2).subs(alpha, 1 / 2 * acos(1 / sqrt(5))).evalf().equals(
    (2 + phi).subs(phi, (1 + sqrt(5)) / 2).evalf()
)

True

from sympy import cot


z_n = -cot(alpha)
z_n.subs(alpha, 1 / 2 * acos(1 / sqrt(5))).evalf()
\[\displaystyle -1.61803398874989\]
j = -csc(alpha)
(j).subs(alpha, 1 / 2 * acos(1 / sqrt(5))).evalf()
\[\displaystyle -1.90211303259031\]
(z_n + j).subs(alpha, 1 / 2 * acos(1 / sqrt(5))).evalf()
\[\displaystyle -3.5201470213402\]
1 / _
\[\displaystyle -0.284079043840412\]
(z_n - j).subs(alpha, 1 / 2 * acos(1 / sqrt(5))).evalf()
\[\displaystyle 0.284079043840412\]
(z_n - j).subs(z_n, -(1 + sqrt(5)) / 2).subs(j, -csc(1 / 2 * acos(1 / sqrt(5)))).evalf()
\[\displaystyle 0.284079043840412\]
from sympy import exp


e = exp(2 * I * pi / 5)
(e + e**4).evalf()
\[\displaystyle 0.618033988749895 + 7.0 \cdot 10^{-21} i\]

From here on out, it becomes quite straightforward to compute the invariant polynomial: for example, the antipodal point to \(v_1\) is

\[\begin{split}\begin{aligned} \hat v_1 &= -\frac1{v_1} = -\frac1{z_n+j}\times\frac{(z_n-j)}{(z_n-j)} \\ &=\frac{z_n-j}{2-z_n-z_n^2} = z_n -j \,. \end{aligned}\end{split}\]

The next rung of midpoints can be expressed in terms of \(v_1\) by applying a five-fold rotation, a reflection in \({\bf \hat n}_2\), and then three more five-fold rotations (yielding the edge midpoint lying on the positive real axis, which we refer to as \(v_9\)):

\[\begin{split}\begin{aligned} v_9 &=\epsilon^3 \times\left(\frac{\epsilon\,v_1\,z_n+1}{\epsilon\,v_1-z_n}\right) \\ &=\frac{\epsilon^3(\epsilon\,v_1\,z_n+1)(\epsilon^4\,v_1-z_n)}{v_1^2-z_n\,z_m\,v_1+z_n^2} \\ &=\frac{\epsilon^3(\epsilon\,v_1\,z_n+1)(\epsilon^4\,v_1-z_n)}{j(2\,z_n+1)-2\,z_n+4} \\ &=\frac{\epsilon^3(\epsilon\,v_1\,z_n+1)(\epsilon^4\,v_1-z_n)[j(2\,z_n+1)-4+2\,z_n]}{15\,z_n-10} \\ &=-\frac{\epsilon^3(\epsilon\,v_1\,z_n+1)(\epsilon^4\,v_1-z_n)[j(2\,z_n+1)-4+2\,z_n][15(1+z_n)+10]}{25} \\ &=(j-1)(1+z_n) \,. \end{aligned}\end{split}\]
# Check of the above:
v_1, j, z_n = symbols("v_1 j z_n")
num = (
    -(e**3)
    * (e * v_1 * z_n + 1)
    * (e**4 * v_1 - z_n)
    * (j * (2 * z_n + 1) - 4 + 2 * z_n)
    * (15 * (1 + z_n) + 10)
)
num = (
    num.expand()
    .subs(v_1**2, z_n**2 + 2 * z_n * j + 2 - z_n)
    .subs(v_1, j + z_n)
    .expand()
    .subs(j**2, 2 - z_n)
    .expand()
)
num
\[\displaystyle 30 j z_{n}^{5} e^{- \frac{2 i \pi}{5}} - 30 j z_{n}^{5} e^{- \frac{4 i \pi}{5}} + 95 j z_{n}^{4} e^{- \frac{2 i \pi}{5}} - 95 j z_{n}^{4} e^{- \frac{4 i \pi}{5}} + 30 j z_{n}^{3} e^{- \frac{4 i \pi}{5}} - 30 j z_{n}^{3} e^{\frac{4 i \pi}{5}} + 15 j z_{n}^{3} e^{- \frac{2 i \pi}{5}} + 160 j z_{n}^{2} e^{- \frac{4 i \pi}{5}} - 95 j z_{n}^{2} e^{\frac{4 i \pi}{5}} - 100 j z_{n}^{2} e^{- \frac{2 i \pi}{5}} - 15 j z_{n} e^{\frac{4 i \pi}{5}} - 25 j z_{n} e^{- \frac{4 i \pi}{5}} + 100 j e^{\frac{4 i \pi}{5}} + 30 z_{n}^{5} e^{- \frac{4 i \pi}{5}} + 50 z_{n}^{4} e^{- \frac{4 i \pi}{5}} - 15 z_{n}^{4} e^{- \frac{2 i \pi}{5}} + 5 z_{n}^{3} e^{- \frac{2 i \pi}{5}} - 150 z_{n}^{3} e^{- \frac{4 i \pi}{5}} + 50 z_{n}^{2} e^{- \frac{2 i \pi}{5}} + 15 z_{n}^{2} e^{\frac{4 i \pi}{5}} - 190 z_{n}^{2} e^{- \frac{4 i \pi}{5}} + 100 z_{n} e^{- \frac{4 i \pi}{5}} - 5 z_{n} e^{\frac{4 i \pi}{5}} - 50 e^{\frac{4 i \pi}{5}}\]
def z_simp(expr, z_n):
    """
    Reduces the polynomial in z_n using the identity z_n^2 = 1-z_n.
    """
    poly = Poly(expr, z_n)
    coeffs = poly.all_coeffs()
    while len(coeffs) > 2:
        top_factor = coeffs[0]
        coeffs[1] -= top_factor
        coeffs[2] += top_factor
        coeffs.pop(0)
    return Poly.from_list(coeffs, z_n).as_expr()


e_simp(z_simp(num, z_n).subs(z_n, e**2 + e**3), e).factor()

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/.local/share/mise/installs/python/3.13.7/lib/python3.13/site-packages/sympy/polys/polyutils.py:231, in _parallel_dict_from_expr_if_gens(exprs, opt)
    229         base, exp = decompose_power_rat(factor)
--> 231     monom[indices[base]] = exp
    232 except KeyError:

KeyError: exp(I*pi/5)

During handling of the above exception, another exception occurred:

PolynomialError                           Traceback (most recent call last)
Cell In[39], line 15
     11         coeffs.pop(0)
     12     return Poly.from_list(coeffs, z_n).as_expr()
---> 15 e_simp(z_simp(num, z_n).subs(z_n, e**2 + e**3), e).factor()

Cell In[7], line 2, in e_simp(expr, e, full, cleanup)
      1 def e_simp(expr, e, full=True, cleanup=True):
----> 2     result = Poly(expr, e)
      3     e_coeffs = reversed(result.all_coeffs())
      4     simp_coeffs = [0] * 5

File ~/.local/share/mise/installs/python/3.13.7/lib/python3.13/site-packages/sympy/polys/polytools.py:190, in Poly.__new__(cls, rep, *gens, **args)
    188     return cls._from_poly(rep, opt)
    189 else:
--> 190     return cls._from_expr(rep, opt)

File ~/.local/share/mise/installs/python/3.13.7/lib/python3.13/site-packages/sympy/polys/polytools.py:319, in Poly._from_expr(cls, rep, opt)
    316 @classmethod
    317 def _from_expr(cls, rep, opt):
    318     """Construct a polynomial from an expression. """
--> 319     rep, opt = _dict_from_expr(rep, opt)
    320     return cls._from_dict(rep, opt)

File ~/.local/share/mise/installs/python/3.13.7/lib/python3.13/site-packages/sympy/polys/polyutils.py:388, in _dict_from_expr(expr, opt)
    385         expr = expand_mul(expr)
    387 if opt.gens:
--> 388     rep, gens = _dict_from_expr_if_gens(expr, opt)
    389 else:
    390     rep, gens = _dict_from_expr_no_gens(expr, opt)

File ~/.local/share/mise/installs/python/3.13.7/lib/python3.13/site-packages/sympy/polys/polyutils.py:327, in _dict_from_expr_if_gens(expr, opt)
    325 def _dict_from_expr_if_gens(expr, opt):
    326     """Transform an expression into a multinomial form given generators. """
--> 327     (poly,), gens = _parallel_dict_from_expr_if_gens((expr,), opt)
    328     return poly, gens

File ~/.local/share/mise/installs/python/3.13.7/lib/python3.13/site-packages/sympy/polys/polyutils.py:236, in _parallel_dict_from_expr_if_gens(exprs, opt)
    234                 coeff.append(factor)
    235             else:
--> 236                 raise PolynomialError("%s contains an element of "
    237                                       "the set of generators." % factor)
    239 monom = tuple(monom)
    241 if monom in poly:

PolynomialError: exp(2*I*pi/5) contains an element of the set of generators.

We can also compute the antipodal point \(\hat v_9\) by inverting the expression above:

\[\begin{split}\begin{aligned} \hat v_9 &= -\frac1{v_9}= -\frac1{(j-1)(1+z_n)} \\ &=-\frac{(1+j)z_n}{1-z_n} \\ &=-\frac{(j+1)z_n(2+z_n)}{(1-z_n)(2+z_n)} \\ &=-(1+j)(1+z_n) \,. \end{aligned}\end{split}\]

Since the three pairs of antipodal points \((i,-i)\), \((v_1, \hat v_1)\) and \((v_9, \hat v_9)\) compromise one of the five octohedra, we now have a second way to compute the polynomial invariant \(t_I(u,v)\): since pairs of real antipodal points \((z, \hat z)\) contribute a factor \((u^2 - (z+\hat z)uv - v^2)\), it’s sufficient to compute the sums

\[\begin{split}\begin{aligned} v_1 + \hat v_1 &= 2\,z_n \\ v_9 + \hat v_9 &= - 2(1+z_n) = 2\,z_m \,. \end{aligned}\end{split}\]

Therefore, the factor for this octohedron is given as before:

\[t_0(u,v) = (u^2+v^2)(u^2-2z_n\,uv-v^2)(u^2-2z_m\,uv-v^2) \,.\]

The contributions to \(T(u,v)\) from the other 24 edge midpoints can be determined by applying successive five-fold rotations about the \(\zeta\) axis:

\[\begin{split}\begin{aligned} t_r(u,v) &= \epsilon^{15 -3r}t_0(u,v\,\epsilon^{r}) \\ &= \epsilon^{15 -3r} u^6 + 2\,\epsilon^{10 -2r}\,u^5\,v - 5\,\epsilon^{5-r}\,u^4\,v^2 - 5\,\epsilon^r \,u^2\,v^4 -2\,\epsilon^{2r}\,u\,v^5 + \epsilon^{3r}\,v^6\,. \end{aligned}\end{split}\]
# Check of the above
from utils import 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))
t_r
\[\displaystyle e^{15 - 3 r} u^{6} + 2 e^{15 - 2 r} u^{5} v - 5 e^{15 - r} u^{4} v^{2} - 5 e^{r + 15} u^{2} v^{4} - 2 e^{2 r + 15} u v^{5} + e^{3 r + 15} v^{6}\]
e_simp(t_r.subs(r, 3).expand().coeff(v**6), e).subs(e**3, -1 - e - e**2 - e**4).equals(
    e**4
)

True

For completeness, we also record for later use the corresponding Hessian of the octahedral vertex function \(t_r(u,v)\), representing the polynomial invariant for the faces of each of the five inscribed octahedra:

\[\begin{split}\begin{aligned} W_r(u,v) := -\frac{1}{400}\cdot \mathrm{Hess}[t_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}\end{split}\]
# Check of the above
from utils import hessian

W_r = simplify(hessian(t_r, u, v).expand() / (-400))
W_r
\[\displaystyle e^{30 - 4 r} u^{8} - e^{30 - 3 r} u^{7} v + 7 e^{30 - 2 r} u^{6} v^{2} + 7 e^{30 - r} u^{5} v^{3} - 7 e^{r + 30} u^{3} v^{5} + 7 e^{2 r + 30} u^{2} v^{6} + e^{3 r + 30} u v^{7} + e^{4 r + 30} v^{8}\]
e_simp(reduce(lambda x, y: x * y, [t_r.subs(r, x) for x in range(5)]), e)
\[\displaystyle u^{30} + 522 u^{25} v^{5} - 10005 u^{20} v^{10} - 10005 u^{10} v^{20} - 522 u^{5} v^{25} + v^{30}\]
z = symbols("z")


expr = (f**5).subs(u, z).subs(v, 1)
expr
\[\displaystyle z^{5} \left(z^{10} + 11 z^{5} - 1\right)^{5}\]
from sympy import series

from utils import j_simp, k_simp, z_simp

z_n, j = symbols("z_n j")

z_simp(j_simp(series(expr, z, z_n - j, n=2).removeO(), j, z_n), z_n)

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polyutils.py:229, in _parallel_dict_from_expr_if_gens(exprs, opt)
    227         base, exp = decompose_power_rat(factor)
--> 229     monom[indices[base]] = exp
    230 except KeyError:

KeyError: 1/(j**11 - 11*j**10*z_n + 55*j**9*z_n**2 - 165*j**8*z_n**3 + 330*j**7*z_n**4 - 462*j**6*z_n**5 - 11*j**6 + 462*j**5*z_n**6 + 66*j**5*z_n - 330*j**4*z_n**7 - 165*j**4*z_n**2 + 165*j**3*z_n**8 + 220*j**3*z_n**3 - 55*j**2*z_n**9 - 165*j**2*z_n**4 + 11*j*z_n**10 + 66*j*z_n**5 - j - z_n**11 - 11*z_n**6 + z_n)

During handling of the above exception, another exception occurred:

PolynomialError                           Traceback (most recent call last)
Cell In[54], line 7
      3 from utils import j_simp, k_simp, z_simp
      5 z_n, j = symbols("z_n j")
----> 7 z_simp(j_simp(series(expr, z, z_n - j, n=2).removeO(), j, z_n), z_n)

File ~/Documents/coding/ag_notes/notes/utils.py:137, in j_simp(expr, j, z_n)
    136 def j_simp(expr, j, z_n):
--> 137     poly = Poly(expr, j)
    138     coeffs = poly.all_coeffs()
    139     while len(coeffs) > 2:

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polytools.py:186, in Poly.__new__(cls, rep, *gens, **args)
    184     return cls._from_poly(rep, opt)
    185 else:
--> 186     return cls._from_expr(rep, opt)

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polytools.py:315, in Poly._from_expr(cls, rep, opt)
    312 @classmethod
    313 def _from_expr(cls, rep, opt):
    314     """Construct a polynomial from an expression. """
--> 315     rep, opt = _dict_from_expr(rep, opt)
    316     return cls._from_dict(rep, opt)

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polyutils.py:386, in _dict_from_expr(expr, opt)
    383         expr = expand_mul(expr)
    385 if opt.gens:
--> 386     rep, gens = _dict_from_expr_if_gens(expr, opt)
    387 else:
    388     rep, gens = _dict_from_expr_no_gens(expr, opt)

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polyutils.py:325, in _dict_from_expr_if_gens(expr, opt)
    323 def _dict_from_expr_if_gens(expr, opt):
    324     """Transform an expression into a multinomial form given generators. """
--> 325     (poly,), gens = _parallel_dict_from_expr_if_gens((expr,), opt)
    326     return poly, gens

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polyutils.py:234, in _parallel_dict_from_expr_if_gens(exprs, opt)
    232                 coeff.append(factor)
    233             else:
--> 234                 raise PolynomialError("%s contains an element of "
    235                                       "the set of generators." % factor)
    237 monom = tuple(monom)
    239 if monom in poly:

PolynomialError: 1/(j**11 - 11*j**10*z_n + 55*j**9*z_n**2 - 165*j**8*z_n**3 + 330*j**7*z_n**4 - 462*j**6*z_n**5 - 11*j**6 + 462*j**5*z_n**6 + 66*j**5*z_n - 330*j**4*z_n**7 - 165*j**4*z_n**2 + 165*j**3*z_n**8 + 220*j**3*z_n**3 - 55*j**2*z_n**9 - 165*j**2*z_n**4 + 11*j*z_n**10 + 66*j*z_n**5 - j - z_n**11 - 11*z_n**6 + z_n) contains an element of the set of generators.

z_simp(j_simp(series(expr, z, z_n - j, n=2).expand().removeO(), j, z_n), z_n)

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polyutils.py:229, in _parallel_dict_from_expr_if_gens(exprs, opt)
    227         base, exp = decompose_power_rat(factor)
--> 229     monom[indices[base]] = exp
    230 except KeyError:

KeyError: 1/(j**11 - 11*j**10*z_n + 55*j**9*z_n**2 - 165*j**8*z_n**3 + 330*j**7*z_n**4 - 462*j**6*z_n**5 - 11*j**6 + 462*j**5*z_n**6 + 66*j**5*z_n - 330*j**4*z_n**7 - 165*j**4*z_n**2 + 165*j**3*z_n**8 + 220*j**3*z_n**3 - 55*j**2*z_n**9 - 165*j**2*z_n**4 + 11*j*z_n**10 + 66*j*z_n**5 - j - z_n**11 - 11*z_n**6 + z_n)

During handling of the above exception, another exception occurred:

PolynomialError                           Traceback (most recent call last)
Cell In[55], line 1
----> 1 z_simp(j_simp(series(expr, z, z_n - j, n=2).expand().removeO(), j, z_n), z_n)

File ~/Documents/coding/ag_notes/notes/utils.py:137, in j_simp(expr, j, z_n)
    136 def j_simp(expr, j, z_n):
--> 137     poly = Poly(expr, j)
    138     coeffs = poly.all_coeffs()
    139     while len(coeffs) > 2:

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polytools.py:186, in Poly.__new__(cls, rep, *gens, **args)
    184     return cls._from_poly(rep, opt)
    185 else:
--> 186     return cls._from_expr(rep, opt)

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polytools.py:315, in Poly._from_expr(cls, rep, opt)
    312 @classmethod
    313 def _from_expr(cls, rep, opt):
    314     """Construct a polynomial from an expression. """
--> 315     rep, opt = _dict_from_expr(rep, opt)
    316     return cls._from_dict(rep, opt)

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polyutils.py:386, in _dict_from_expr(expr, opt)
    383         expr = expand_mul(expr)
    385 if opt.gens:
--> 386     rep, gens = _dict_from_expr_if_gens(expr, opt)
    387 else:
    388     rep, gens = _dict_from_expr_no_gens(expr, opt)

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polyutils.py:325, in _dict_from_expr_if_gens(expr, opt)
    323 def _dict_from_expr_if_gens(expr, opt):
    324     """Transform an expression into a multinomial form given generators. """
--> 325     (poly,), gens = _parallel_dict_from_expr_if_gens((expr,), opt)
    326     return poly, gens

File ~/Documents/coding/ag_notes/.venv/lib/python3.13/site-packages/sympy/polys/polyutils.py:234, in _parallel_dict_from_expr_if_gens(exprs, opt)
    232                 coeff.append(factor)
    233             else:
--> 234                 raise PolynomialError("%s contains an element of "
    235                                       "the set of generators." % factor)
    237 monom = tuple(monom)
    239 if monom in poly:

PolynomialError: 1/(j**11 - 11*j**10*z_n + 55*j**9*z_n**2 - 165*j**8*z_n**3 + 330*j**7*z_n**4 - 462*j**6*z_n**5 - 11*j**6 + 462*j**5*z_n**6 + 66*j**5*z_n - 330*j**4*z_n**7 - 165*j**4*z_n**2 + 165*j**3*z_n**8 + 220*j**3*z_n**3 - 55*j**2*z_n**9 - 165*j**2*z_n**4 + 11*j*z_n**10 + 66*j*z_n**5 - j - z_n**11 - 11*z_n**6 + z_n) contains an element of the set of generators.

expr2 = (H**3).subs(u, z).subs(v, 1)
expr2
\[\displaystyle \left(z^{20} - 228 z^{15} + 494 z^{10} + 228 z^{5} + 1\right)^{3}\]
from sympy import series

from utils import k_simp, z_simp, j_simp

z_n, j = symbols("z_n j")

z_simp(j_simp(series(-expr2 / 1728, z, z_n - j, n=1).removeO(), j, z_n), z_n)
\[\displaystyle - 75413619148164581909179687500 j + z_{n} \left(122021798996369555969238281250 j + 232099254131109928588867187500\right) - 143445227816524921356201171875\]
from sympy import exp


edge_val = (
    (z_n + j)
    .subs(j, sqrt(2 - z_n))
    .subs(z_n, e**2 + e**3)
    .subs(e, exp(2 * pi * I / 5))
    .evalf()
)
edge_val
\[\displaystyle 0.284079043840412 + 5.0 \cdot 10^{-26} i\]
from sympy import exp


edge_val = (
    (z_n + j)
    .subs(j, sqrt(2 - z_n))
    .subs(z_n, e**2 + e**3)
    .subs(e, exp(2 * pi * I / 5))
    .evalf()
)
edge_val
\[\displaystyle 0.284079043840412 + 5.0 \cdot 10^{-26} i\]
e_simp(e**8, e)
\[\displaystyle e^{3}\]
series(-expr2 / (1728 / expr), z, 1, n=1)
\[\displaystyle - \frac{307063701824}{27} + O\left(z - 1; z\rightarrow 1\right)\]