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.

Last time I introduced the idea of divisors on a curve; but an observant reader may have noticed that along the way the idea of rational points seemed to be lost. Further, whilst the Riemann-Roch theorem guarantees that a divisor from the jacobian will have a reduced representative, no indication was given as to how that representative is to be found. In this post I’ll try to clear up both of these issues.

Recall that a semi-reduced divisor D from the jacobian takes the form (Σ_i^rP_i) - r∞ where the P_i are points (x_i,y_i) of C. We will represent this by a pair of polynomials D=div(a(u),b(u)) with the following properties:

such that b has degree less than that of a, and the appropriate multiplicity for repeated points- i.e., if P_i occurs k times in the semi-reduced representation of D, then (u-x_i)^k divides b-y_i. This ensures uniqueness.

The empty divisor (zero element of the jacobian) is denoted by div(1,0); if a is linear (and hence b constant) then the divisor corresponds to a single point of C; for the point (x,y) the divisor is div(u-x,y). The degree of a is described as the weight; ‘most’ reduced divisors will be of weight g. Recall that the co-ordinates of a point were only required to be in A rather than K; we describe a divisor as rational over K if the coefficients of a and b are from K. Beware that a rational divisor may therefore be the sum of points which are not K-rational points of the curve; however, a weight 1 rational divisor obviously corresponds to a K-rational point.

The K-rational principal divisors of C are a subgroup of P₀ and their image in J, J_K, the subgroup of rational divisors of the Jacobian, is the object of computational interest. The K-rational points of C are then identified with a subset of J_K; namely the divisors of weight at most 1; except in genus 1 (where they are J_K), this isn’t a subgroup due to the lack of closure.

So we wish to work with rational divisors from J_K.Given such divisors of the form div(a,b), it is undesirable to construct their sum by ‘unpacking’ the points P_i and forming new polynomials as the co-ordinates (roots of a, evaluations of b) might not be from K but A. Fortunately, it is also unnecessary: the only complications are connected to repeated points or the combination of a point and its negative; careful manipulation of gcds allows for direct computation of the semi-reduced form. See ¹ for details, which also describes moving from semi-reduced to reduced form. To do this, note that for a divisor D=div(a,b), the divisor E=-((b-v)-D)=div(a’,b’) is equivalent to D with deg(a’)=max(2g+1,2)-deg(a); thus by repeated iteration we can move to a reduced representative (the explicit formulae are a’=(f-b²)/a, b’=-b mod a’).

Maple procedures to do all this will be provided in the next post.

Reference
¹Computing in the Jacobian of a Hyperelliptic Curve D.Gantor Mathematics of Comptuation Vol.48.No.177 (Jan., 1987).

One of the most celebrated properties of elliptic curves is that the set of rational points is a group, with a highly geometric explanation of the group law: the ‘chord and tangent’ process. Two points and their sum are linked by consideration of the intersection of straight lines with the curve: as the curve is a cubic, there are three intersections (subject to some technical book-keeping with repeated points and the point ‘at infinity’). Such an approach clearly won’t transfer immediately to curves given by higher degree polynomials, as there will be more intersections, but as there are only finitely many, one would still hope to be able to define relationships between points. For instance, on an elliptic curve, if A, B and ∞ are colinear then B is -A; thus if on a more complicated curve we had A,B,C,∞ colinear it might make sense to think of C as -(A+B), and then A+C as -B and so on. That is to say, there may be relationships between groups of points rather than individual points.

In algebraic geometry, this (and much more) is captured by the notion of a divisor; rather than present them here in full generality, I will consider the specific case of divisors on (hyperelliptic) curves. These will then serve as the building blocks for a group structure connected to the curve which reduces in the special case of an elliptic curve to the familiar group of rational points.

To fix ideas, let K be a field of characteristic other than 2 with algebraic closure A. A curve C is described as a hyperelliptic curve of genus g if there is some degree 2g+1 polynomial f with distinct roots such that v²=f(u) is a model of C: so the familiar elliptic curves are the special cases with genus 1.

