Modulo Errors

On Friday I gave a geometry club seminar, speaking about some of the computational aspects of discrete-logarithm cryptography in general and as implemented for elliptic curves. My notes supplement rather than completely describe the talk, being heavier on the formalities and lighter on the narrative.

The topics covered are: Diffie-Hellman and one-way functions for key exchange; the generic Discrete Logarithm Problem and BSGS algorithm; scalar multiplication- addition chains, fast exponentiation, m-ary methods and windowing; group law implementations, Side-channel attacks and the Edwards form.

I’ve discussed several of these ideas elsewhere on this blog, as well as the cryptanalysis ideas I mentioned on the day but which are not in the notes. I also refered to a recent real-world example of a side-channel attack; see this story from The Register for details.

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

In Shipsey’s algorithm¹ for EDS, the septuples centred on terms mk and (m+1)k can be used to construct the septuple about 2mk, (2m+1)k or (2m+2)k from that about k. With k=1 this gives rise to a double-and-add algorithm, such that term n can be obtained in logarithmic time; knowledge of larger k introduces an additional speedup.

However, Shipsey’s approach makes use of various inverses, potentially giving division-by-zero errors for a sequence corresponding to a point of finite order (any point, in the finite field cases of interest for pairing computations). Further, the recurrence relation for a rank 2 elliptic net is more complicated than for an EDS, and we ultimately wish to find terms (m,0) and (m,1). Stange provides an algorithm to achieve this, replacing the septuple with an 11-element block centred on k as follows:

Notice that the blocks centred at any of m, m-1 or m+1 will contain the information needed to compute the pairing. Thus we generate such a block by a double-and-add algorithm.

Algorithm

Double-and-add Algorithm

INPUT: Integer n and block centred at 1.

OUTPUT: Block centred at n.

1. Compute binary digits d_i of n such that n=(d_kd_k-1…d₁)₂ with d_k=1.
2. Set c=1 (centre)
3. For i=k-1 down to 1 do:
   • If d_i=0 Compute block with centre 2c; Set c to 2c
   • Else, Compute block with centre 2c+1; Set c to 2c+1
4. Return final block.

Clearly, this requires procedures to generate new blocks from old- Stange² gives formulae for these ((12)-(17)). A speed-up (Id. S5.1) can be achieved by computing commonly used squares/products at each step, denoted by the A(i) (squares) and B(i) (products) below.

Double

DoubleAdd

The inversions in the recurrence formulae are independent of k, so are precomputed and made available to Double and DoubleAdd as multipliers. Including multiplication by these values, the total cost of a block double/double-and-add is 35 multiplications, 7 squarings, 11 additions and no inversions.

A minor refinement

Define a 10-element short block centred at k by

shortDouble

Notice from the Double dependencies that the (k+4,0) term is necessary only to generate (2k+4,0). Thus given a short block centred at k, we can double it to obtain the short block centred at 2k, with the following dependencies:

The inversions in the recurrence formulae are independent of k, so are precomputed and made available to shortDouble as multipliers. Including multiplication by these values, the total cost of a short block double is 31 multiplications, 6 squarings, 10 additions and no inversions.

shortDoubleAdd

From the DoubleAdd dependencies it can be seen that the missing (k+4,0) is required for both (2k+4,0) and (2k+5,0)- although the latter can be dropped for a short block, the former cannot. Thus it is not possible to adapt this choice of recurrence formulae to give a DoubleAdd for short blocks. Nonetheless, the basic algorithm can be ammended to exploit the simpler short block Double once no further DoubleAdds are required; as the short block is sufficient for finding the Tate Pairing.

Algorithm

Amended Double-and-add Algorithm

INPUT: Integer n and block centred at 1.

OUTPUT: Block centred at n.

1. Compute s,p such that p is odd and n=2^sp.
2. Compute binary digits d_i of p such that p=(d_kd_k-1…d₁)₂ with d_k=1.
3. Set c=1 (centre)
4. For i=k-1 down to 1 do:
   • If d_i=0 Compute block with centre 2c using Double; Set c to 2c
   • Else, Compute block with centre 2c+1 using DoubleAdd; Set c to 2c+1
5. This gives the block with centre p, and thus the short block with centre p.
6. For i from 1 to s do:
   • Compute short block with centre 2c using shortDouble; Set c to 2c
7. Return final block.

