Modulo Errors

Update 4/v/07: I’ve switched from Cantor’s definitions for a curve of the form y²=f(x) to a more general form, following the notation of a paper by Tanja Lange; that also describes many efficiency gains for these calculations, none of which I have yet adopted… I’m also implementing these procedures in SAGE, which seems a more natural environment. So consider all genus 2 stuff as work in progress!

jac is an implementation of the group law on the jacobian of a genus 2 hyperelliptic curve over a finite field, to work with the generic_group procedures described previously. Standard version is for Maple 10; you can also get a version for Maple 9, but this may not be updated as frequently.

An arbitrary divisor D is now either a list [a(u),b(u)] or the identity element zero. Addition of two such divisors D,E is given by g2JacGroupLaw(D,E) whilst g2JacMinus(D) gives the inverse. So these functions can be used as arguments for ncopies and so on. To set up the worksheet, specify a characteristic p; a degree five monic squarefree polynomial f(u) and a polynomial h(u) of degree at most 2. Only rational divisors and prime fields seem to work: working mod p generates sufficiently ugly Maple code to discourage me from trying extension fields there.

Also included are a couple of ways to get random divisors to compute with. randomDiv is incredibly slow as it naively tests random choices of a (monic quadratic) and b (linear) for suitability (that is, a dividing b^2+bh-f). randPoint(f,h) is smarter (transforming to y²=g(x), using the legendre symbol to test random choices of x for square g(x) then finding a root and transforming back to a suitable y) and of course you can combine two points into a weight 2 divisor using the group law.

Order computation, even with BSGS, becomes very slow for less than staggering values of p: this is of course the point cryptographically! For instance, it took my university workstation about 11 hours to find the order of a randomly constructed divisor from a curve over a field with around a hundred-thousand elements.

ella is an update of gla, the elliptic curve group law procedures, to work with the generic_group procedures described previously.

An arbitrary point P is now either a list [x,y] or the identity element zero. Addition of two such points P,Q is given by ellGroupLaw(P,Q) whilst ellMinus(P) gives the inverse. So these functions can be used as arguments for ncopies and so on. Setting up the curve is as described for gla: specify variables a_1 through a_6 and optionally set workModM and a modulus M for computation over prime fields. The getQuantities function is also included for convenience; order calculations and functions like the old mnadd should make use of the generic_group procedures and thus are omitted.

I’m keeping gla available for use with torsion_tools, most of the functions of the latter depend strongly on the underlying group being that of an elliptic curve so cannot be translated to generic form, so I can’t be bothered to update the notation from x,y to [x,y] !

Overview

For the last month or so I’ve been working with hyperelliptic rather than elliptic curves. This requires working in the Jacobian of the curve rather than the curve itself, but the basic computational tasks are of course the same: adding arbitrary points, multiplication by m, finding orders. I wrote Maple code for all of these in the elliptic curve case, but the simplicity of both the group law and the examples I worked on meant I could get away with brute force computation of these. But for higher genus (even genus 2), this becomes crippling very early on, so I needed some smarter methods. Having discovered that you can pass functions as arguments in Maple, it seemed best to divorce the group methods from the specifics of the group being computed with. To that end, I present generic_group, a set of Maple procedures for generic group algorithms.

A generic algorithm is one which only performs the group operation, inversion, or testing for equality. To that end, these procedures may take as arguments: elements; a function grouplaw that performs the group operation on an input of two elements; or another function groupminus that inverts its input element. Since I’m usually working with additive groups, the neutral element is given as zero, and so the supplied functions must be able to handle this.

The idea is that other sets of procedures will specify those functions: for instance, I’ve been writing procedures for the group of divisors of an algebraic curve: as a special case I have a function g2JacGroupLaw which will sum two rational divisors from the Jacobian of a genus 2 hyperelliptic curve. Rather than writing a multiplication by m function for such divisors, I just plug that function into the generic procedure for finding n copies of an element. I’ve updated my earlier elliptic curve procedures to this format as well.

The Procedures

Multiplication by n (ncopies and brutencopies)

A call of ncopies(n,g,grouplaw,groupminus) computes [n]g; it catches the special cases of n=0 (giving zero), n=1 (giving g) and n negative (giving [|n|](-g)).
This procedure uses around 2ceil(log₂n) instances of the group law by exploiting fast exponentiation, rather than naively adding successive terms. If out of boredom or to compare performance you want naive addition, you can use brutencopies(n,G,grouplaw,groupminus) instead.

Discrete Logarithm Problem (BSGSDL)