A point P on C is a pair (x,y) of elements of A (not K) satisfying y=f(x); or the point at infinity ∞. Then a divisor D of C is a finite formal sum Σ_im_iP_i for integers m_i and points P_i on C. D is described as having degree Σ_im_i; if all the m_i≥0 then we write D≥0. Formal (that is, pointwise) addition of divisors gives the additive group D of divisors; its identity is the divisor consisting of summing no points and it has a subgroup D₀ consisting of divisors with degree zero.

Any polynomial p(u,v) can be considered as a function on C of the form p=a(u)+b(u)v, since v²=f(u). If p vanishes at (x,y) then the order of the zero (x,y) of p is the exponent of the highest power of (u-x) which divides a²-b²f.

Thus we can define functions on C as h=p/q for p,q polynomials from K[u,v] such that v²-f(u) does not divide q(u,v): that is, q is not everywhere zero on C. Then h will have a finite set of zeros (those of p) and of poles (those of q); we associate to h a divisor, (h) = Σ_im_iP_i where the P_i are those zeros and poles and m_i their multiplicities:

If there is a nonzero function h on C such that D a divisor is (h), then D is described as principal. the principal divisors form a subgroup P of D₀ and hence D: the jacobian J of C is then the quotient D₀/P. That is, two divisors correspond to the same element of the jacobian if they differ by a principal divisor. This gives some idea as to how to simplify arbitrary divisors- we work in the jacobian and seek a simplest representative; that is, one comprised of the minimal number of points.

Consider that if (x,y) is a point P of C, then so is P’=(x,-y). The function u-x has zeros P and P’ with a double pole at ∞ so P+P’+∞ = (u-x) is principal and hence equivalent to zero mod P. Hence -P’ is equivalent to P-2∞ so we can rewrite divisors to only feature positive multiples of points other than ∞ Thus in J, where the degree is necessarily 0, any element has a representation

such that if P_i appears in D, then no P_j=P’ for any j different to i. Hence, any point of the form (x,0) will appear at most once. Such a representation is called semi-reduced; if r≤g then it is called reduced.

Remarkably, (by the Riemann-Roch theorem) any divisor in the Jacobian will have a unique, reduced representative (in other words, any divisor is the sum of a reduced divisor and a principal divisor). Now we can see what’s really going on with the elliptic curve group law: as a reduced divisor will have r≤1, it takes the form P-&infty; so there is an obvious isomorphism between the set of rational points and the Jacobian. Hence adding two points A,B on the curve gives rise to another point of the curve, by reducing the divisor A+B-2∞ to some representative C-∞ and setting A+B=C.

But with hyperelliptic curves, this needn’t be the case: the sum of two points is a perfectly good reduced divisor in the next simplest case of genus 2, for instance, so we can’t add two points and expect the answer to be a point. Hence we need to consider the divisors corresponding to rational points in the broader setting of the jacobian; to extract useful information about those points, we’ll need to consider the rational divisors. This motivates an alternative notation for divisors, more suitable to computation: I leave all these issues to the next post.

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] !

The discrete logarithm problem is the vital part of elliptic curve cryptography, but can be defined (to varying cryptographic strength) for any cyclic group:

Discrete Logarithm Problem (DLP)
Let G be a cyclic group, with operation ⊕. Let [n] represent (n-1)-fold application of ⊕ i.e., [2]g=g⊕g, [3]g=g⊕g⊕g etc.

Given g,h from G, the discrete logarithm problem is to find t such that [t]g=h.

The baby-step, giant step (BSGS) algorithm is a generic algorithm for solving the DLP- that is to say, it makes no appeal to properties of the group involved, merely calculating abstractly with ⊕. There is a theoretical upper bound to the effectiveness of generic algorithms, and BSGS approaches are of that order of magnitude.

The simplest demonstration of BSGS (the original by Shanks) assumes that g generates G, with both of known order n: recall that the group order is simply the number of elements it contains, whereas the order of an element g is the least n such that [n]g=id_G, should such a value exist. The order n of an element always divides the order of the group N, with equality when g generates G. Between the two there is also the exponent, the least value e such that [e]g=id_G for all g in G; for cyclic groups this is N, but for products of cyclic groups (the structure of groups of rational points) it is the lowest common multiple of their respective orders, and thus, although a divisor of N, may be significantly smaller than it.