The efficiency gain of this variant depends on n- if n is a power of 2, then we may proceed entirely by shortDouble with 47 instead of 53 operations, a saving of around 11 percent (this will vary depending on the precise cost of multiplication, squaring and addition). If, however, n is odd then there is no gain. Thus pairing computation can be improved in the case of points with order a high power of two.

Implementation

A SAGE ellnet class is available, which implements all these ideas for rank 2 (Stange’s algorithm, with precomputation of inverses, caching of commonly-used factors at each step, and finishing by short blocks). See this post for how to generate nets from points and compute pairings.

References

1: Shipsey- Elliptic Divisibility Sequences (Download)

2: Stange- The Tate Pairing via Elliptic Nets (Preprint)

As always in mathematics, when something turns out to be useful, it’s worth asking whether it can be generalised to a broader setting. For Elliptic Divisibility Sequences, that generalisation is Elliptic Nets. For a free abelian group of finite rank, this is a function W satisfying the condition:

for any p,q,r,s from the group.
(EDS are then recovered as the special case where the chosen group is the integers).

The theory of elliptic nets has been developed extensively by Katherine Stange, and the papers/presentations from her website provide a wealth of information. I’ve been writing a SAGE implementation of her results, from which we need a few key ideas.

What’s remarkable is that it is not just that sequences can be generalised to arrays: the connection to elliptic curves is also preserved. That is, given n distinct points of an elliptic curve, there is a rank n elliptic net corresponding to those points. The terms of the EDS corresponding to each point will appear in the array, but more complicated interactions of the points are also captured. In the rank 1 case, it was possible to estimate heights purely by manipulation of terms of an EDS; in a similar way, routine manipulation of elliptic nets allows for computation of the Tate pairing.

The use of pairings is currently a hot topic in cryptography, with both destructive and constructive applications. The Tate Pairing allows pairs of elements of an elliptic curve group to be mapped to a finite field. In certain cases, and provided the pairing can be computed efficiently, this transplants the Discrete Logarithm Problem (DLP) into an easier environment, thus weakening the cryptosystem. On the other hand, the gap in difficulty between the two DLPs can itself be used to design cryptosystems and protocols, such as Identity based cryptography.

For the theoretical background, I recommend again Stange’s work. The take-home message however, is as follows:

• For a point P of order m on an elliptic curve E, to compute the pairing of P with itself it suffices to recover terms m+1 and m+2 of the EDS associated with P.
• For distinct points P,Q on an elliptic curve E with P of order m, to compute the pairing of P with Q it suffices to recover elements (m,0) and (m,1) of the elliptic net associated with P and Q.

Thus I have been interested in the sequence/net calculations necessary to recover these quantities. Notice that it is not necessary to compute general (i,j) elements, but only the rows (i,0) and (i,1). The former is in fact the EDS associated with P. Like Shipsey’s algorithm, Stange presents a double-and-add algorithm for the general case of reaching elements (m,0) and (m,1) via successive blocks of elements. For details on this, a possible refinement of my own, and a SAGE implementation, see the next post. This would solve the first problem as well, by choosing a random Q, but the (i,1) terms introduce an unnecessary overhead. Eliminating this gives rise to an alternative to Shipsey’s algorithm for EDS, so I’ve implemented this case separately with an equivalent interface to my existing eds class. It’s unclear which approach is preferable where (although Stange’s is almost definitely the one to use over finite fields), so I hope to furnish both with equivalent functionality and compare their performance for height computation (details to follow). Finally, the wrappers to move from points to nets and compute pairings have been implemented, although I need a lot more test cases to identify potential shortcomings! Ultimately, I’d like to bundle all these ideas together into a more friendly, rank-aware elliptic net class for both roles (heights and pairings).

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

Tomorrow I give my first year presentation, which determines whether I’m allowed to continue my studies. At 20-30 minutes, it’s a rather condensed version of my recent geometry club talk on the point counting problem, although this time I’m skipping more quickly through the fundamentals so that I can discuss some of the algorithms in depth, and taking a (hopefully) clearer route.

Both the report itself and the OHP slides version are available (pdf). Content covered: hyperelliptic curves, points, divisors, mumford polynomials and the Picard group/Jacobian; the discrete logarithm problem; explicit group law computation; characteristic polynomial of Frobenius and Weil theorems/interval; group-theoretic approaches; Schoof’s algorithm, SEA in genus 1, genus 2 hybrid algorithms.

Today I spoke at the Geometry Club about the use of hyperelliptic curves in public key cryptography. You can find my slides here, although they were supplemented by some boardwork that you’ll have to figure out from my other postings!