In the previous post I described the Baby Step, Giant Step algorithm for computing discrete logarithms, in the case of unknown order of the generator (Terr’s variant). BSGSDL(g,h,grouplaw) will use this to return an integer t such that [t]g=h; this assumes that there is such a t, i.e., that h is in <g>. beware that if it isn’t, the procedure will wander off to deep space and never terminate.

Order finding (BSGSOrder and bruteOrder)

Terr’s variant can be used to establish the order of an element by solving the DLP with h=zero; however you have to check that g itself is not the identity first. So BSGSOrder(g,grouplaw) will find the least t such that [t]g=zero. If you want to try this by checking successive multiples of g, you can use bruteOrder(g,grouplaw) instead; again, knowing nothing about the group neither of these procedures knows when to give up so can churn forever if g is of infinite order.

I’ve spent the last few weeks getting my head around the SEA algorithm- a series of techniques for computing cardinality of elliptic curves over finite fields. The basic Schoof algorithm (the S of SEA) is based on the Weil conjectures, from which to find the cardinality it suffices to determine the trace of Frobenius, t. This can be shown to fall within certain bounds and thus (by the chinese remainder theorem) the problem can be solved by establishing t modulo l for a series of smaller primes l. The contributions by Elkies and Atkin make these subcalculations feasible, but trade computational complexity for mathematical difficulty, and it’s these I’ve been grappling with recently- particularly the Elkies approach.

Fortunately there is an excellent book - Blake, Seroussi and Smart’s Elliptic Curves in Cryptography which covers this topic thoroughly. As usual, attempting to implement the procedures myself has uncovered (and repaired) many gaps in my comprehension of the algorithms involved. I’m finally getting correct answers, albeit for small examples, so I’ve written up some summary notes on the book and a worked example in Maple. It still requires a high degree of interaction to get to the final result; my understanding is that SAGE (on linux) implements the SEA algorithm, so it doesn’t seem worth trying to tie it all together into a standalone Maple procedure- this would be a programming challenge, but wouldn’t add anything to my mathematical understanding!

Notes: [SEA.dvi]
Appendix (Maple procedures): [SEA_tools.mpl]
Example (Maple worksheet): [SEA.mw]

Ask SAGE for the cardinality (that is, number of points) of an elliptic curve over a finite field and, unless you happen to have a prime field, it’ll warn you that it’s resorting to the “very, very stupid” algorithm of testing every point. Can we do better?

Background: The Zeta function

Recall that the zeta function of a curve encodes information about the number of points over an infinite family of fields. For a curve C over the field F_q This is defined as

Thus, we can recover the cardinality over F_q^r by differentiating r times with respect to T, evaluating for T=0 and dividing through by (r-1)!. This is all well and good, provided that computing the zeta function is less work than simply testing points in the field directly. For this, we need a couple of results:

Weil Conjectures for curve C over F_q
We have

such that

with the α_i algebraic integers such that |α_i|=q^1/2.

For such α_i we then have

For F_q to be a field, q=p^r for some prime p. Thus we need only compute the zeta function over F_p to be able to retrieve #C(F_q) or indeed any higher #C(F_q^r). The question is, therefore, whether computing the zeta function for F_p is any easier than the naive point search. With Weil’s result, we can see that knowledge of #C(F_p),…,#C(F_p^2g) would, via some Taylor series trickery, be enough to determine P(T) and hence Z(T). If q is a high power of p, or we are interested in high powers of q, then this is already progress. But if q was itself a (potentially very large) prime, this approach would be useless for finding #C(F_q) since it would require knowledge of #C(F_q) (!) and at least #C(F_q²) too! Fortunately SAGE does not complain about searching over prime fields- there are a couple of algorithms - so we can restrict our attention to the original problem, that of composite q.

Elliptic Curves

Let us consider then the problem of finding #E(F_q^r) where q=p^s for a prime p with rs≥2 and E an Elliptic curve with coefficients from F_p. Then by the above we need only find #E(F_p) and #E(F_p²), which by brute force is rarely more work than finding #E(F_q^r)=#E(F_p^rs). But, since the genus is 1 the various results allow us to cut down our workload still further:

Weil Conjectures for Elliptic curve E over F_p

We have

such that

with the α_i algebraic integers such that |α_i|=p^1/2.

For such α_i we then have

It follows that α₁, α₂ are conjugates, so over F_p we have

and

So,

meaning that the zeta function can be constructed from just #E(F_p) as it is simply:

This then allows us to recover #E(F_q^r)=#E(F_p^rs) as desired, by repeated differentiation to recover the rs‘th coefficient.