So, if we can find an order m rational point on a curve, we know that the cardinality of the group of rational points is zero modulo m. By testing random points and taking the lowest common multiple of their order we can usually find the exponent e of the group. With luck, but not always, this will be large enough that when combined with bounds on the cardinality the latter is established exactly (as with modular information in the SEA algorithm).

But it is no good attempting to do so using an algorithm that requires the order of the point! There are BSGS algorithms for the DLP which can handle unknown order, for a given g we can apply these to solve for h=id_G in the cyclic group generated by g. Provided the algorithm gives the minimal t such that [t]g=h=id_G in G=<g>, then t is the order of g.

One such algorithm, by Terr, is given in Cohen and Frey’s Handbook of Elliptic and Hyperelliptic curve cryptography (my current bible!) It relies on the following observation:

Let T_n be the nth triangular number: that is, defined by the recursion T₁=0, T_n+1=T_n+n.

Then any non-negative integer t, there are unique integers j and k with t=T_j+1-k with 0≤k<j.

To see how, notice that there is a unique j such that t lies in the interval (T_j , … , T_j+1]=(T_j+1-j , … , T_j+1-0] and that this interval has width j so there is an appropriate choice of k.

We then suppose that t satisfies [t]g=h. Then we have [T_j+1]g=h⊕[k]g. So, instead of testing each [t] in turn until we hit equality, for a given j we need only test the ‘giant step’ [T_j+1]g against the set of ‘baby steps’ β={β_i}={h⊕[0]g,…,h⊕[i]g,…,h⊕[j-1]g. If that j yields no matches, we move to j’=j+1: by keeping track of [j]g at each iteration, the new T_j’+1 and the extra h⊕[j’-1]g are found by a single group operation each. Hence this is attractive from both storage and time complexity perspectives: we need only record the baby steps β, [j]g and a single giant step at any given iteration j, whilst the time complexity is of order t^1/2.

Terr’s BSGS variant for the DLP
Finds t such that [t]g=h for g generating G of unknown order, h from G.

Initialise β={β₀=h}, γ=δ=id_G, j=0. ((For each iteration, we have γ=[T_j+1]g and δ=[j]g)
Then loop over j as follows:

• Increment j by 1, then update:
• Set δ=δ⊕g=[j-1]g⊕g=[j]g.
• Set γ=γ⊕δ=[T_j]g+[j]g=[T_j+1]g.
• If j≥2 then
  • For s from 0 to j, if γ=β_s then return T_j+1-s
• Add β_j=β_j-1⊕g to β

I’ve implemented this generic DLP algorithm and an order-finding version (which requires you to catch g=id_G) as Maple procedures for arbitrary groups- details on that will be in the next post! BSGS-based approaches become impractical at finite field sizes well within the grasp of SEA for counting on elliptic curves; but are of interest in the higher genus case which lacks an Elkies procedure.

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.

View as:

Since the restart of term I’ve been studying rational points of curves over finite fields, trying to see how the geometry and number theory can mesh together to say more than either approach offers in isolation. The Weil Bounds give an upper limit on the number of points in terms of the size of the field and the genus of the curve. There is practical interest in finding curves with a large number of points to get good error-correcting codes; the shortfall between the bound and the actual number of points is described as the defect of the curve.

The ‘types’ of zeta function can be characterised in terms of their roots; number-theoretic arguments demonstrate that only finitely many types give curves of small defect. Geometrically this ‘type’ specifys the eigenvalues of the Frobenius endomorphism on the Jacobian of the curve; and thus under certain conditions some more types can also be eliminated, which may render particular defect levels impossible. Conversely, eliminating small defects by better number-theory bounds means the corresponding types cannot arise from the geometry.

This summary article describes the Weil bounds; some results on algebraic integers and their application to determining the possible types of zeta function (following Serre, with a slight refinement); a worked example for genus 3 plus the types for defect 0,1 and 2; and a brief sketch of some of the geometric arguments/results that arise.

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