1package org.bouncycastle.math.ec;
2
3import java.math.BigInteger;
4
5/**
6 * Class holding methods for point multiplication based on the window
7 * τ-adic nonadjacent form (WTNAF). The algorithms are based on the
8 * paper "Improved Algorithms for Arithmetic on Anomalous Binary Curves"
9 * by Jerome A. Solinas. The paper first appeared in the Proceedings of
10 * Crypto 1997.
11 */
12class Tnaf
13{
14    private static final BigInteger MINUS_ONE = ECConstants.ONE.negate();
15    private static final BigInteger MINUS_TWO = ECConstants.TWO.negate();
16    private static final BigInteger MINUS_THREE = ECConstants.THREE.negate();
17
18    /**
19     * The window width of WTNAF. The standard value of 4 is slightly less
20     * than optimal for running time, but keeps space requirements for
21     * precomputation low. For typical curves, a value of 5 or 6 results in
22     * a better running time. When changing this value, the
23     * <code>&alpha;<sub>u</sub></code>'s must be computed differently, see
24     * e.g. "Guide to Elliptic Curve Cryptography", Darrel Hankerson,
25     * Alfred Menezes, Scott Vanstone, Springer-Verlag New York Inc., 2004,
26     * p. 121-122
27     */
28    public static final byte WIDTH = 4;
29
30    /**
31     * 2<sup>4</sup>
32     */
33    public static final byte POW_2_WIDTH = 16;
34
35    /**
36     * The <code>&alpha;<sub>u</sub></code>'s for <code>a=0</code> as an array
37     * of <code>ZTauElement</code>s.
38     */
39    public static final ZTauElement[] alpha0 = {
40        null,
41        new ZTauElement(ECConstants.ONE, ECConstants.ZERO), null,
42        new ZTauElement(MINUS_THREE, MINUS_ONE), null,
43        new ZTauElement(MINUS_ONE, MINUS_ONE), null,
44        new ZTauElement(ECConstants.ONE, MINUS_ONE), null
45    };
46
47    /**
48     * The <code>&alpha;<sub>u</sub></code>'s for <code>a=0</code> as an array
49     * of TNAFs.
50     */
51    public static final byte[][] alpha0Tnaf = {
52        null, {1}, null, {-1, 0, 1}, null, {1, 0, 1}, null, {-1, 0, 0, 1}
53    };
54
55    /**
56     * The <code>&alpha;<sub>u</sub></code>'s for <code>a=1</code> as an array
57     * of <code>ZTauElement</code>s.
58     */
59    public static final ZTauElement[] alpha1 = {null,
60        new ZTauElement(ECConstants.ONE, ECConstants.ZERO), null,
61        new ZTauElement(MINUS_THREE, ECConstants.ONE), null,
62        new ZTauElement(MINUS_ONE, ECConstants.ONE), null,
63        new ZTauElement(ECConstants.ONE, ECConstants.ONE), null
64    };
65
66    /**
67     * The <code>&alpha;<sub>u</sub></code>'s for <code>a=1</code> as an array
68     * of TNAFs.
69     */
70    public static final byte[][] alpha1Tnaf = {
71        null, {1}, null, {-1, 0, 1}, null, {1, 0, 1}, null, {-1, 0, 0, -1}
72    };
73
74    /**
75     * Computes the norm of an element <code>&lambda;</code> of
76     * <code><b>Z</b>[&tau;]</code>.
77     * @param mu The parameter <code>&mu;</code> of the elliptic curve.
78     * @param lambda The element <code>&lambda;</code> of
79     * <code><b>Z</b>[&tau;]</code>.
80     * @return The norm of <code>&lambda;</code>.
81     */
82    public static BigInteger norm(final byte mu, ZTauElement lambda)
83    {
84        BigInteger norm;
85
86        // s1 = u^2
87        BigInteger s1 = lambda.u.multiply(lambda.u);
88
89        // s2 = u * v
90        BigInteger s2 = lambda.u.multiply(lambda.v);
91
92        // s3 = 2 * v^2
93        BigInteger s3 = lambda.v.multiply(lambda.v).shiftLeft(1);
94
95        if (mu == 1)
96        {
97            norm = s1.add(s2).add(s3);
98        }
99        else if (mu == -1)
100        {
101            norm = s1.subtract(s2).add(s3);
102        }
103        else
104        {
105            throw new IllegalArgumentException("mu must be 1 or -1");
106        }
107
108        return norm;
109    }
110
111    /**
112     * Computes the norm of an element <code>&lambda;</code> of
113     * <code><b>R</b>[&tau;]</code>, where <code>&lambda; = u + v&tau;</code>
114     * and <code>u</code> and <code>u</code> are real numbers (elements of
115     * <code><b>R</b></code>).
116     * @param mu The parameter <code>&mu;</code> of the elliptic curve.
117     * @param u The real part of the element <code>&lambda;</code> of
118     * <code><b>R</b>[&tau;]</code>.
119     * @param v The <code>&tau;</code>-adic part of the element
120     * <code>&lambda;</code> of <code><b>R</b>[&tau;]</code>.
121     * @return The norm of <code>&lambda;</code>.
122     */
123    public static SimpleBigDecimal norm(final byte mu, SimpleBigDecimal u,
124            SimpleBigDecimal v)
125    {
126        SimpleBigDecimal norm;
127
128        // s1 = u^2
129        SimpleBigDecimal s1 = u.multiply(u);
130
131        // s2 = u * v
132        SimpleBigDecimal s2 = u.multiply(v);
133
134        // s3 = 2 * v^2
135        SimpleBigDecimal s3 = v.multiply(v).shiftLeft(1);
136
137        if (mu == 1)
138        {
139            norm = s1.add(s2).add(s3);
140        }
141        else if (mu == -1)
142        {
143            norm = s1.subtract(s2).add(s3);
144        }
145        else
146        {
147            throw new IllegalArgumentException("mu must be 1 or -1");
148        }
149
150        return norm;
151    }
152
153    /**
154     * Rounds an element <code>&lambda;</code> of <code><b>R</b>[&tau;]</code>
155     * to an element of <code><b>Z</b>[&tau;]</code>, such that their difference
156     * has minimal norm. <code>&lambda;</code> is given as
157     * <code>&lambda; = &lambda;<sub>0</sub> + &lambda;<sub>1</sub>&tau;</code>.
158     * @param lambda0 The component <code>&lambda;<sub>0</sub></code>.
159     * @param lambda1 The component <code>&lambda;<sub>1</sub></code>.
160     * @param mu The parameter <code>&mu;</code> of the elliptic curve. Must
161     * equal 1 or -1.
162     * @return The rounded element of <code><b>Z</b>[&tau;]</code>.
163     * @throws IllegalArgumentException if <code>lambda0</code> and
164     * <code>lambda1</code> do not have same scale.
165     */
166    public static ZTauElement round(SimpleBigDecimal lambda0,
167            SimpleBigDecimal lambda1, byte mu)
168    {
169        int scale = lambda0.getScale();
170        if (lambda1.getScale() != scale)
171        {
172            throw new IllegalArgumentException("lambda0 and lambda1 do not " +
173                    "have same scale");
174        }
175
176        if (!((mu == 1) || (mu == -1)))
177        {
178            throw new IllegalArgumentException("mu must be 1 or -1");
179        }
180
181        BigInteger f0 = lambda0.round();
182        BigInteger f1 = lambda1.round();
183
184        SimpleBigDecimal eta0 = lambda0.subtract(f0);
185        SimpleBigDecimal eta1 = lambda1.subtract(f1);
186
187        // eta = 2*eta0 + mu*eta1
188        SimpleBigDecimal eta = eta0.add(eta0);
189        if (mu == 1)
190        {
191            eta = eta.add(eta1);
192        }
193        else
194        {
195            // mu == -1
196            eta = eta.subtract(eta1);
197        }
198
199        // check1 = eta0 - 3*mu*eta1
200        // check2 = eta0 + 4*mu*eta1
201        SimpleBigDecimal threeEta1 = eta1.add(eta1).add(eta1);
202        SimpleBigDecimal fourEta1 = threeEta1.add(eta1);
203        SimpleBigDecimal check1;
204        SimpleBigDecimal check2;
205        if (mu == 1)
206        {
207            check1 = eta0.subtract(threeEta1);
208            check2 = eta0.add(fourEta1);
209        }
210        else
211        {
212            // mu == -1
213            check1 = eta0.add(threeEta1);
214            check2 = eta0.subtract(fourEta1);
215        }
216
217        byte h0 = 0;
218        byte h1 = 0;
219
220        // if eta >= 1
221        if (eta.compareTo(ECConstants.ONE) >= 0)
222        {
223            if (check1.compareTo(MINUS_ONE) < 0)
224            {
225                h1 = mu;
226            }
227            else
228            {
229                h0 = 1;
230            }
231        }
232        else
233        {
234            // eta < 1
235            if (check2.compareTo(ECConstants.TWO) >= 0)
236            {
237                h1 = mu;
238            }
239        }
240
241        // if eta < -1
242        if (eta.compareTo(MINUS_ONE) < 0)
243        {
244            if (check1.compareTo(ECConstants.ONE) >= 0)
245            {
246                h1 = (byte)-mu;
247            }
248            else
249            {
250                h0 = -1;
251            }
252        }
253        else
254        {
255            // eta >= -1
256            if (check2.compareTo(MINUS_TWO) < 0)
257            {
258                h1 = (byte)-mu;
259            }
260        }
261
262        BigInteger q0 = f0.add(BigInteger.valueOf(h0));
263        BigInteger q1 = f1.add(BigInteger.valueOf(h1));
264        return new ZTauElement(q0, q1);
265    }
266
267    /**
268     * Approximate division by <code>n</code>. For an integer
269     * <code>k</code>, the value <code>&lambda; = s k / n</code> is
270     * computed to <code>c</code> bits of accuracy.
271     * @param k The parameter <code>k</code>.
272     * @param s The curve parameter <code>s<sub>0</sub></code> or
273     * <code>s<sub>1</sub></code>.
274     * @param vm The Lucas Sequence element <code>V<sub>m</sub></code>.
275     * @param a The parameter <code>a</code> of the elliptic curve.
276     * @param m The bit length of the finite field
277     * <code><b>F</b><sub>m</sub></code>.
278     * @param c The number of bits of accuracy, i.e. the scale of the returned
279     * <code>SimpleBigDecimal</code>.
280     * @return The value <code>&lambda; = s k / n</code> computed to
281     * <code>c</code> bits of accuracy.
282     */
283    public static SimpleBigDecimal approximateDivisionByN(BigInteger k,
284            BigInteger s, BigInteger vm, byte a, int m, int c)
285    {
286        int _k = (m + 5)/2 + c;
287        BigInteger ns = k.shiftRight(m - _k - 2 + a);
288
289        BigInteger gs = s.multiply(ns);
290
291        BigInteger hs = gs.shiftRight(m);
292
293        BigInteger js = vm.multiply(hs);
294
295        BigInteger gsPlusJs = gs.add(js);
296        BigInteger ls = gsPlusJs.shiftRight(_k-c);
297        if (gsPlusJs.testBit(_k-c-1))
298        {
299            // round up
300            ls = ls.add(ECConstants.ONE);
301        }
302
303        return new SimpleBigDecimal(ls, c);
304    }
305
306    /**
307     * Computes the <code>&tau;</code>-adic NAF (non-adjacent form) of an
308     * element <code>&lambda;</code> of <code><b>Z</b>[&tau;]</code>.
309     * @param mu The parameter <code>&mu;</code> of the elliptic curve.
310     * @param lambda The element <code>&lambda;</code> of
311     * <code><b>Z</b>[&tau;]</code>.
312     * @return The <code>&tau;</code>-adic NAF of <code>&lambda;</code>.
313     */
314    public static byte[] tauAdicNaf(byte mu, ZTauElement lambda)
315    {
316        if (!((mu == 1) || (mu == -1)))
317        {
318            throw new IllegalArgumentException("mu must be 1 or -1");
319        }
320
321        BigInteger norm = norm(mu, lambda);
322
323        // Ceiling of log2 of the norm
324        int log2Norm = norm.bitLength();
325
326        // If length(TNAF) > 30, then length(TNAF) < log2Norm + 3.52
327        int maxLength = log2Norm > 30 ? log2Norm + 4 : 34;
328
329        // The array holding the TNAF
330        byte[] u = new byte[maxLength];
331        int i = 0;
332
333        // The actual length of the TNAF
334        int length = 0;
335
336        BigInteger r0 = lambda.u;
337        BigInteger r1 = lambda.v;
338
339        while(!((r0.equals(ECConstants.ZERO)) && (r1.equals(ECConstants.ZERO))))
340        {
341            // If r0 is odd
342            if (r0.testBit(0))
343            {
344                u[i] = (byte) ECConstants.TWO.subtract((r0.subtract(r1.shiftLeft(1))).mod(ECConstants.FOUR)).intValue();
345
346                // r0 = r0 - u[i]
347                if (u[i] == 1)
348                {
349                    r0 = r0.clearBit(0);
350                }
351                else
352                {
353                    // u[i] == -1
354                    r0 = r0.add(ECConstants.ONE);
355                }
356                length = i;
357            }
358            else
359            {
360                u[i] = 0;
361            }
362
363            BigInteger t = r0;
364            BigInteger s = r0.shiftRight(1);
365            if (mu == 1)
366            {
367                r0 = r1.add(s);
368            }
369            else
370            {
371                // mu == -1
372                r0 = r1.subtract(s);
373            }
374
375            r1 = t.shiftRight(1).negate();
376            i++;
377        }
378
379        length++;
380
381        // Reduce the TNAF array to its actual length
382        byte[] tnaf = new byte[length];
383        System.arraycopy(u, 0, tnaf, 0, length);
384        return tnaf;
385    }
386
387    /**
388     * Applies the operation <code>&tau;()</code> to an
389     * <code>ECPoint.F2m</code>.
390     * @param p The ECPoint.F2m to which <code>&tau;()</code> is applied.
391     * @return <code>&tau;(p)</code>
392     */
393    public static ECPoint.F2m tau(ECPoint.F2m p)
394    {
395        return p.tau();
396    }
397
398    /**
399     * Returns the parameter <code>&mu;</code> of the elliptic curve.
400     * @param curve The elliptic curve from which to obtain <code>&mu;</code>.
401     * The curve must be a Koblitz curve, i.e. <code>a</code> equals
402     * <code>0</code> or <code>1</code> and <code>b</code> equals
403     * <code>1</code>.
404     * @return <code>&mu;</code> of the elliptic curve.
405     * @throws IllegalArgumentException if the given ECCurve is not a Koblitz
406     * curve.
407     */
408    public static byte getMu(ECCurve.F2m curve)
409    {
410        if (!curve.isKoblitz())
411        {
412            throw new IllegalArgumentException("No Koblitz curve (ABC), TNAF multiplication not possible");
413        }
414
415        if (curve.getA().isZero())
416        {
417            return -1;
418        }
419
420        return 1;
421    }
422
423    /**
424     * Calculates the Lucas Sequence elements <code>U<sub>k-1</sub></code> and
425     * <code>U<sub>k</sub></code> or <code>V<sub>k-1</sub></code> and
426     * <code>V<sub>k</sub></code>.
427     * @param mu The parameter <code>&mu;</code> of the elliptic curve.
428     * @param k The index of the second element of the Lucas Sequence to be
429     * returned.
430     * @param doV If set to true, computes <code>V<sub>k-1</sub></code> and
431     * <code>V<sub>k</sub></code>, otherwise <code>U<sub>k-1</sub></code> and
432     * <code>U<sub>k</sub></code>.
433     * @return An array with 2 elements, containing <code>U<sub>k-1</sub></code>
434     * and <code>U<sub>k</sub></code> or <code>V<sub>k-1</sub></code>
435     * and <code>V<sub>k</sub></code>.
436     */
437    public static BigInteger[] getLucas(byte mu, int k, boolean doV)
438    {
439        if (!((mu == 1) || (mu == -1)))
440        {
441            throw new IllegalArgumentException("mu must be 1 or -1");
442        }
443
444        BigInteger u0;
445        BigInteger u1;
446        BigInteger u2;
447
448        if (doV)
449        {
450            u0 = ECConstants.TWO;
451            u1 = BigInteger.valueOf(mu);
452        }
453        else
454        {
455            u0 = ECConstants.ZERO;
456            u1 = ECConstants.ONE;
457        }
458
459        for (int i = 1; i < k; i++)
460        {
461            // u2 = mu*u1 - 2*u0;
462            BigInteger s = null;
463            if (mu == 1)
464            {
465                s = u1;
466            }
467            else
468            {
469                // mu == -1
470                s = u1.negate();
471            }
472
473            u2 = s.subtract(u0.shiftLeft(1));
474            u0 = u1;
475            u1 = u2;
476//            System.out.println(i + ": " + u2);
477//            System.out.println();
478        }
479
480        BigInteger[] retVal = {u0, u1};
481        return retVal;
482    }
483
484    /**
485     * Computes the auxiliary value <code>t<sub>w</sub></code>. If the width is
486     * 4, then for <code>mu = 1</code>, <code>t<sub>w</sub> = 6</code> and for
487     * <code>mu = -1</code>, <code>t<sub>w</sub> = 10</code>
488     * @param mu The parameter <code>&mu;</code> of the elliptic curve.
489     * @param w The window width of the WTNAF.
490     * @return the auxiliary value <code>t<sub>w</sub></code>
491     */
492    public static BigInteger getTw(byte mu, int w)
493    {
494        if (w == 4)
495        {
496            if (mu == 1)
497            {
498                return BigInteger.valueOf(6);
499            }
500            else
501            {
502                // mu == -1
503                return BigInteger.valueOf(10);
504            }
505        }
506        else
507        {
508            // For w <> 4, the values must be computed
509            BigInteger[] us = getLucas(mu, w, false);
510            BigInteger twoToW = ECConstants.ZERO.setBit(w);
511            BigInteger u1invert = us[1].modInverse(twoToW);
512            BigInteger tw;
513            tw = ECConstants.TWO.multiply(us[0]).multiply(u1invert).mod(twoToW);
514//            System.out.println("mu = " + mu);
515//            System.out.println("tw = " + tw);
516            return tw;
517        }
518    }
519
520    /**
521     * Computes the auxiliary values <code>s<sub>0</sub></code> and
522     * <code>s<sub>1</sub></code> used for partial modular reduction.
523     * @param curve The elliptic curve for which to compute
524     * <code>s<sub>0</sub></code> and <code>s<sub>1</sub></code>.
525     * @throws IllegalArgumentException if <code>curve</code> is not a
526     * Koblitz curve (Anomalous Binary Curve, ABC).
527     */
528    public static BigInteger[] getSi(ECCurve.F2m curve)
529    {
530        if (!curve.isKoblitz())
531        {
532            throw new IllegalArgumentException("si is defined for Koblitz curves only");
533        }
534
535        int m = curve.getM();
536        int a = curve.getA().toBigInteger().intValue();
537        byte mu = curve.getMu();
538        int h = curve.getH().intValue();
539        int index = m + 3 - a;
540        BigInteger[] ui = getLucas(mu, index, false);
541
542        BigInteger dividend0;
543        BigInteger dividend1;
544        if (mu == 1)
545        {
546            dividend0 = ECConstants.ONE.subtract(ui[1]);
547            dividend1 = ECConstants.ONE.subtract(ui[0]);
548        }
549        else if (mu == -1)
550        {
551            dividend0 = ECConstants.ONE.add(ui[1]);
552            dividend1 = ECConstants.ONE.add(ui[0]);
553        }
554        else
555        {
556            throw new IllegalArgumentException("mu must be 1 or -1");
557        }
558
559        BigInteger[] si = new BigInteger[2];
560
561        if (h == 2)
562        {
563            si[0] = dividend0.shiftRight(1);
564            si[1] = dividend1.shiftRight(1).negate();
565        }
566        else if (h == 4)
567        {
568            si[0] = dividend0.shiftRight(2);
569            si[1] = dividend1.shiftRight(2).negate();
570        }
571        else
572        {
573            throw new IllegalArgumentException("h (Cofactor) must be 2 or 4");
574        }
575
576        return si;
577    }
578
579    /**
580     * Partial modular reduction modulo
581     * <code>(&tau;<sup>m</sup> - 1)/(&tau; - 1)</code>.
582     * @param k The integer to be reduced.
583     * @param m The bitlength of the underlying finite field.
584     * @param a The parameter <code>a</code> of the elliptic curve.
585     * @param s The auxiliary values <code>s<sub>0</sub></code> and
586     * <code>s<sub>1</sub></code>.
587     * @param mu The parameter &mu; of the elliptic curve.
588     * @param c The precision (number of bits of accuracy) of the partial
589     * modular reduction.
590     * @return <code>&rho; := k partmod (&tau;<sup>m</sup> - 1)/(&tau; - 1)</code>
591     */
592    public static ZTauElement partModReduction(BigInteger k, int m, byte a,
593            BigInteger[] s, byte mu, byte c)
594    {
595        // d0 = s[0] + mu*s[1]; mu is either 1 or -1
596        BigInteger d0;
597        if (mu == 1)
598        {
599            d0 = s[0].add(s[1]);
600        }
601        else
602        {
603            d0 = s[0].subtract(s[1]);
604        }
605
606        BigInteger[] v = getLucas(mu, m, true);
607        BigInteger vm = v[1];
608
609        SimpleBigDecimal lambda0 = approximateDivisionByN(
610                k, s[0], vm, a, m, c);
611
612        SimpleBigDecimal lambda1 = approximateDivisionByN(
613                k, s[1], vm, a, m, c);
614
615        ZTauElement q = round(lambda0, lambda1, mu);
616
617        // r0 = n - d0*q0 - 2*s1*q1
618        BigInteger r0 = k.subtract(d0.multiply(q.u)).subtract(
619                BigInteger.valueOf(2).multiply(s[1]).multiply(q.v));
620
621        // r1 = s1*q0 - s0*q1
622        BigInteger r1 = s[1].multiply(q.u).subtract(s[0].multiply(q.v));
623
624        return new ZTauElement(r0, r1);
625    }
626
627    /**
628     * Multiplies a {@link org.bouncycastle.math.ec.ECPoint.F2m ECPoint.F2m}
629     * by a <code>BigInteger</code> using the reduced <code>&tau;</code>-adic
630     * NAF (RTNAF) method.
631     * @param p The ECPoint.F2m to multiply.
632     * @param k The <code>BigInteger</code> by which to multiply <code>p</code>.
633     * @return <code>k * p</code>
634     */
635    public static ECPoint.F2m multiplyRTnaf(ECPoint.F2m p, BigInteger k)
636    {
637        ECCurve.F2m curve = (ECCurve.F2m) p.getCurve();
638        int m = curve.getM();
639        byte a = (byte) curve.getA().toBigInteger().intValue();
640        byte mu = curve.getMu();
641        BigInteger[] s = curve.getSi();
642        ZTauElement rho = partModReduction(k, m, a, s, mu, (byte)10);
643
644        return multiplyTnaf(p, rho);
645    }
646
647    /**
648     * Multiplies a {@link org.bouncycastle.math.ec.ECPoint.F2m ECPoint.F2m}
649     * by an element <code>&lambda;</code> of <code><b>Z</b>[&tau;]</code>
650     * using the <code>&tau;</code>-adic NAF (TNAF) method.
651     * @param p The ECPoint.F2m to multiply.
652     * @param lambda The element <code>&lambda;</code> of
653     * <code><b>Z</b>[&tau;]</code>.
654     * @return <code>&lambda; * p</code>
655     */
656    public static ECPoint.F2m multiplyTnaf(ECPoint.F2m p, ZTauElement lambda)
657    {
658        ECCurve.F2m curve = (ECCurve.F2m)p.getCurve();
659        byte mu = curve.getMu();
660        byte[] u = tauAdicNaf(mu, lambda);
661
662        ECPoint.F2m q = multiplyFromTnaf(p, u);
663
664        return q;
665    }
666
667    /**
668    * Multiplies a {@link org.bouncycastle.math.ec.ECPoint.F2m ECPoint.F2m}
669    * by an element <code>&lambda;</code> of <code><b>Z</b>[&tau;]</code>
670    * using the <code>&tau;</code>-adic NAF (TNAF) method, given the TNAF
671    * of <code>&lambda;</code>.
672    * @param p The ECPoint.F2m to multiply.
673    * @param u The the TNAF of <code>&lambda;</code>..
674    * @return <code>&lambda; * p</code>
675    */
676    public static ECPoint.F2m multiplyFromTnaf(ECPoint.F2m p, byte[] u)
677    {
678        ECCurve.F2m curve = (ECCurve.F2m)p.getCurve();
679        ECPoint.F2m q = (ECPoint.F2m) curve.getInfinity();
680        for (int i = u.length - 1; i >= 0; i--)
681        {
682            q = tau(q);
683            if (u[i] == 1)
684            {
685                q = (ECPoint.F2m)q.addSimple(p);
686            }
687            else if (u[i] == -1)
688            {
689                q = (ECPoint.F2m)q.subtractSimple(p);
690            }
691        }
692        return q;
693    }
694
695    /**
696     * Computes the <code>[&tau;]</code>-adic window NAF of an element
697     * <code>&lambda;</code> of <code><b>Z</b>[&tau;]</code>.
698     * @param mu The parameter &mu; of the elliptic curve.
699     * @param lambda The element <code>&lambda;</code> of
700     * <code><b>Z</b>[&tau;]</code> of which to compute the
701     * <code>[&tau;]</code>-adic NAF.
702     * @param width The window width of the resulting WNAF.
703     * @param pow2w 2<sup>width</sup>.
704     * @param tw The auxiliary value <code>t<sub>w</sub></code>.
705     * @param alpha The <code>&alpha;<sub>u</sub></code>'s for the window width.
706     * @return The <code>[&tau;]</code>-adic window NAF of
707     * <code>&lambda;</code>.
708     */
709    public static byte[] tauAdicWNaf(byte mu, ZTauElement lambda,
710            byte width, BigInteger pow2w, BigInteger tw, ZTauElement[] alpha)
711    {
712        if (!((mu == 1) || (mu == -1)))
713        {
714            throw new IllegalArgumentException("mu must be 1 or -1");
715        }
716
717        BigInteger norm = norm(mu, lambda);
718
719        // Ceiling of log2 of the norm
720        int log2Norm = norm.bitLength();
721
722        // If length(TNAF) > 30, then length(TNAF) < log2Norm + 3.52
723        int maxLength = log2Norm > 30 ? log2Norm + 4 + width : 34 + width;
724
725        // The array holding the TNAF
726        byte[] u = new byte[maxLength];
727
728        // 2^(width - 1)
729        BigInteger pow2wMin1 = pow2w.shiftRight(1);
730
731        // Split lambda into two BigIntegers to simplify calculations
732        BigInteger r0 = lambda.u;
733        BigInteger r1 = lambda.v;
734        int i = 0;
735
736        // while lambda <> (0, 0)
737        while (!((r0.equals(ECConstants.ZERO))&&(r1.equals(ECConstants.ZERO))))
738        {
739            // if r0 is odd
740            if (r0.testBit(0))
741            {
742                // uUnMod = r0 + r1*tw mod 2^width
743                BigInteger uUnMod
744                    = r0.add(r1.multiply(tw)).mod(pow2w);
745
746                byte uLocal;
747                // if uUnMod >= 2^(width - 1)
748                if (uUnMod.compareTo(pow2wMin1) >= 0)
749                {
750                    uLocal = (byte) uUnMod.subtract(pow2w).intValue();
751                }
752                else
753                {
754                    uLocal = (byte) uUnMod.intValue();
755                }
756                // uLocal is now in [-2^(width-1), 2^(width-1)-1]
757
758                u[i] = uLocal;
759                boolean s = true;
760                if (uLocal < 0)
761                {
762                    s = false;
763                    uLocal = (byte)-uLocal;
764                }
765                // uLocal is now >= 0
766
767                if (s)
768                {
769                    r0 = r0.subtract(alpha[uLocal].u);
770                    r1 = r1.subtract(alpha[uLocal].v);
771                }
772                else
773                {
774                    r0 = r0.add(alpha[uLocal].u);
775                    r1 = r1.add(alpha[uLocal].v);
776                }
777            }
778            else
779            {
780                u[i] = 0;
781            }
782
783            BigInteger t = r0;
784
785            if (mu == 1)
786            {
787                r0 = r1.add(r0.shiftRight(1));
788            }
789            else
790            {
791                // mu == -1
792                r0 = r1.subtract(r0.shiftRight(1));
793            }
794            r1 = t.shiftRight(1).negate();
795            i++;
796        }
797        return u;
798    }
799
800    /**
801     * Does the precomputation for WTNAF multiplication.
802     * @param p The <code>ECPoint</code> for which to do the precomputation.
803     * @param a The parameter <code>a</code> of the elliptic curve.
804     * @return The precomputation array for <code>p</code>.
805     */
806    public static ECPoint.F2m[] getPreComp(ECPoint.F2m p, byte a)
807    {
808        ECPoint.F2m[] pu;
809        pu = new ECPoint.F2m[16];
810        pu[1] = p;
811        byte[][] alphaTnaf;
812        if (a == 0)
813        {
814            alphaTnaf = Tnaf.alpha0Tnaf;
815        }
816        else
817        {
818            // a == 1
819            alphaTnaf = Tnaf.alpha1Tnaf;
820        }
821
822        int precompLen = alphaTnaf.length;
823        for (int i = 3; i < precompLen; i = i + 2)
824        {
825            pu[i] = Tnaf.multiplyFromTnaf(p, alphaTnaf[i]);
826        }
827
828        p.getCurve().normalizeAll(pu);
829
830        return pu;
831    }
832}
833