Comparison

For an elliptic curve over a finite field, using SAGE to calculate the cardinality over a field of p^r elements by the “very, very stupid” algorithm rapidly gets impractical. For instance, for a given elliptic curve y²=x³+Ax+b with A,B from F₅ determining the cardinality takes about 0.04s over F₅; about a second over F₂₅, around 3.7s for F₁₂₅ and a tedious 19s over F₆₂₅. Implementing the approach above, asking for the cardinality of three such curves over F₆₂₅ is rapid enough in Maple to not register on its timer. This is despite my program lazily using brute force for the #E(F₅) calculations! (The same machine, with a 2.6ghz celeron processor, was used for each run; with Maple on Windows XP and SAGE in Xubuntu Linux; SAGE timings were for CPU rather than wall time.)

Smarter Still

In fact, if we only wish to determine a single #E(F_q^r) we can sidestep the construction of the zeta function by working with the α_p, since

Thus if α_p=a+bi then

and

So

will suffice, and again this requires only knowledge of #E(F_p).

Maple Source

The file ellCount contains Maple procedures for finding the cardinality of elliptic curves over finite fields (with coefficients from the prime subfield). Simply put, CountPoints(F,q) will find the number of points over the field of q elements for an Elliptic curve F=0. For instance, CountPoints(y^2-x^3-x-2,625) gives the cardinality of E: y²=x³+x+2 over F_5⁴ (should be 640). The route taken is to find points over F₅, determine α then directly calculate from the formula in the previous section. You can retrieve the zeta function with ellZeta(F,q), or ellZetap(F,p), if you know that q is prime (which saves Maple trying to determine the prime factor). Thus the number of points over F_q^r can be retrieved with getZetaCoeffr(Z,r) which performs the appropriate differentiation and scaling; this is useful to avoid having to calculate the number of solutions over the prime field, or indeed determine the prime, every time. ZetaCountPoints(F,q) behaves as CountPoints except it follows this alternative route.

The procedure I described in my previous post for computing types of zeta functions is informative the first time you try it and tedious thereafter; thus I’ve cobbled together some maple code to do the job for some simple cases.

A call to totalTrace(d,t) will scurry off and determine all degree d polynomials of trace t, provided neither your chosen trace or degree are too high (data tables only exist up to a certain point, and I wasn’t patient enough to implement much of what is known, either!). If you’re confident that the trace is sufficiently small to guarantee a building block from the set of exceptional polynomials S, you can use the fractionally faster totalTraceS(d,t) instead- it’ll tell you if you’re wrong!

Calling smallDefectPol(g) will, for a genus g, display the possible polynomials Q (whose roots are the β_i) corresponding to small defect curves (where small means at most 0.780022g, using the exceptional set S). You can get the types (of the form used e.g. by Serre) instead by calling smallDefectTypes(g).

Some defects that don’t meet the bound for small can nonetheless be computed with totalTrace(d,t), use its output as the argument for zetaTypes(X) to recover their corresponding types. This will become more useful if I add additional cases from the data tables.

As I continue to try and get a grasp of heights on elliptic curves, I’ve been looking into the various techniques available for computing them. Since the literal definition as a limit is unworkable, it is common to follow the approach of Silverman and consider the archimedean and nonarchimedean heights in turn. This isn’t too difficult to implement (it’s given in Cremona’s book as well as Silverman’s paper¹), but requires a switching process to avoid certain technical wrinkles. An unpublished paper by my supervisor gives a variation on the theme which avoids the switching; I managed to get both working - and agreeing - in Maple fairly easily.

However, there is a radically different approach possible using Elliptic Divisibility Sequences. These a priori needn’t have anything to do with elliptic curves: they’re simply a type of sequence governed by a not-too-unpleasant rule for iteration: plenty of information on them is available from Graham Everest’s homepage. What’s remarkable, however, is that given an elliptic curve a sequence can be defined which allows you to recover the height of points on that curve. There are numerous advantages to this: it doesn’t matter whether you’re working over the rationals or some other number field; you needn’t have a minimal equation for the curve; and computation is essentially mechanical, as the algebraic geometry is hidden from view. For instance, to follow Silverman’s approach (or variations) it is necessary to factorise the discrinant of the curve to establish primes of bad reduction- this may take up the lions’ share of computation time, but isn’t needed when working via EDS.

