for when the margin is too small
As a companion to my previous post, here’s the list of valid forms of a connected maximal cyclotomic graph with an entry from the ring of integers of Q(√ -7):
A chain of the form
for any integer k.
or
A chain of the form
for any integer k.
To recap: I’ve been trying to completely classify the possible matrices/graphs subject to a constraint on their eigenvalues we’re describing as cyclotomicity. This is a problem that can be posed in the ring of integers of any imaginary quadratic extension field, but for all but finitely many of them reduces to the problem in the rational-integer case which has been solved in a paper by my supervisor.
For a couple of the remaining fields, the problem is easy: there’s only a finite supply of graphs featuring a non-rational integer label, which can be found simply by running a growing process to termination. But once you move to fields with norm 2 integers, there’s enough freedom for things to get interesting: I’ve been working in the rings of integers of Q(√ -2) and Q(√ -7), where I can demonstrate an infinite family of such graphs and so the growing algorithm can never terminate. Nonetheless, in the simpler, ‘uncharged’ version of the problem, I have (proven) a complete classification in both of these fields. With a bit more work I’ve now settled the general case in Q(√ -2) and expect the logic of the argument (although not the precise results) to carry over to Q(√ -7).
That argument is essentially a lengthy case analysis; rather than get into the details I thought I’d just present the ‘zoo’ of possible graphs. The forms presented are necessary conditions (any cyclotomic graph will take one of these forms) but not sufficient (there may be non-cyclotomic graphs satisfying the form). However, no form contains no cyclotomic graphs - for a given form, including any instance of the infinite ones, I can exhibit at least one class of cyclotomic graphs.
After applying a numbering, the visual styling of an edge between two nodes i<j indicates the norm of the edge label (entry [i,j] of the matrix; take the conjugate for entry [j,i]); uncharged nodes are indicated by a point whilst [C] denotes a charged node (value of ±1 for entry [i.i] if node i is charged, otherwise zero). The precise choice of labels and charges requires some care over signs to ensure that the matrix has minimal polynomial x2-4; working with forms saves tracking such details, choosing equivalence class representatives, etc.
Then for the ring of integers of Q(√ -2) any connected maximal cyclotomic graph with a non-rational integer label must take one of the following forms:
A chain of the form
for any integer k.
or
A chain of the form
for any integer k.
As I said, I expect this to easily generalise to Q(√-7); the remaining fields Q(√-1) and Q(√-3) present more of a computational challenge, but meeting that challenge will hopefully be rewarded with more interesting behaviour!
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.
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:
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.
(updated 10/i/08)
In an earlier post I described Stange’s algorithm for efficiently finding terms in elliptic nets (with a view to pairing computation). I also made the observation that a shorter block structure could be used for doubling- but once employed, it was not possible to perform a double-and-add. This meant that unless the desired term had a higher power of two as a factor then savings would be minor.
However, for a cost it is possible to ‘upgrade’ these short blocks to long blocks, since they contain enough information to recover the missing (k+4,0) term:
Better still, this only introduces an additional two multiplications and one inversion, since some of the terms feature in the precomputation (and assuming (2,0)2 is computed once and stored for subsequent use):
Thus, given a short block centred at k, we can obtain the short block centred at 2k+1 (that is, perform a shortDoubleAdd). The dependencies are as follows:
Notice that a shortDoubleAdd is more expensive than DoubleAdd, even though it gives a short rather than long block! Thus a purely short-block algorithm would perform worse than the standard algorithm for binary strings with a high Hamming weight, since for each Double-and-add an inversion is introduced in place of a multiplication. However, when the Hamming weight is low, then the occasional cost of an inversion is balanced by the savings accrued during short doublings. To exploit this, whilst guarding against too many inversions, we introduce an algorithm that uses a mixture of standard (‘long’) and short blocks. Since only a single inversion is required to switch to long blocks, doing so allows us to more efficiently compute long runs of 1s in the binary representation by spreading the cost across several DoubleAdds.
We consider the generation of long or short blocks with centre 2k (double) or 2k + 1 (double-and-add) from long or short blocks of centre k. The cheapest such operation is the generation of the short block with centre 2k from a short or long block with centre k, at a cost of 31 multiplications, 6 squarings and no inversions. Using this as a base line, each procedure introduces the following additional operations:
| Procedure | M | S | I |
|---|---|---|---|
| DoubleShortFromShort | 0 | 0 | 0 |
| DoubleLongFromShort | 6 | 1 | 1 |
| DoubleAddShortFromShort | 4 | 1 | 1 |
| DoubleAddLongFromShort | 6 | 1 | 1 |
| DoubleShortFromLong | 0 | 0 | 0 |
| DoubleLongFromLong | 4 | 1 | 0 |
| DoubleAddShortFromLong | 2 | 1 | 0 |
| DoubleAddLongFromLong | 4 | 1 | 0 |
We adopt a windowing approach with two-bit windows: that is, bit bi informs whether we are to double or double-and-add, but bi-1 is also examined to determine whether we should generate a long or short block.
Hence the approach is to always go to (or stay with, if already the case) short blocks, unless bibi-1=11 in which case one should go to (or stay with) long blocks. Thus a 2-bit window is sufficient to determine appropriate block length, leading to the following algorithm.
Double-and-add Mixed-blocks Algorithm
INPUT: Integer n and long block centred at 1.
OUTPUT: Block centred at n.
As described before, the maximum possible gain is when n is a power of two, in which case the algorithm proceeds entirely by short doubles. In this case, there is a 12 percent reduction in the number of multiplications/squarings performed, with no inversions required.
Brute-force analysis of all possible 16-bit strings gives an average reduction of around 9 percent in the number of multiplications/squarings performed. Costing each inversion at 1.6 multiplications (based on average performance in SAGE for a 256-bit prime field), this leads to an average reduction of around 5 percent in the number of multiplications required for such strings. Testing several hundred 256-bit strings gives a similar figure.
Clearly, inversion is not viable if it will lead to a division-by-zero error. However, since the first zero along (i,0) will arise at (m,0), no such error will occur when performing Tate pairing computations.
A summary of this post and the earlier one on Stange’s algorithm is available as pdf, containing (hopefully) clearer copies of the dependency graphs.
Suppose you have a rule for addition of some objects (most generally, elements of a semigroup), and you wish to compute the sum of n copies of the same object. How could you achieve this?
The primary school answer is to do just that: given g, compute 2g=g+g; then compute 3g=2g+g and so on, for n steps. But this rapidly becomes tedious, and isn’t (I hope) how we’d perform such a calculation in our heads: instead, we’d try to take bigger jumps. For instance, to compute the number 16g, we can find 2g, add that to itself to reach 4g, and again for 8g, hitting 16g after just 4 steps. Had we wanted 17, we’d then just add g again.
The sequence of intermediate values correspond to an addition chain for n: this is a sequence starting at a0=1 and ending at al=n with the property that for any term ak there exist terms ai,aj with ak=ai+aj. That is, every term is the sum of two earlier terms (not necessarily distinct), so we can reach ng by computing each akg in turn without needing anything fancier than our addition law.
As a side effect, an addition chain also tells us how to exponentiate if we know how to multiply, since xa+b=xaxb so we can build up to xn by finding the powers xak in turn.
Thus our tedious sequence for 16 is 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16, with the more sophisticated one being 1,2,4,8,16. If I asked you to compute 16x for a number x, you’d probably do something different: compute 10x (easy in base 10!) and 6x, then add those. With access to a multiplication procedure, this is definitely more efficient: but without that luxury it can still be used to give an addition chain for 16, by finding chains for 10 and 6. These are 1,2,4,5,10 and 1,2,3,6 respectively, so we can merge them into the chain 1,2,3,4,5,6,10,16.
We describe the length of the chain as l: that is, there are l terms after the 1. Notice that at worst, therefore, we need n-1 terms. If we move as quickly as possible by taking ak+1=ak+ak=2ak then the greatest value we can reach after m terms is 2m, as demonstrated with the chain for 24=16. However, if we can compute n and m in ln,lm terms, it needn’t take as many as ln+lm+1 terms to compute n+m, since there may be common terms in the chains. Finally, as defined an addition chain needn’t make use of all it’s terms: for example, 1,2,4,6,8,16 is still a chain for 16, although the 6 isn’t needed. We will obviously be interested in chains without any such wasted terms.
The observation for powers of two motivates a first attempt at an efficient algorithm for building addition chains. if n is 1, then we have nothing to do: the chain 1 suffices. Otherwise, we can work recursively: for an even number n=2k, we ask for the chain for k and then add 2k to the end; whereas for an odd number n=2k+1, we ask for the chain for k and add 2k then 2k+1 to the end. This will require at most log2n calls, and each step adds 1 or 2 entries to the chain.
This can also be captured by considering the binary representation of n, which will also indicate whether we hit an even or odd number at each stage. Suppose n is of the form
where each di is either 0 or 1, and dk=1 so there are no leading zeros. Then we start the chain with 1, and read from left to right, ignoring the leading 1: if we see a 0, then we double the last term of the chain; if we see a 1, add 1 to this new term as well. More formally:
Input: the binary expansion of n.
Output: an addition chain for n.
Set a0=1
Set t=1
for i from k-1 down to 0:
Set at=2at-1 and t=t+1.
If d_i=1:
Set at=at-1+1 and t=t+1
return a0…at
For instance, with n a power of 2 we will see only a run of zeros, and thus perform the appropriate number of doublings, for a chain of length log2n. For n one less than a power of two, the expansion consists entirely of 1s, costing 2log2n. For an ‘average’ number where each digit has around a 50% chance of being a 1 or a 0, the chain will therefore be of approximate length (3/2)log2n. More formally, we can define the Hamming Weight v(n) to be the number of 1s in the binary expansion of n; then the chain will have length floor(log2n) + v(n) -1.
Is this optimal? A counterexample suffices to show it is not, the first being 15, which is 1111 in binary. Using the binary expansion we get the sequence 1,2,3,6,7,14,15 of length 6 (=3+4-1 as claimed). But the sequence would be 1,2,3,5,10,15 of length 5 and hence 1 shorter.
Nonetheless, our binary chain gives a much better upper bound than the naive chain of (n-1) steps, and is very simple to implement both in logic and memory. The idea also generalises to arbitrary bases, as follows.
Suppose n is of the form
where each di is in the range 0…m-1, dk non-zero. Then we can compute an m-ary addition chain for n:
Input: the m-ary expansion of n.
Output: an addition chain for n.
Set a=0
Start the chain with 1,2,3,…,m-1.
for i from k down to 0:
Extend the chain to m*a (if not already present) and set a=m*a.
Add a+di to the chain (if not already present), and set a=a+di.
return chain.
Notice that a tracks the last term computed, and that extending the chain to m*a is itself an addition chain problem: this makes some 2k a good choice for m, since this extension can then be accomplished in k additions.
The shorter chain for 15 arises from its ternary expansion: considered as (120)3 1*9+2*3+0*1, we have an initial chain of 1,2 then apply the loop. With a=0,d2=1 we set a=3a=0 then a=a+1=1; neither of these need to be added to the chain; for a=1, d1=2 we set a=3a=3 which requires the addition of 3 to the chain; then a=a+2=5 which extends the chain to 1,2,3,5; finally for a=5,d2=0 we set a=3a=15, adding first 10 then 15 to the chain; then a=a+0=15 so we terminate with chain 1,2,3,5,10,15.
The optimal choice of m depends on n, since although there will be less steps (logmn) for greater m, each step will require more additions in the computation of m*a. Nor are m-ary expansions the end of the story- there are various other tricks available to finesse addition chains. However, there is no known algorithm for finding shortest addition chains- and a generalisation of the problem to finding shortest addition sequences (chains that contain a desired list of values n1,n2,…,nk) is NP-complete. Searching for short chains is therefore an interesting challenge, and is constrained by the fact that unless you are precomputing them for future use, any shorter chain must offer a speedup in computation greater than the time taken to find it.
There are, however, bounds for the length of the chain:
Where v(n) is the Hamming weight as before.
However, if you have access to basic operations other than addition, all bets are off. In some contexts (such as elliptic curves), finding the inverse of an element is easy and so we can achieve subtraction, giving rise to addition-subtraction chains. For the pairing computations I’ve been interested in, I’ve never actually described a block addition. It is possible, following an algorithm of Shipsey, to do so, although it is much more expensive than the block double. But even if it cost no more, and we had an oracle for optimal chains, such an approach would only rarely beat the current algorithm. This is because access to the double-add procedure makes the Hamming weight irrelevant: for instance, with 15 we could form the chain 1,3,7,15. This means that we can compute in log2n block operations, and comparing to the lower bound for addition chains above shows that this is almost always better.
Still, ideas from addition chains can motivate construction of other operations: access to triple, tripleadd and tripleadd2 procedures would allow for chain of log3n length, although it is likely that these block operations would be more expensive than the double and doubleadds. Even if they don’t turn out to be useful for my work, the subtleties of this innocent-looking problem still makes for interesting mathematics!
References/Further Reading
Knuth- The Art of Computer Programming Vol. 2: Seminumerical Algorithms.
Bos and Coster- Addition Chain Heuristics in Advances in Cryptology- Crypto’89 (Lecture Notes in Computer Science) (Available via SpringerLink.)
Gordon- A Survey of Fast Exponentiation Methods. (Available via citeseer.)
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.
Graeme Taylor
PhD Student,
Number Theory
About this blog