Modulo Errors

So it’s been over two months since a post; more attentive readers will have noticed that there was one, but now there isn’t. I’ve moved away from thinking about cryptography to generalising some number/graph theoretic results of my supervisor, concerning matrices with constrained eigenvalues. However, this creates a problem: unless I ‘blog every up and down of the research process (which could be interesting, but would slow me down!) information on here becomes decreasingly accurate or relevant as I revise my thinking on the topic. Certainly it would be premature to present firm results at the moment.

But I can at least set the stage for more technically-minded readers (a friendlier explanation/illustration will hopefully follow once I truly understand all this!). Chris has characterised all symmetric integer matrices with the property that their eigenvalues are at most 2 in modulus; under a suitable transformation of their characteristic polynomials, these give cyclotomic polynomials and thus are referred to as cyclotomic matrices. Conveniently, any submatrix of a cyclotomic matrix is itself cyclotomic, so it suffices to find maximal examples. Although there are infinite families of these matrices, there are only a few ‘types’ possible.

These types are best understood by considering not the matrix, but an associated graph, where values in the matrix determine the weights on edges and nodes of the graph. This introduces a notion of equivalence, since many matrices will correspond to the same graph or certain well-defined variations on it. Further, we can adjoin nodes and edges to the corresponding graph to try and ‘grow’ towards maximal examples.

The motivation comes from finding polynomials of small Mahler measure- whilst a cyclotomic polynomial has measure 1, all others seem to be pushed away, with the smallest known value being 1.176… The question is how to generate small examples, and these matrices provide a way: by adjoining a single extra node to a maximal cyclotomic graph, a non-cyclotomic graph/matrix is obtained and thus a non-cyclotomic polynomial. The minimal graphs with this property (non-cyclotomic, but all subgraphs cyclotomic) often correspond to polynomials with some of the smallest known Mahler measures.

But some examples are not generated in this way, which is where I’ve stepped in. There is no reason to restrict attention to integer matrices, and I’ve established which imaginary quadratic extensions of the rationals give rise to rings of integers over which suitable matrices can be found. For a couple of fields, there are very few new (non rational-integer) cyclotomic matrices, and I have a complete description of them, but in others there are again infinite families as well as occasional examples that don’t generalise.

So I explore this behaviour by growing graphs/matrices, and try to spot patterns as they emerge from the fragments. I use the university’s parallel computing cluster Eddie for brute force work in SAGE, but such is the nature of the combinatorial explosion that even this doesn’t suffice without some mathematical insight along the way, as I try to refine my growing algorithms and capture equivalence as early as possible. I’m hopefully nearing the point where all examples fit into known families, at which point I’ll need to switch into serious mathematician mode and try and prove why this should be so. But for now I need to make sure that nothing unexpected tumbles out of each batch of calculations!

On a completely unrelated note, I’ve dragged modulo errors up to date with wordpress 2.5 and switched themes; please shout if you find I’ve broken something along the way.

The latest version of my ellnet class is ellnet2d_lowmem.spyx. It combines all the tricks I know of:

• The use of precomputed inverses for all steps, and precomputed squares/products for each step, as described by Stange,
• computation with a local vector to avoid overhead from function calls to keep the dictionary up-to-date,
• mixed block lengths as described in the previous post,
• and compilation to pyrex.

Thus it’s the fastest implementation I currently have for finding Tate pairings in SAGE (about twice as fast as accessing Stange’s PARI script from SAGE). Attach it in the normal way; example calculations are here.

A couple of comments now that I’ve tested my Tate pairing procedures on some random input.

First, it turns out to be not too difficult to make use of PARI scripts from within SAGE. I spent ages trying to figure out how to call them through the gp interface, until I discovered that new functions simply become additional methods of the gp object; although some SAGE objects have to be converted to their PARI equivalent before you can use them as input. For instance, assuming you have Stange’s tate_via_nets script in your working directory, you can attach it with gp.read(”tate_via_nets.gp”). Then for SAGE points P,Q on a curve E with the order of P stored as n, the following will compute the Tate pairing via PARI:

Egp=gp(E)
Pgp=gp([gp(P.xy()[0]),gp(P.xy()[1])])
Qgp=gp([gp(Q.xy()[0]),gp(Q.xy()[1])])
gp.tate_pairing_alg(Egp,Pgp,Qgp,n)

It’s hard to get accurate measures of the time taken by a PARI subprocess (it’s invisible to SAGE’s cputime, and walltime will be influenced by any other running processes), and these calculations are pretty quick, but it seemed the PARI implementation was faster. So I’ve developed a SAGE version which avoids the overhead of dictionary read/writes by not caching intermediate terms; I’ve also set it up as a pyrex file, which despite containing purely python code (I’ve not ported any of the internals to C) is faster than both the old ellnet2d and the above PARI approach after a one-time compilation. You can download that here; then use attach “ellnet2d_lowmem.spyx” instead of attach “ellnet2d.sage” before proceeding as before.

The good news is that in every test I’ve run (over both prime and prime-power fields), the various implementations all agree. The bad news is that it doesn’t seem possible to compute elliptic curve point orders beyond about 64bits in SAGE (or rather in PARI, upon which it depends), and before that stage the order calculation is significantly slower than the pairing computation. It also seems that a random curve tends not to feature particularly high powers of 2 in its order, so the speed-up I suggested earlier doesn’t offer much advantage at all. Hopefully I can find a class of curves where it is of merit!

Eventually, then, we’re in a position to compute Tate pairings in SAGE. The necessary ingredients are an elliptic net class (the more general one or a faster version for pairings only), an eds class (this one based on Stange’s algorithm is recommended and used in the examples below; the old one should also work), and supporting procedures tate.

Load these in the usual fashion:

sage: attach "tate.sage"
sage: attach "ellnet1d.sage"
sage: attach "ellnet2d.sage"

We now need to introduce a curve and some points. Stange gives a simple example of y^2+y=x^3+x^2-2x over the field with 5 elements corresponding to points P=(0,0) and Q=(1,0). We enter the field, curve and points and set up an elliptic net:

sage: K=GF(5)
sage: E=EllipticCurve(K,[0,1,1,-2,0])
sage: E
Elliptic Curve defined by y^2 + y = x^3 + x^2 + 3*x over Finite Field of size 5
sage: P=E([0,0])
sage: Q=E([1,0])
sage: X=netFromPoints(P,Q)
sage: X

Elliptic net with (1,0) block:
                1       1       1
4       4       0       1       1       2       1       3

The block centred at 1 is displayed; this has been generated from the initial data, and suffices to compute any term of the form (x,0) or (x,1). If such a term is known, it is simply retrieved from the hash table; otherwise, the algorithm is employed. Requests for unreachable points give a ValueError:

sage: X.has_element((5,0))
True
sage: X[(5,0)]
3
sage: X.has_element((100,1))
False
sage: time X[(100,1)]
Computing terms
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.01
1
sage: X.has_element((100,1))
True
sage: time X[(100,1)]
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 0.00
1
sage: X[(100,2)]
---------------------------------------------------------------------------
<type 'exceptions.ValueError'>            Traceback (most recent call last)

/home/s0677951/SAGE/eds/<ipython console> in <module>()

/home/s0677951/SAGE/eds/ellnet2d.py in __getitem__(self, A)
     41                                 # Term not entered and cannot be computed; give an error.
     42                                 raise ValueError, \
---> 43                                 "Term unavailable, cannot compute unless of form (x,0) or (x,1)"
     44
     45         def has_element(self, A):

<type 'exceptions.ValueError'>: Term unavailable, cannot compute unless of form (x,0) or (x,1)

You can also test for and print entire (small)blocks; this illustrates the route taken during computation of elements, for instance in computing the 18th term, the blocks with centres 1, 2, 4, 9 are constructed, then the small block centred on 18:

sage: X.has_block((9,0))
False
sage: X[(18,0)]
Computing terms
0
sage: X.has_block((9,0))
True
sage: X.print_block(9,0)
                1       1       3
4       3       2       0       3       2       1       2
sage: X.has_block((18,0))
False
sage: X.has_small_block((18,0))
True
sage: X.print_small_block(18,0)
                2       4       3
3       4       4       0       1       1       2

To compute the pairing of P and Q, use TatePairing(m,P,Q), where m is the order of P (determining this left as an exercise!):

sage: TatePairing(9,P,Q)
Computing terms
3

Finally, a call to TatePairing with P=Q will construct an eds rather than an ellnet. We can verify an example of Stange’s (see slides from ECC):

sage: K=GF(73)
sage: E=EllipticCurve(K,[0,1,1,-2,0])
sage: P=E([0,0])
sage: TatePairing(9,P,P)
Computing terms
24

Since EDS are the rank 1 case of elliptic nets, Stange’s approach to the rank 2 case can be specialised to them: we are only concerned with terms of the form (x,0), so only the bottom row of a block need be calculated. Such sequences are considered when pairing a point with itself.

Thus I’ve simplified the ellnet class to give another eds class. I’ve changed the internals so that this class can be initialized and interacted with in the same way as my previous eds class based on Shipsey’s algorithm. (You retrieve the n term from a rank 1 net X using X[n] rather than the X[(n,0)] required for a general elliptic net, for instance.) There will probably be a few discrepencies to be ironed out, but I’ve thrown in the ward_point methods necessary for height calculation, and hope to have the two versions interchangeable in due course.

This requires only 1 initial inversion, and continues to use the step-by-step precomputation speedups given by Stange. It doesn’t support the larger jumps possible with Shipsey’s algorithm, although my existing code doesn’t exploit these as fully as it could anyway. Final doublings are still performed with smaller blocks, so for height estimation purposes it’s probably wise to choose a power of two for the term used.

The next post illustrates pairing computation with rank 1 nets.

I’ve made a few cosmetic and minor functional tweaks to some of my SAGE offerings. The latest version of eds suppresses some diagnostic information, and the representation of an eds object now gives the defining terms rather than the entire dictionary (which runs to scores of pages for sequences over the integers). Further, heights will now catch the zero division exception that sometimes arises when you use Shipsey’s algorithm on a sequence corresponding to a point of finite order. It will now give a warning to this effect, and return a height of zero as required.
This, coupled with the singular curve catch introduced last time should mean that batch testing of large numbers of curves shouldn’t now fail- but I’m sure I’ll find some new exception in due course!

I’ve taken stock of my various coding efforts listed on my university space and trimmed out the least useful / developed / efficient ones from the past year. What’s left generally offers features not otherwise available for the given CAS- generic group operations and their use for [hyper]elliptic curve groups in Maple; EDS/height calculations for SAGE. The bug I found in hyperelliptic curve divisor class addition over finite fields has since been fixed, making my alternative class redundant; and the NTL interface seems to have changed, rendering my faster version broken; so those have been dropped from the site.

Finally, I’ll be back in the West Country for the SAGE Days 6 talks this weekend.

Last time, I gave some procedures to generate an elliptic divisibility sequence from a point of an elliptic curve, the motivation being that the growth rate of the sequence corresponds to the height of the point- and thus we can compute heights by computing terms of the sequence.

But what about going the other way? Given that we seek points of low height, can we search more intelligently by crafting sequences with a small growth rate, rather than simply mapping from (typically uninteresting) random points to their sequence? The easiest way to achieve this would be to start from the sequence, rather than from the point- but that requires a means to recover the point from the sequence.

Fortunately, Ward’s Memoir¹ provides just such a technique, although presented in terms of elliptic functions (the complex analysis way to view elliptic curves). Recall that to define a sequence, we need only specify the second, third and fourth terms h₂, h₃ and h₄. Then the corresponding point is (X,Y) on curve E: Y²=4X³ -g₂X -g₃, given by the following unpleasant formulae:

However, these needn’t be integral, and it’s preferable to rescale to eliminate that X³ coefficient. Fortunately both of these problems can be solved by clearing denominators, leading to a more pleasant representation of the point and curve. My basic algebra being as rusty as it is, this procedure me a surprisingly long time to figure out in the general case, but I have now added suitable features to my eds SAGE class. For such an object X, a call to X.ward_point() will return a point of an elliptic curve: you can extract the curve from a point P in SAGE by P.curve(). There’s a slight issue in that an EDS may define a singular curve: I’ve caught this exception, and in such a case you’ll get back the integer zero instead. I needed this to prevent large searches over EDS from crashing when such a sequence arose from the parameters, but it’s a bit of an ugly solution.

Given a method to turn sequences into points, the obvious question is whether it acts as the inverse to the creation of sequences from points. The answer turns out to be, almost! Extra scaling factors creep in to the sequence terms each time you loop through a cycle of EDS→point→EDS. The larger the terms, the greater the time/space requirements for manipulation of the sequence; this also has a knock-on effect on the approximate height computed. However, this can be fixed for curves over number fields: considering field elements as polynomials, rescaling each of h₂,h₃,h₄ to have content 1 recovers the original EDS. I offer no mathematical justification for this other than that it works! Thus heights now also includes EDSfromPoint_nf to generate EDS by just such a scaling.

Armed with these techniques, it’s now easy to generate curves of moderately small height that would otherwise be presumably difficult to stumble upon: for instance, with u=1+√2, the curve over the number field Q(u) given by y²=x³ + (-432u-864)*x + (-1728u+6480) has a the point (24-12u,108) of height around 0.0107. This unlikely looking curve arises from the much simpler EDS with defining terms -1,u-1,2u-2, found after a search over sequences from K with small terms that took about an hour. Considering how easily it was found, this compares favourably to the lowest known height over Q of about 0.0045. However, in the context of the Elliptic Lehmer problem heights should be scaled by the degree of the field extension, so to set a record a value of around 0.002 is needed for quadratic extensions such as K. So the search continues…

References

1: Morgan Ward- Memoir on Elliptic Divisibility Sequences (MathSciNet entry)

Way back in December I spent some time looking at Elliptic Divisibility Sequences in connection with computing heights of points on elliptic curves. At the cryptography conference I recently attended, one of the talks presented a generalisation of EDS to higher dimensions described as elliptic nets. These too find application to elliptic curves- in the two dimensional case, allowing for computation of pairings. I’m hoping to look into these in more depth, and so they’ll probably crop up in later posts - but first I felt I should tie up some loose ends related to EDS.

At the time I’d been working with Maple, which proved unwieldy for working over number fields such as those discussed in Everest and Ward’s paper². Further, I only got as far as implementing a special case of Shipsey’s algorithm for computing terms of such a sequence efficiently. SAGE’s object-oriented approach provides a vastly superior environment for algebraic geometry to Maple, which helped solve the first problem; and my mathematical understanding has progressed sufficiently to get to grips with the full version of Shipsey’s algorithm (and fill in the blanks necessary for an implementation).

So eds is a SAGE class for arbitrary Elliptic Divisibility Sequences. It’s pretty basic at present but fit for purpose- terms of the sequence are stored in a hash table, and when an as-yet uncomputed term n is requested the factors of n are compared to the known septuples to establish an efficient chain to the desired term from existing data. There’s still room for improvement here, however, as there are seven tuples containing the nth term of varying difficulty to compute, and there’s no attempt to compute some earlier tuples to speed things up further yet. These are both approaches I’d like to attempt before delving into the more general elliptic nets.

These are implementation concerns however, and shouldn’t alter interaction with the class, which I’ve hopefully made as simple as possible. As before, the working assumption is of a ‘proper’ EDS with leading terms 0, 1. Supplying terms 2, 3 and 4 is then sufficient to specify the entire sequence: initially the object is populated with the first eight terms. You also need to supply a ring of definition for the terms of the sequence so that SAGE doesn’t trip up- at least, until I find a way to determine this unambiguously from the input.