There are however sacrifices to be made regarding efficiency or accuracy. Searching for points of very small height clearly only requires a few decimal places of accuracy to indicate whether you’re on the right track; but working purely from the definition of an EDS even this can be cripplingly slow to obtain. Thus, guided by Rachel Shipsey’s PhD thesis², I’ve been trying to create a useful implementation in Maple. The fruits of my labour can be obtained here; they’re still not perfect, as a couple of the examples given by Everest and Ward still cause things to grind to a halt.

Nonetheless, they chew through the others quite happily, and are a lot simpler to use than they were to write! To get started, simply run eds_init with three arguments- the second, third and fourth terms of your sequence. Terms 0 and 1 are assumed to be 0 and 1, for a ‘proper’ EDS. Then, a call to termk(n) will give you the nth term of the sequence. You’ll find that globally there’s a vector h storing all known terms, and another h_known indicating just which those are: typically (and preferably) not all of them, since with Shipsey’s techniques determining the six neighbours of a term h_n allows for a jump to the terms h_2n (or h_2n+1) and its neighbours. In this way, a double and add approach can be employed to obtain a desired term efficiently. There’s still room for improvement here, since this is only the case k=1 from Shipsey’s algorithm: knowing larger factors k allows yet greater jumps through the sequence (to h_2kn from h_k and h_n etc). I’m also tempted to try coding them in PARI instead; I think Maple is stumbling more on performing calculations such as finding norms and gcds than determining the terms themselves, and PARI may be better suited to these.

To get an estimate of the height (as described by Everest and Ward³) , use height_everest(n,d,acc,x,y): n is term of the sequence you wish to work with (the higher the better but slower); d should be the degree of the field etension over the rationals; acc is the level of accuracy (in decimal places) passed to Maple’s eval function during the approximation of the logarithm; whilst x and y are the coordinates of the point of interest.

References

1: Silverman- Computing heights on elliptic curves. (MathSciNet entry)

2: Shipsey- Elliptic Divisibility Sequences (Download)

3: Everest, Ward- The canonical height of an algebraic point on an elliptic curve. (MathSciNet entry) (Download [PS])

One of the central results in the study of Elliptic curves is the Mordell-Weil theorem, which asserts that the group E(K) is finitely generated. Thus it consists of a finite part- the torsion subgroup - and a free abelian part, the rank of which is notoriously difficult to compute. However, the torsion subgroup is relatively accessible, and this is something I’ve been playing with for a while. It covers a range of techniques and ideas and attempting a concrete implementation in Maple has helped considerably in my understanding of those, even if it is effectively reinventing the wheel given the existence of John Cremona’s Algorithms for elliptic curves. The procedures themselves and worked examples are after the cut; first, some theory.

Mazur’s Theorem

Let E/Q be an elliptic curve. Then the torsion subgroup E_tors(Q) is one of the following fifteen groups:

Z/nZ for 1≤n≤10 or n=12;

Z/2Z X Z/2nZ 1≤n≤4.

Further, each of these groups does occur as an E_tors(Q).

This result is particularly handy as it allows for an experimental approach to be taken, gathering enough computational evidence to determine which form the torsion subgroup takes; knowledge of the order of points being especially useful. For instance, the presence of an order 7 element instantly shows that E_tors(Q) is Z/7Z. Better still, there are results which aid in finding such points:

Nagell-Lutz Theorem

Let E/Q be an elliptic curve of the form

(that is, with the usual labelling of coefficients, a_1=a_3=0) with a,b,c integers. If P an element of E(Q) has finite order then x(P), y(P) are also integers.

Further, For such a point either y(P)=0 or y(P) divides

Hence for such curves it is sufficient to look for integer points; and only finitely many such points are suitable candidates for being torsion points.

Good and Bad Reduction

What of Elliptic curves not in the above form? It is possible to bound the number of torsion elements (and generate candidates) by working over finite fields (which I’ve coincidentally considered before). Save for finitely many primes of bad reduction - those which divide the discriminant of the elliptic curve - it transpires that the torsion subgroup maps injectively to E(F_p). For small primes, this is readily found without anything more sophisticated than brute force. Testing a number of primes can give an upper bound whilst naive searches for integer points can provide a lower bound: appealing again to Mazur’s theorem then usually settles the question.

Read the rest of this entry »

I’ve added a couple of extra toys to my Maple procedures for elliptic curves. The major change is that it now supports calculation over some finite fields; that is, the integers modulo some prime. To activate this, set workModM to true and specify a modulus M. Then the usual commands ella, ellm, ncopies and mnadd will compute answers mod M instead.

This also makes it much more likely that you’ll be interested in the order of a point, so a procedure modgetorder is included to calculate this by brute force- that is, repeated addition until the zero element is reached.

