## Hilbert Class Polynomials

These require floating point arithmetic to compute [Cohen p415, BSS p151]. Routines for complex numbers need to be written. For multiplication, using a trick similar to Karatsuba multiplication

$(a+b i)(c+d i) = a c - b d + ((a+b)(c+d)-a c-b d)i$

saves one multiplication. Similarly, for squaring, one should use

$(a+b i)^2 = (a+b)(a-b) + 2a b i$

if multiplication is sufficiently less efficient than addition or subtraction.

To minimize the possbility of underflow, overflow or precision loss, the complex modulus of $a+bi$ should be computed as [Numerical Recipes, 5.4]

$|a+i b|= \begin{cases} |a|\sqrt{1 + (b/a)^2} & |a| \ge |b| \\ |b|\sqrt{1 + (a/b)^2} & |a| \lt |b| \end{cases}$

Similarly, division should be implemented as

$\frac{a+i b}{c + i d} = \begin{cases} \frac{(a + b(d/c)) + i(b - a(d/c))}{c+d(d/c)} & |c| \ge |d| \\ \frac{(a(c/d) + b) + i(b(c/d) - a)}{c(c/d)+d} & |c| \lt |d| \end{cases}$

We will need to compute the function $j(\tau)$ for $\tau$ in the upper complex plane. Set $q = e^{2i \pi \tau}$. Define

$\Delta(\tau)=q{({ 1 + \sum_{n\ge 1} (-1)^n {({ q^{n(3n-2)/2} + q^{n(3n+1)/2} })} })}^{24}$

Then compute

$j(\tau) = \frac{(256f(\tau)+1)^3}{f(\tau)}$

where

$f(\tau)=\frac{\Delta(2\tau)}{\Delta(\tau)}$

The Hilbert class polynomial for the discriminant $-D$ is given by $H_D(x) = \prod (x-j(\alpha))$ where $\alpha$ runs over all complex numbers such that

$\alpha = \frac{-b+\sqrt{-D}}{2a}$

where $a x^2 + b x y + c y^2$ is a primitive reduced positive definite binary quadratic form of discriminant $-D$. In other words, $b^2 - 4a c = -D$, $|b|\le a\le \sqrt{|D|/3}$, $a\le c$, $\gcd(a,b,c)=1$ and if $|b|=a$ or $a=c$ then $b\ge 0$.

Thus $H_D(x)$ can be computed by the following procedure. Let $P$ be a polynomial variable. Assume $D$ is $0$ or $-1$ $mod 4$.

1. Set $P = 1, b = D mod 2, B = {\lfloor \sqrt{|D| / 3} \rfloor}$

2. Set $t = (b^2 + D)/4$ and $a = max(b, 1)$

3. If $a$ does not divide $t$ go to the next step. Otherwise compute $j = j((-b+\sqrt{-D})/(2a))$. If $a = b$ or $a^2 = t$ or $b = 0$ set $P = P(X-j)$ otherwise set $P = P(X^2 - 2 \Re (j) X + |j|^2)$.

4. Set $a = a + 1$. If $a^2 \le t$ go to the previous step.

5. Set $b = b + 2$. If $b \le B$ go to step 2, otherwise round the coefficients of $P$ to the nearest integer and output $P$.