Using the class

The simplest example of an EDS is the integers themselves. So we can construct this EDS X with arguments 2,3,4:

sage: attach "eds.sage"
sage: X=eds(ZZ,2,3,4)
sage: X
EDS over Integer Ring with terms {0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6, 7: 7}

Here I’ve attached the class, specified an EDS and asked for it back - at the moment, such a request gives you the ring of definition and all the currently-stored terms. Individual terms are extracted in the obvious way:

sage: X[7]
7
sage: X[25]
Calculating via Shipsey
Already have septuple with centre 1 , binary-chaining from there.
25
sage: X[32]
Calculating via Shipsey
Already have septuple with centre 8 , binary-chaining from there.
32

If a term is already known (such as the 7 here), it’s simply returned from the hash table. Otherwise, Shipsey’s algorithm (Thm 3.4.1 of her thesis¹) is employed, after determination of a factor k to use as a starting block: for 25, as the tuple about 5 wasn’t known, the basic double-and-add has been employed; but for 32, the tuple about 8 allows for faster computation.
You can test for the existence of both individual terms and entire septuples within the hash table:

sage: X.has_element(32)
True
sage: X.has_septuple(32)
True
sage: X.has_element(35)
True
sage: X.has_septuple(35)
False
sage: X.has_element(38)
False

In this case, both the 32nd term and it’s surrounding tuple are known; term 35 is known but not the tuple; and term 38 is not known at all.

EDS from elliptic curves- finding heights

Much of the interest in EDS arises from the fact that for a point P of an elliptic curve E, the division polynomials φ_n(P) are just such a sequence. These polynomials have the property that the zeroes of φ_n are the x co-ordinates of the n-torsion points, and feature in the explicit formulae for the group law on an elliptic curve; their computation is also of relevance to the SEA algorithm.

heights contains supplementary functions for constructing EDS with terms corresponding to the φ_n(P); technically you can specify an EDS with terms the φ_n themselves then evaluate for chosen P, but this is much slower and thus it is generally preferable to attach sequences to points rather than curves. The basic formula (21) from Everest and Ward² then gives a means to estimate the height of the point; SAGE can easily make sense at runtime of the norm, degree of field extension and so on that are required.

For instance, their example 4 gives the curve y²+y=x³-x² over K=Q(√ -2) and the point Q=(2+√-2,1+2√-2). We construct the curve and EDS (using EDSfromPoint) and recover a height estimate (based here on the 100th term of the sequence, using heightn) as follows:

sage: K=NumberField(x^2+2,'a')
sage: a=K.0
sage: E=EllipticCurve(K,[0,-1,1,0,0])
sage: E
Elliptic Curve defined by y^2 + y = x^3 + (-1)*x^2
over Number Field in a with defining polynomial x^2 + 2
sage: Q=E([2+a,1+2*a])
sage: attach "heights.sage"
sage: X=EDSfromPoint(Q)
sage: X
EDS over Number Field in a with defining polynomial x^2 + 2 with terms
{0: 0, 1: 1, 2: 4*a + 3, 3: 11*a - 63, 4: -612*a + 689, 5: 17637*a,
6: -5055750*a - 4675301, 7: 2521168375*a - 4167504129}
sage: time heightn(Q,X,100)
Calculating via Shipsey
Already have septuple with centre 4 , binary-chaining from there.
CPU times: user 0.19 s, sys: 0.01 s, total: 0.20 s
Wall time: 0.20
0.457448706299596

This compares favourably with Silverman’s accurate value of 0.45754…

Note that due to the scaling convention adopted, heights computed in this way will be half those computed by SAGE via PARI. However, PARI can only be used for curves defined over the rationals - with these procedures, it’s possible (to a low degree of accuracy) to determine heights for curves defined over number fields.

References

1: Shipsey- Elliptic Divisibility Sequences (Download)

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

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.