This makes questions of the type I faced in MA40188: Algebraic Curves much easier. For instance, consider the curve

Over the field with 37 elements, and with a suitable dehomogenisation, the point P: (x,y)=(0,23) is easily verified as an element of E. Then we may easily determine the point Q=-2P, the third intersection of E with the tangent to E at P:

>read “gla.mpl”;

> a_1:=0;a_2:=0;a_4:=-9;a_3:=0;a_6:=11;

>workModM:=true;

>

>M:=37;

>Q:=ncopies(-2,0,23);

1,22

So Q=(1,22). Further, Q is an inflexion point: that is, the tangent to E at Q meets E three times at Q. In terms of the group law, this means -2Q=Q, or equivalently 3Q=0. We can verify this in a couple of ways:

> ncopies(3,Q);

zero

> modgetorder(Q);

3

Since Q=-2P and 3Q=0, it should follow that 6P=0. Which, fortunately, it does:

> modgetorder(0,23);

6

Here are some applications of the procedures developed in the previous post.

Example 1

We consider example 2.4/problem 3.4 from Silverman;

By inspection we can identify some integer points, such as P₁=(-2,3) and P₃=(2,5). A brute force search for x in the range -1000 to 1000 generates the following results-

> #naive point search;

> for k from -1000 to 1000 do

> if(type(simplify((k^3+17)^(1/2)),integer)) then print(k,simplify((k^3+17)^(1/2))); end;

> end;

-2, 3

-1, 4

2, 5

4, 9

8, 23

43, 282

52, 375

Silverman tells us there is a further point, P₈=(5234,378661), plus we have missed all the inverses of our points (since only the positive square root was computed). Brute force on the range -6000 to 6000 of course uncovers P₈, but this computation takes 70.6 seconds, 10.24mb of memory and produces alarming sounds from my new office computer. Silverman observes that (due to a result of Nagell) the rational points are generated by integer combinations of P₁, P₃, so we can proceed by testing some of these instead:

> read “gla.mpl”;

> a_1:=0;a_3:=0;a_2:=0;a_4:=0;a_6:=17; #setting
up the curve;

>#smarter search

> for i from -5 to 5 do

> for j from -5 to 5 do

> if(type(mnadd(i,j,-2,3,-1,4)[1],integer) and type(mnadd(i,j,-2,3,-1,4)[2],integer)) then print(i,j,mnadd(i,j,-2,3,-1,4));

> end if;

> end:

> end:

-3, -2, 43, -282

-2, -3, 5234, -378661

-2, -1, 2, -5

-2, 0, 8, 23

-1, -1, 4, 9

-1, 0, -2, -3

-1, 1, 52, -375

0, -1, -1, -4

0, 1, -1, 4

1, -1, 52, 375

1, 0, -2, 3

1, 1, 4, -9

2, 0, 8, -23

2, 1, 2, 5

2, 3, 5234, 378661

3, 2, 43, 282

All sixteen points are recovered in 0.02 seconds, consuming merely 0.31mb of memory!

Of course, for this example I’m cheating somewhat because I know where I’d like to get to in that I know this list is complete; although a priori there’s no indication of how large the arguments m,n needed to be to generate points such as P₇ or P₈. Nonetheless, this indicates that the procedures allow for more rapid exploration of points on the curve, even if they don’t prove anything (besides existence) by themselves.

Example 2

Maple’s own algcurves package can also be useful to tackle problems given in projective terms. For instance, we can rapidly demonstrate the first result claimed in Exercise 3.3b. Here we are concerned with the curve

which homogenizing away from Z=0 gives

However, this is not of Weierstrass form; but we can retrieve this from Maple:

> with(algcurves):

> f:=x^3+y^3-A

> Weierstrassform(f,x,y,x0,y0)

This yields

But this is not quite of the Weierstrass form as used in Silverman; we substitute -x₀ for x₀ to arrive at

That is, we have that the curve coefficients a_i are all zero except a₆=27A²/4; we also have an isomorphism φ between E and its Weierstrass form given coordinatewise. We can verify with the procedure j_invariant that these are indeed the same curve (it turns out to have j invariant zero, too). Moreover, we can show the desired result, that

For this, let P=(x_0,y_0) a point on the curve in Weierstrass form. Then we compute -P:

> a_1:=0;a_3:=0;a_2:=0;a_4:=0,a_6=-27*A^2/4:

> read “gla.mpl”:

> ellm(x_0,y_0);

x_0, -y_0

Then, identifying P with a projective point via the isomorphism, we find

Which is the desired result.