Subject: Genus 2 CM construction for class number 20,016

Dear number theorists,

We are pleased to announce the successful construction of a genus 2
hyperelliptic Jacobian over a 128-bit prime field, using the complex
multiplication method, for an input CM field of class number 20,016.

The computation was performed using the complex analytic method, using
the implementation [3], which is described in [2]. Important
characteristics of this implementation is the ability to compute
theta-constants in quasi-linear complexity in the input size. Earlier
work by Dupont was used as a starting point for this project [1].

We are aware of previous existing computation using this method. To our
knowledge, the largest of these is for a CM field of class number 530 [4].

The input of the computation is the CM field K given by the defining
equation x^4+1357*x^2+3299s. The class group of K has structure
Z/2*Z/4*Z/5004. The real quadratic subfield K0 is Q(sqrt(1828253)), whose
class number is 2. The Shimura class group C(K) is isomorphic to
Z/2*Z/2*Z/5004, and thus has size 20,016. The size of the Shimura class
group corresponds to the degree of the irreducible factors (over K0r, the
real quadratic subfield of the reflex field Kr of K) of the Igusa class
polynomials.

The most important part of the computation was the computation of the
Igusa class polynomials (H1, \hat H2, \hat H3). These were computed
relative to the choice of invariants proposed in [5].

Computing (H1, \hat H2, \hat H3) required computing complex
approximations of 10,008 pairs of complex conjugate triples of Igusa
invariants, associated to a set of period matrices.

The breakdown of the computation times goes as follows.
 - Computation of the period matrices defined over the reflex field of K:
    388 s
 - First steps of Newton lifting up to precision roughly 2,000,000 bits:
    610 s per period matrix (1,700 hours total).
 - Last two steps of Newton lifting up to precision roughly 8,000,000 bits:
    3,040 s per period matrix. 8,451 hours total.  2 days wall-clock time.
 - Computation of floating point polynomials H1, \hat H2, \hat H3:
    3 days wall-clock time.
 - Recognition of polynomial coefficients as elements of K0r:
    980 s per coefficient. 16,350 hours total. 2 days wall-clock time.

All computation steps have been performed in parallel over up to several
hundreds of CPU cores, to the extent allowed by the constraints of the
various steps of the computation. A variety of hardware platforms have
been used, the most common being Intel Nehalem-based and Sandy
Bridge-based processors with 2 GHz clock frequency (see [3] for details).

The uncompressed storage size of the three resulting polynomials in base
10 is about 95 GB. The common denominator of the coefficients of H1 has
8884 distinct prime factors, the largest one being 1506803839.

In order to validate the computation, the construction of a genus 2 curve
was done. We picked the prime number p=2^128+5399685, and the following
Weil number:
    pi:=2587584949432298*y^3 + 598749326588980*y^2 + 3489110163205995872*y - 17626367557116479015;
    where y is a root of t^4+1357*t^2+3299;

The evaluation of the Igusa class polynomials at a single root of H1 gave
a triple of Igusa invariants, to which we applied Mestre's algorithm for
reconstructing a curve from its invariants. The hyperelliptic curve
obtained has the following equation over GF(p):

    y^2 = 329105434147215182703081697774190891717*x^5 + 219357712933218699650940059644263138156*x^4 + 94773520721686083389380651745963315116*x^3 + 13612280714446818104030347122109215819*x^2 + 224591198286067822213326173663420732292*x + 62350272396394045327709463978232206155;

The complex multiplication theory predicts the cardinal of the Jacobian,
which we have been able to validate as being equal to:
    Norm(1-pi) = 115792089237316195448115714075359184294817543055097597784192266245863261539824

This validation of the computation took roughly 20 minutes, mostly
dominated by the input-output time required to reduce the class
polynomial modulo p (part of the computation was offloaded to a nearby
cluster).

All computations have been performed by cmh-1.0, which is programmed in C
and pari/gp, and can be freely downloaded.

Questions regarding characteristics of the data computed (prime numbers
appearing in the denominator, etc) may be directed to the authors.

Andreas Enge,
Emmanuel Thomé.


[1] R. Dupont. Moyenne arithmético-géométrique, suites de Borchardt et
    applications. Thèse, École Polytechnique, 2006.
    http://www.lix.polytechnique.fr/Labo/Regis.Dupont/these_soutenance.pdf
[2] A. Enge, E. Thomé, Computing class polynomials for abelian surfaces.
    To appear in Experimental Mathematics, 2014.
    http://hal.inria.fr/hal-00823745/
[3] A. Enge, E. Thomé, cmh: Computation of Igusa Class Polynomials,
    version 1.0, 2014.
    http://cmh.gforge.inria.fr/
[4] T. Houtmann, Complex multiplication in genus 2, 2007.
    http://tiresias.dyndns-at-home.com/cm/
[5] M. Streng, Complex multiplication of abelian surfaces. PhD thesis,
    Universiteit Leiden, 2010.
    http://hdl.handle.net/1887/15572


Verification script in magma:

p:=340282366920938463463374607431773611141;
A:=1357;
B:=3299;
K:=GF(p);
KP<x>:=PolynomialRing(GF(p));
C:=HyperellipticCurve(329105434147215182703081697774190891717*x^5 +
219357712933218699650940059644263138156*x^4 +
94773520721686083389380651745963315116*x^3 +
13612280714446818104030347122109215819*x^2 +
224591198286067822213326173663420732292*x +
62350272396394045327709463978232206155);
J:=Jacobian(C);

ZP<t>:=PolynomialRing(Integers());
K<y>:=NumberField(t^4+A*t^2+B);
pi:=2587584949432298*y^3 + 598749326588980*y^2 + 3489110163205995872*y
- 17626367557116479015;
assert Norm(pi) eq p^2;
assert IsZero(Integers()!Norm(1-pi)*Random(J));
chi:=CharacteristicPolynomial(pi);
s1:=-Coefficient(chi,3);
s2:=Coefficient(chi,2);
assert Coefficient(chi,1) eq -p*s1;
assert Coefficient(chi,0) eq p^2;
/* The curve C has p+1-s1 points over GF(p) */

assert IsZero(Discriminant(ZP!chi) mod Discriminant(MaximalOrder(K)));

J`EulerFactor:=ReciprocalPolynomial(chi);
J`Order:=Norm(1-pi);

    /* base-extend to GF(p^3), for a check */
    n:=3;
    Jn:=BaseChange(J,n);
    ZP1:=PolynomialRing(Integers());
    ZP2:=PolynomialRing(ZP1);
    E1:=ZP1!J`EulerFactor;
    E:=Evaluate(E1,ZP2.1);
    P:=ZP2.1^n-ZP1.1;
    En:=Resultant(E,P);
    assert(IsZero(Evaluate(En,1)*Random(Jn)));
    Jn`EulerFactor:=En;
    Jn`Order:=Evaluate(En,1);