Make BN_GF2m_mod_arr more constant time
[openssl.git] / crypto / bn / bn_gf2m.c
1 /*
2  * Copyright 2002-2017 The OpenSSL Project Authors. All Rights Reserved.
3  * Copyright (c) 2002, Oracle and/or its affiliates. All rights reserved
4  *
5  * Licensed under the OpenSSL license (the "License").  You may not use
6  * this file except in compliance with the License.  You can obtain a copy
7  * in the file LICENSE in the source distribution or at
8  * https://www.openssl.org/source/license.html
9  */
10
11 #include <assert.h>
12 #include <limits.h>
13 #include <stdio.h>
14 #include "internal/cryptlib.h"
15 #include "bn_lcl.h"
16
17 #ifndef OPENSSL_NO_EC2M
18
19 /*
20  * Maximum number of iterations before BN_GF2m_mod_solve_quad_arr should
21  * fail.
22  */
23 # define MAX_ITERATIONS 50
24
25 # define SQR_nibble(w)   ((((w) & 8) << 3) \
26                        |  (((w) & 4) << 2) \
27                        |  (((w) & 2) << 1) \
28                        |   ((w) & 1))
29
30
31 /* Platform-specific macros to accelerate squaring. */
32 # if defined(SIXTY_FOUR_BIT) || defined(SIXTY_FOUR_BIT_LONG)
33 #  define SQR1(w) \
34     SQR_nibble((w) >> 60) << 56 | SQR_nibble((w) >> 56) << 48 | \
35     SQR_nibble((w) >> 52) << 40 | SQR_nibble((w) >> 48) << 32 | \
36     SQR_nibble((w) >> 44) << 24 | SQR_nibble((w) >> 40) << 16 | \
37     SQR_nibble((w) >> 36) <<  8 | SQR_nibble((w) >> 32)
38 #  define SQR0(w) \
39     SQR_nibble((w) >> 28) << 56 | SQR_nibble((w) >> 24) << 48 | \
40     SQR_nibble((w) >> 20) << 40 | SQR_nibble((w) >> 16) << 32 | \
41     SQR_nibble((w) >> 12) << 24 | SQR_nibble((w) >>  8) << 16 | \
42     SQR_nibble((w) >>  4) <<  8 | SQR_nibble((w)      )
43 # endif
44 # ifdef THIRTY_TWO_BIT
45 #  define SQR1(w) \
46     SQR_nibble((w) >> 28) << 24 | SQR_nibble((w) >> 24) << 16 | \
47     SQR_nibble((w) >> 20) <<  8 | SQR_nibble((w) >> 16)
48 #  define SQR0(w) \
49     SQR_nibble((w) >> 12) << 24 | SQR_nibble((w) >>  8) << 16 | \
50     SQR_nibble((w) >>  4) <<  8 | SQR_nibble((w)      )
51 # endif
52
53 # if !defined(OPENSSL_BN_ASM_GF2m)
54 /*
55  * Product of two polynomials a, b each with degree < BN_BITS2 - 1, result is
56  * a polynomial r with degree < 2 * BN_BITS - 1 The caller MUST ensure that
57  * the variables have the right amount of space allocated.
58  */
59 #  ifdef THIRTY_TWO_BIT
60 static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
61                             const BN_ULONG b)
62 {
63     register BN_ULONG h, l, s;
64     BN_ULONG tab[8], top2b = a >> 30;
65     register BN_ULONG a1, a2, a4;
66
67     a1 = a & (0x3FFFFFFF);
68     a2 = a1 << 1;
69     a4 = a2 << 1;
70
71     tab[0] = 0;
72     tab[1] = a1;
73     tab[2] = a2;
74     tab[3] = a1 ^ a2;
75     tab[4] = a4;
76     tab[5] = a1 ^ a4;
77     tab[6] = a2 ^ a4;
78     tab[7] = a1 ^ a2 ^ a4;
79
80     s = tab[b & 0x7];
81     l = s;
82     s = tab[b >> 3 & 0x7];
83     l ^= s << 3;
84     h = s >> 29;
85     s = tab[b >> 6 & 0x7];
86     l ^= s << 6;
87     h ^= s >> 26;
88     s = tab[b >> 9 & 0x7];
89     l ^= s << 9;
90     h ^= s >> 23;
91     s = tab[b >> 12 & 0x7];
92     l ^= s << 12;
93     h ^= s >> 20;
94     s = tab[b >> 15 & 0x7];
95     l ^= s << 15;
96     h ^= s >> 17;
97     s = tab[b >> 18 & 0x7];
98     l ^= s << 18;
99     h ^= s >> 14;
100     s = tab[b >> 21 & 0x7];
101     l ^= s << 21;
102     h ^= s >> 11;
103     s = tab[b >> 24 & 0x7];
104     l ^= s << 24;
105     h ^= s >> 8;
106     s = tab[b >> 27 & 0x7];
107     l ^= s << 27;
108     h ^= s >> 5;
109     s = tab[b >> 30];
110     l ^= s << 30;
111     h ^= s >> 2;
112
113     /* compensate for the top two bits of a */
114
115     if (top2b & 01) {
116         l ^= b << 30;
117         h ^= b >> 2;
118     }
119     if (top2b & 02) {
120         l ^= b << 31;
121         h ^= b >> 1;
122     }
123
124     *r1 = h;
125     *r0 = l;
126 }
127 #  endif
128 #  if defined(SIXTY_FOUR_BIT) || defined(SIXTY_FOUR_BIT_LONG)
129 static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
130                             const BN_ULONG b)
131 {
132     register BN_ULONG h, l, s;
133     BN_ULONG tab[16], top3b = a >> 61;
134     register BN_ULONG a1, a2, a4, a8;
135
136     a1 = a & (0x1FFFFFFFFFFFFFFFULL);
137     a2 = a1 << 1;
138     a4 = a2 << 1;
139     a8 = a4 << 1;
140
141     tab[0] = 0;
142     tab[1] = a1;
143     tab[2] = a2;
144     tab[3] = a1 ^ a2;
145     tab[4] = a4;
146     tab[5] = a1 ^ a4;
147     tab[6] = a2 ^ a4;
148     tab[7] = a1 ^ a2 ^ a4;
149     tab[8] = a8;
150     tab[9] = a1 ^ a8;
151     tab[10] = a2 ^ a8;
152     tab[11] = a1 ^ a2 ^ a8;
153     tab[12] = a4 ^ a8;
154     tab[13] = a1 ^ a4 ^ a8;
155     tab[14] = a2 ^ a4 ^ a8;
156     tab[15] = a1 ^ a2 ^ a4 ^ a8;
157
158     s = tab[b & 0xF];
159     l = s;
160     s = tab[b >> 4 & 0xF];
161     l ^= s << 4;
162     h = s >> 60;
163     s = tab[b >> 8 & 0xF];
164     l ^= s << 8;
165     h ^= s >> 56;
166     s = tab[b >> 12 & 0xF];
167     l ^= s << 12;
168     h ^= s >> 52;
169     s = tab[b >> 16 & 0xF];
170     l ^= s << 16;
171     h ^= s >> 48;
172     s = tab[b >> 20 & 0xF];
173     l ^= s << 20;
174     h ^= s >> 44;
175     s = tab[b >> 24 & 0xF];
176     l ^= s << 24;
177     h ^= s >> 40;
178     s = tab[b >> 28 & 0xF];
179     l ^= s << 28;
180     h ^= s >> 36;
181     s = tab[b >> 32 & 0xF];
182     l ^= s << 32;
183     h ^= s >> 32;
184     s = tab[b >> 36 & 0xF];
185     l ^= s << 36;
186     h ^= s >> 28;
187     s = tab[b >> 40 & 0xF];
188     l ^= s << 40;
189     h ^= s >> 24;
190     s = tab[b >> 44 & 0xF];
191     l ^= s << 44;
192     h ^= s >> 20;
193     s = tab[b >> 48 & 0xF];
194     l ^= s << 48;
195     h ^= s >> 16;
196     s = tab[b >> 52 & 0xF];
197     l ^= s << 52;
198     h ^= s >> 12;
199     s = tab[b >> 56 & 0xF];
200     l ^= s << 56;
201     h ^= s >> 8;
202     s = tab[b >> 60];
203     l ^= s << 60;
204     h ^= s >> 4;
205
206     /* compensate for the top three bits of a */
207
208     if (top3b & 01) {
209         l ^= b << 61;
210         h ^= b >> 3;
211     }
212     if (top3b & 02) {
213         l ^= b << 62;
214         h ^= b >> 2;
215     }
216     if (top3b & 04) {
217         l ^= b << 63;
218         h ^= b >> 1;
219     }
220
221     *r1 = h;
222     *r0 = l;
223 }
224 #  endif
225
226 /*
227  * Product of two polynomials a, b each with degree < 2 * BN_BITS2 - 1,
228  * result is a polynomial r with degree < 4 * BN_BITS2 - 1 The caller MUST
229  * ensure that the variables have the right amount of space allocated.
230  */
231 static void bn_GF2m_mul_2x2(BN_ULONG *r, const BN_ULONG a1, const BN_ULONG a0,
232                             const BN_ULONG b1, const BN_ULONG b0)
233 {
234     BN_ULONG m1, m0;
235     /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */
236     bn_GF2m_mul_1x1(r + 3, r + 2, a1, b1);
237     bn_GF2m_mul_1x1(r + 1, r, a0, b0);
238     bn_GF2m_mul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1);
239     /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */
240     r[2] ^= m1 ^ r[1] ^ r[3];   /* h0 ^= m1 ^ l1 ^ h1; */
241     r[1] = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0; /* l1 ^= l0 ^ h0 ^ m0; */
242 }
243 # else
244 void bn_GF2m_mul_2x2(BN_ULONG *r, BN_ULONG a1, BN_ULONG a0, BN_ULONG b1,
245                      BN_ULONG b0);
246 # endif
247
248 /*
249  * Add polynomials a and b and store result in r; r could be a or b, a and b
250  * could be equal; r is the bitwise XOR of a and b.
251  */
252 int BN_GF2m_add(BIGNUM *r, const BIGNUM *a, const BIGNUM *b)
253 {
254     int i;
255     const BIGNUM *at, *bt;
256
257     bn_check_top(a);
258     bn_check_top(b);
259
260     if (a->top < b->top) {
261         at = b;
262         bt = a;
263     } else {
264         at = a;
265         bt = b;
266     }
267
268     if (bn_wexpand(r, at->top) == NULL)
269         return 0;
270
271     for (i = 0; i < bt->top; i++) {
272         r->d[i] = at->d[i] ^ bt->d[i];
273     }
274     for (; i < at->top; i++) {
275         r->d[i] = at->d[i];
276     }
277
278     r->top = at->top;
279     bn_correct_top(r);
280
281     return 1;
282 }
283
284 /*-
285  * Some functions allow for representation of the irreducible polynomials
286  * as an int[], say p.  The irreducible f(t) is then of the form:
287  *     t^p[0] + t^p[1] + ... + t^p[k]
288  * where m = p[0] > p[1] > ... > p[k] = 0.
289  */
290
291 /* Performs modular reduction of a and store result in r.  r could be a. */
292 int BN_GF2m_mod_arr(BIGNUM *r, const BIGNUM *a, const int p[])
293 {
294     int j, k;
295     int n, dN, d0, d1;
296     BN_ULONG zz, *z;
297
298     bn_check_top(a);
299
300     if (!p[0]) {
301         /* reduction mod 1 => return 0 */
302         BN_zero(r);
303         return 1;
304     }
305
306     /*
307      * Since the algorithm does reduction in the r value, if a != r, copy the
308      * contents of a into r so we can do reduction in r.
309      */
310     if (a != r) {
311         if (!bn_wexpand(r, a->top))
312             return 0;
313         for (j = 0; j < a->top; j++) {
314             r->d[j] = a->d[j];
315         }
316         r->top = a->top;
317     }
318     z = r->d;
319
320     /* start reduction */
321     dN = p[0] / BN_BITS2;
322     for (j = r->top - 1; j > dN;) {
323         zz = z[j];
324         if (z[j] == 0) {
325             j--;
326             continue;
327         }
328         z[j] = 0;
329
330         for (k = 1; p[k] != 0; k++) {
331             /* reducing component t^p[k] */
332             n = p[0] - p[k];
333             d0 = n % BN_BITS2;
334             d1 = BN_BITS2 - d0;
335             n /= BN_BITS2;
336             z[j - n] ^= (zz >> d0);
337             if (d0)
338                 z[j - n - 1] ^= (zz << d1);
339         }
340
341         /* reducing component t^0 */
342         n = dN;
343         d0 = p[0] % BN_BITS2;
344         d1 = BN_BITS2 - d0;
345         z[j - n] ^= (zz >> d0);
346         if (d0)
347             z[j - n - 1] ^= (zz << d1);
348     }
349
350     /* final round of reduction */
351     while (j == dN) {
352
353         d0 = p[0] % BN_BITS2;
354         zz = z[dN] >> d0;
355         if (zz == 0)
356             break;
357         d1 = BN_BITS2 - d0;
358
359         /* clear up the top d1 bits */
360         if (d0)
361             z[dN] = (z[dN] << d1) >> d1;
362         else
363             z[dN] = 0;
364         z[0] ^= zz;             /* reduction t^0 component */
365
366         for (k = 1; p[k] != 0; k++) {
367             BN_ULONG tmp_ulong;
368
369             /* reducing component t^p[k] */
370             n = p[k] / BN_BITS2;
371             d0 = p[k] % BN_BITS2;
372             d1 = BN_BITS2 - d0;
373             z[n] ^= (zz << d0);
374             if (d0 && (tmp_ulong = zz >> d1))
375                 z[n + 1] ^= tmp_ulong;
376         }
377
378     }
379
380     bn_correct_top(r);
381     return 1;
382 }
383
384 /*
385  * Performs modular reduction of a by p and store result in r.  r could be a.
386  * This function calls down to the BN_GF2m_mod_arr implementation; this wrapper
387  * function is only provided for convenience; for best performance, use the
388  * BN_GF2m_mod_arr function.
389  */
390 int BN_GF2m_mod(BIGNUM *r, const BIGNUM *a, const BIGNUM *p)
391 {
392     int ret = 0;
393     int arr[6];
394     bn_check_top(a);
395     bn_check_top(p);
396     ret = BN_GF2m_poly2arr(p, arr, OSSL_NELEM(arr));
397     if (!ret || ret > (int)OSSL_NELEM(arr)) {
398         BNerr(BN_F_BN_GF2M_MOD, BN_R_INVALID_LENGTH);
399         return 0;
400     }
401     ret = BN_GF2m_mod_arr(r, a, arr);
402     bn_check_top(r);
403     return ret;
404 }
405
406 /*
407  * Compute the product of two polynomials a and b, reduce modulo p, and store
408  * the result in r.  r could be a or b; a could be b.
409  */
410 int BN_GF2m_mod_mul_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
411                         const int p[], BN_CTX *ctx)
412 {
413     int zlen, i, j, k, ret = 0;
414     BIGNUM *s;
415     BN_ULONG x1, x0, y1, y0, zz[4];
416
417     bn_check_top(a);
418     bn_check_top(b);
419
420     if (a == b) {
421         return BN_GF2m_mod_sqr_arr(r, a, p, ctx);
422     }
423
424     BN_CTX_start(ctx);
425     if ((s = BN_CTX_get(ctx)) == NULL)
426         goto err;
427
428     zlen = a->top + b->top + 4;
429     if (!bn_wexpand(s, zlen))
430         goto err;
431     s->top = zlen;
432
433     for (i = 0; i < zlen; i++)
434         s->d[i] = 0;
435
436     for (j = 0; j < b->top; j += 2) {
437         y0 = b->d[j];
438         y1 = ((j + 1) == b->top) ? 0 : b->d[j + 1];
439         for (i = 0; i < a->top; i += 2) {
440             x0 = a->d[i];
441             x1 = ((i + 1) == a->top) ? 0 : a->d[i + 1];
442             bn_GF2m_mul_2x2(zz, x1, x0, y1, y0);
443             for (k = 0; k < 4; k++)
444                 s->d[i + j + k] ^= zz[k];
445         }
446     }
447
448     bn_correct_top(s);
449     if (BN_GF2m_mod_arr(r, s, p))
450         ret = 1;
451     bn_check_top(r);
452
453  err:
454     BN_CTX_end(ctx);
455     return ret;
456 }
457
458 /*
459  * Compute the product of two polynomials a and b, reduce modulo p, and store
460  * the result in r.  r could be a or b; a could equal b. This function calls
461  * down to the BN_GF2m_mod_mul_arr implementation; this wrapper function is
462  * only provided for convenience; for best performance, use the
463  * BN_GF2m_mod_mul_arr function.
464  */
465 int BN_GF2m_mod_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
466                     const BIGNUM *p, BN_CTX *ctx)
467 {
468     int ret = 0;
469     const int max = BN_num_bits(p) + 1;
470     int *arr = NULL;
471     bn_check_top(a);
472     bn_check_top(b);
473     bn_check_top(p);
474     if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
475         goto err;
476     ret = BN_GF2m_poly2arr(p, arr, max);
477     if (!ret || ret > max) {
478         BNerr(BN_F_BN_GF2M_MOD_MUL, BN_R_INVALID_LENGTH);
479         goto err;
480     }
481     ret = BN_GF2m_mod_mul_arr(r, a, b, arr, ctx);
482     bn_check_top(r);
483  err:
484     OPENSSL_free(arr);
485     return ret;
486 }
487
488 /* Square a, reduce the result mod p, and store it in a.  r could be a. */
489 int BN_GF2m_mod_sqr_arr(BIGNUM *r, const BIGNUM *a, const int p[],
490                         BN_CTX *ctx)
491 {
492     int i, ret = 0;
493     BIGNUM *s;
494
495     bn_check_top(a);
496     BN_CTX_start(ctx);
497     if ((s = BN_CTX_get(ctx)) == NULL)
498         goto err;
499     if (!bn_wexpand(s, 2 * a->top))
500         goto err;
501
502     for (i = a->top - 1; i >= 0; i--) {
503         s->d[2 * i + 1] = SQR1(a->d[i]);
504         s->d[2 * i] = SQR0(a->d[i]);
505     }
506
507     s->top = 2 * a->top;
508     bn_correct_top(s);
509     if (!BN_GF2m_mod_arr(r, s, p))
510         goto err;
511     bn_check_top(r);
512     ret = 1;
513  err:
514     BN_CTX_end(ctx);
515     return ret;
516 }
517
518 /*
519  * Square a, reduce the result mod p, and store it in a.  r could be a. This
520  * function calls down to the BN_GF2m_mod_sqr_arr implementation; this
521  * wrapper function is only provided for convenience; for best performance,
522  * use the BN_GF2m_mod_sqr_arr function.
523  */
524 int BN_GF2m_mod_sqr(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
525 {
526     int ret = 0;
527     const int max = BN_num_bits(p) + 1;
528     int *arr = NULL;
529
530     bn_check_top(a);
531     bn_check_top(p);
532     if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
533         goto err;
534     ret = BN_GF2m_poly2arr(p, arr, max);
535     if (!ret || ret > max) {
536         BNerr(BN_F_BN_GF2M_MOD_SQR, BN_R_INVALID_LENGTH);
537         goto err;
538     }
539     ret = BN_GF2m_mod_sqr_arr(r, a, arr, ctx);
540     bn_check_top(r);
541  err:
542     OPENSSL_free(arr);
543     return ret;
544 }
545
546 /*
547  * Invert a, reduce modulo p, and store the result in r. r could be a. Uses
548  * Modified Almost Inverse Algorithm (Algorithm 10) from Hankerson, D.,
549  * Hernandez, J.L., and Menezes, A.  "Software Implementation of Elliptic
550  * Curve Cryptography Over Binary Fields".
551  */
552 static int BN_GF2m_mod_inv_vartime(BIGNUM *r, const BIGNUM *a,
553                                    const BIGNUM *p, BN_CTX *ctx)
554 {
555     BIGNUM *b, *c = NULL, *u = NULL, *v = NULL, *tmp;
556     int ret = 0;
557
558     bn_check_top(a);
559     bn_check_top(p);
560
561     BN_CTX_start(ctx);
562
563     b = BN_CTX_get(ctx);
564     c = BN_CTX_get(ctx);
565     u = BN_CTX_get(ctx);
566     v = BN_CTX_get(ctx);
567     if (v == NULL)
568         goto err;
569
570     if (!BN_GF2m_mod(u, a, p))
571         goto err;
572     if (BN_is_zero(u))
573         goto err;
574
575     if (!BN_copy(v, p))
576         goto err;
577 # if 0
578     if (!BN_one(b))
579         goto err;
580
581     while (1) {
582         while (!BN_is_odd(u)) {
583             if (BN_is_zero(u))
584                 goto err;
585             if (!BN_rshift1(u, u))
586                 goto err;
587             if (BN_is_odd(b)) {
588                 if (!BN_GF2m_add(b, b, p))
589                     goto err;
590             }
591             if (!BN_rshift1(b, b))
592                 goto err;
593         }
594
595         if (BN_abs_is_word(u, 1))
596             break;
597
598         if (BN_num_bits(u) < BN_num_bits(v)) {
599             tmp = u;
600             u = v;
601             v = tmp;
602             tmp = b;
603             b = c;
604             c = tmp;
605         }
606
607         if (!BN_GF2m_add(u, u, v))
608             goto err;
609         if (!BN_GF2m_add(b, b, c))
610             goto err;
611     }
612 # else
613     {
614         int i;
615         int ubits = BN_num_bits(u);
616         int vbits = BN_num_bits(v); /* v is copy of p */
617         int top = p->top;
618         BN_ULONG *udp, *bdp, *vdp, *cdp;
619
620         if (!bn_wexpand(u, top))
621             goto err;
622         udp = u->d;
623         for (i = u->top; i < top; i++)
624             udp[i] = 0;
625         u->top = top;
626         if (!bn_wexpand(b, top))
627           goto err;
628         bdp = b->d;
629         bdp[0] = 1;
630         for (i = 1; i < top; i++)
631             bdp[i] = 0;
632         b->top = top;
633         if (!bn_wexpand(c, top))
634           goto err;
635         cdp = c->d;
636         for (i = 0; i < top; i++)
637             cdp[i] = 0;
638         c->top = top;
639         vdp = v->d;             /* It pays off to "cache" *->d pointers,
640                                  * because it allows optimizer to be more
641                                  * aggressive. But we don't have to "cache"
642                                  * p->d, because *p is declared 'const'... */
643         while (1) {
644             while (ubits && !(udp[0] & 1)) {
645                 BN_ULONG u0, u1, b0, b1, mask;
646
647                 u0 = udp[0];
648                 b0 = bdp[0];
649                 mask = (BN_ULONG)0 - (b0 & 1);
650                 b0 ^= p->d[0] & mask;
651                 for (i = 0; i < top - 1; i++) {
652                     u1 = udp[i + 1];
653                     udp[i] = ((u0 >> 1) | (u1 << (BN_BITS2 - 1))) & BN_MASK2;
654                     u0 = u1;
655                     b1 = bdp[i + 1] ^ (p->d[i + 1] & mask);
656                     bdp[i] = ((b0 >> 1) | (b1 << (BN_BITS2 - 1))) & BN_MASK2;
657                     b0 = b1;
658                 }
659                 udp[i] = u0 >> 1;
660                 bdp[i] = b0 >> 1;
661                 ubits--;
662             }
663
664             if (ubits <= BN_BITS2) {
665                 if (udp[0] == 0) /* poly was reducible */
666                     goto err;
667                 if (udp[0] == 1)
668                     break;
669             }
670
671             if (ubits < vbits) {
672                 i = ubits;
673                 ubits = vbits;
674                 vbits = i;
675                 tmp = u;
676                 u = v;
677                 v = tmp;
678                 tmp = b;
679                 b = c;
680                 c = tmp;
681                 udp = vdp;
682                 vdp = v->d;
683                 bdp = cdp;
684                 cdp = c->d;
685             }
686             for (i = 0; i < top; i++) {
687                 udp[i] ^= vdp[i];
688                 bdp[i] ^= cdp[i];
689             }
690             if (ubits == vbits) {
691                 BN_ULONG ul;
692                 int utop = (ubits - 1) / BN_BITS2;
693
694                 while ((ul = udp[utop]) == 0 && utop)
695                     utop--;
696                 ubits = utop * BN_BITS2 + BN_num_bits_word(ul);
697             }
698         }
699         bn_correct_top(b);
700     }
701 # endif
702
703     if (!BN_copy(r, b))
704         goto err;
705     bn_check_top(r);
706     ret = 1;
707
708  err:
709 # ifdef BN_DEBUG                /* BN_CTX_end would complain about the
710                                  * expanded form */
711     bn_correct_top(c);
712     bn_correct_top(u);
713     bn_correct_top(v);
714 # endif
715     BN_CTX_end(ctx);
716     return ret;
717 }
718
719 /*-
720  * Wrapper for BN_GF2m_mod_inv_vartime that blinds the input before calling.
721  * This is not constant time.
722  * But it does eliminate first order deduction on the input.
723  */
724 int BN_GF2m_mod_inv(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
725 {
726     BIGNUM *b = NULL;
727     int ret = 0;
728
729     BN_CTX_start(ctx);
730     if ((b = BN_CTX_get(ctx)) == NULL)
731         goto err;
732
733     /* generate blinding value */
734     do {
735         if (!BN_priv_rand(b, BN_num_bits(p) - 1,
736                           BN_RAND_TOP_ANY, BN_RAND_BOTTOM_ANY))
737             goto err;
738     } while (BN_is_zero(b));
739
740     /* r := a * b */
741     if (!BN_GF2m_mod_mul(r, a, b, p, ctx))
742         goto err;
743
744     /* r := 1/(a * b) */
745     if (!BN_GF2m_mod_inv_vartime(r, r, p, ctx))
746         goto err;
747
748     /* r := b/(a * b) = 1/a */
749     if (!BN_GF2m_mod_mul(r, r, b, p, ctx))
750         goto err;
751
752     ret = 1;
753
754  err:
755     BN_CTX_end(ctx);
756     return ret;
757 }
758
759 /*
760  * Invert xx, reduce modulo p, and store the result in r. r could be xx.
761  * This function calls down to the BN_GF2m_mod_inv implementation; this
762  * wrapper function is only provided for convenience; for best performance,
763  * use the BN_GF2m_mod_inv function.
764  */
765 int BN_GF2m_mod_inv_arr(BIGNUM *r, const BIGNUM *xx, const int p[],
766                         BN_CTX *ctx)
767 {
768     BIGNUM *field;
769     int ret = 0;
770
771     bn_check_top(xx);
772     BN_CTX_start(ctx);
773     if ((field = BN_CTX_get(ctx)) == NULL)
774         goto err;
775     if (!BN_GF2m_arr2poly(p, field))
776         goto err;
777
778     ret = BN_GF2m_mod_inv(r, xx, field, ctx);
779     bn_check_top(r);
780
781  err:
782     BN_CTX_end(ctx);
783     return ret;
784 }
785
786 /*
787  * Divide y by x, reduce modulo p, and store the result in r. r could be x
788  * or y, x could equal y.
789  */
790 int BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x,
791                     const BIGNUM *p, BN_CTX *ctx)
792 {
793     BIGNUM *xinv = NULL;
794     int ret = 0;
795
796     bn_check_top(y);
797     bn_check_top(x);
798     bn_check_top(p);
799
800     BN_CTX_start(ctx);
801     xinv = BN_CTX_get(ctx);
802     if (xinv == NULL)
803         goto err;
804
805     if (!BN_GF2m_mod_inv(xinv, x, p, ctx))
806         goto err;
807     if (!BN_GF2m_mod_mul(r, y, xinv, p, ctx))
808         goto err;
809     bn_check_top(r);
810     ret = 1;
811
812  err:
813     BN_CTX_end(ctx);
814     return ret;
815 }
816
817 /*
818  * Divide yy by xx, reduce modulo p, and store the result in r. r could be xx
819  * * or yy, xx could equal yy. This function calls down to the
820  * BN_GF2m_mod_div implementation; this wrapper function is only provided for
821  * convenience; for best performance, use the BN_GF2m_mod_div function.
822  */
823 int BN_GF2m_mod_div_arr(BIGNUM *r, const BIGNUM *yy, const BIGNUM *xx,
824                         const int p[], BN_CTX *ctx)
825 {
826     BIGNUM *field;
827     int ret = 0;
828
829     bn_check_top(yy);
830     bn_check_top(xx);
831
832     BN_CTX_start(ctx);
833     if ((field = BN_CTX_get(ctx)) == NULL)
834         goto err;
835     if (!BN_GF2m_arr2poly(p, field))
836         goto err;
837
838     ret = BN_GF2m_mod_div(r, yy, xx, field, ctx);
839     bn_check_top(r);
840
841  err:
842     BN_CTX_end(ctx);
843     return ret;
844 }
845
846 /*
847  * Compute the bth power of a, reduce modulo p, and store the result in r.  r
848  * could be a. Uses simple square-and-multiply algorithm A.5.1 from IEEE
849  * P1363.
850  */
851 int BN_GF2m_mod_exp_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
852                         const int p[], BN_CTX *ctx)
853 {
854     int ret = 0, i, n;
855     BIGNUM *u;
856
857     bn_check_top(a);
858     bn_check_top(b);
859
860     if (BN_is_zero(b))
861         return BN_one(r);
862
863     if (BN_abs_is_word(b, 1))
864         return (BN_copy(r, a) != NULL);
865
866     BN_CTX_start(ctx);
867     if ((u = BN_CTX_get(ctx)) == NULL)
868         goto err;
869
870     if (!BN_GF2m_mod_arr(u, a, p))
871         goto err;
872
873     n = BN_num_bits(b) - 1;
874     for (i = n - 1; i >= 0; i--) {
875         if (!BN_GF2m_mod_sqr_arr(u, u, p, ctx))
876             goto err;
877         if (BN_is_bit_set(b, i)) {
878             if (!BN_GF2m_mod_mul_arr(u, u, a, p, ctx))
879                 goto err;
880         }
881     }
882     if (!BN_copy(r, u))
883         goto err;
884     bn_check_top(r);
885     ret = 1;
886  err:
887     BN_CTX_end(ctx);
888     return ret;
889 }
890
891 /*
892  * Compute the bth power of a, reduce modulo p, and store the result in r.  r
893  * could be a. This function calls down to the BN_GF2m_mod_exp_arr
894  * implementation; this wrapper function is only provided for convenience;
895  * for best performance, use the BN_GF2m_mod_exp_arr function.
896  */
897 int BN_GF2m_mod_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
898                     const BIGNUM *p, BN_CTX *ctx)
899 {
900     int ret = 0;
901     const int max = BN_num_bits(p) + 1;
902     int *arr = NULL;
903     bn_check_top(a);
904     bn_check_top(b);
905     bn_check_top(p);
906     if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
907         goto err;
908     ret = BN_GF2m_poly2arr(p, arr, max);
909     if (!ret || ret > max) {
910         BNerr(BN_F_BN_GF2M_MOD_EXP, BN_R_INVALID_LENGTH);
911         goto err;
912     }
913     ret = BN_GF2m_mod_exp_arr(r, a, b, arr, ctx);
914     bn_check_top(r);
915  err:
916     OPENSSL_free(arr);
917     return ret;
918 }
919
920 /*
921  * Compute the square root of a, reduce modulo p, and store the result in r.
922  * r could be a. Uses exponentiation as in algorithm A.4.1 from IEEE P1363.
923  */
924 int BN_GF2m_mod_sqrt_arr(BIGNUM *r, const BIGNUM *a, const int p[],
925                          BN_CTX *ctx)
926 {
927     int ret = 0;
928     BIGNUM *u;
929
930     bn_check_top(a);
931
932     if (!p[0]) {
933         /* reduction mod 1 => return 0 */
934         BN_zero(r);
935         return 1;
936     }
937
938     BN_CTX_start(ctx);
939     if ((u = BN_CTX_get(ctx)) == NULL)
940         goto err;
941
942     if (!BN_set_bit(u, p[0] - 1))
943         goto err;
944     ret = BN_GF2m_mod_exp_arr(r, a, u, p, ctx);
945     bn_check_top(r);
946
947  err:
948     BN_CTX_end(ctx);
949     return ret;
950 }
951
952 /*
953  * Compute the square root of a, reduce modulo p, and store the result in r.
954  * r could be a. This function calls down to the BN_GF2m_mod_sqrt_arr
955  * implementation; this wrapper function is only provided for convenience;
956  * for best performance, use the BN_GF2m_mod_sqrt_arr function.
957  */
958 int BN_GF2m_mod_sqrt(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
959 {
960     int ret = 0;
961     const int max = BN_num_bits(p) + 1;
962     int *arr = NULL;
963     bn_check_top(a);
964     bn_check_top(p);
965     if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
966         goto err;
967     ret = BN_GF2m_poly2arr(p, arr, max);
968     if (!ret || ret > max) {
969         BNerr(BN_F_BN_GF2M_MOD_SQRT, BN_R_INVALID_LENGTH);
970         goto err;
971     }
972     ret = BN_GF2m_mod_sqrt_arr(r, a, arr, ctx);
973     bn_check_top(r);
974  err:
975     OPENSSL_free(arr);
976     return ret;
977 }
978
979 /*
980  * Find r such that r^2 + r = a mod p.  r could be a. If no r exists returns
981  * 0. Uses algorithms A.4.7 and A.4.6 from IEEE P1363.
982  */
983 int BN_GF2m_mod_solve_quad_arr(BIGNUM *r, const BIGNUM *a_, const int p[],
984                                BN_CTX *ctx)
985 {
986     int ret = 0, count = 0, j;
987     BIGNUM *a, *z, *rho, *w, *w2, *tmp;
988
989     bn_check_top(a_);
990
991     if (!p[0]) {
992         /* reduction mod 1 => return 0 */
993         BN_zero(r);
994         return 1;
995     }
996
997     BN_CTX_start(ctx);
998     a = BN_CTX_get(ctx);
999     z = BN_CTX_get(ctx);
1000     w = BN_CTX_get(ctx);
1001     if (w == NULL)
1002         goto err;
1003
1004     if (!BN_GF2m_mod_arr(a, a_, p))
1005         goto err;
1006
1007     if (BN_is_zero(a)) {
1008         BN_zero(r);
1009         ret = 1;
1010         goto err;
1011     }
1012
1013     if (p[0] & 0x1) {           /* m is odd */
1014         /* compute half-trace of a */
1015         if (!BN_copy(z, a))
1016             goto err;
1017         for (j = 1; j <= (p[0] - 1) / 2; j++) {
1018             if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1019                 goto err;
1020             if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1021                 goto err;
1022             if (!BN_GF2m_add(z, z, a))
1023                 goto err;
1024         }
1025
1026     } else {                    /* m is even */
1027
1028         rho = BN_CTX_get(ctx);
1029         w2 = BN_CTX_get(ctx);
1030         tmp = BN_CTX_get(ctx);
1031         if (tmp == NULL)
1032             goto err;
1033         do {
1034             if (!BN_priv_rand(rho, p[0], BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ANY))
1035                 goto err;
1036             if (!BN_GF2m_mod_arr(rho, rho, p))
1037                 goto err;
1038             BN_zero(z);
1039             if (!BN_copy(w, rho))
1040                 goto err;
1041             for (j = 1; j <= p[0] - 1; j++) {
1042                 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1043                     goto err;
1044                 if (!BN_GF2m_mod_sqr_arr(w2, w, p, ctx))
1045                     goto err;
1046                 if (!BN_GF2m_mod_mul_arr(tmp, w2, a, p, ctx))
1047                     goto err;
1048                 if (!BN_GF2m_add(z, z, tmp))
1049                     goto err;
1050                 if (!BN_GF2m_add(w, w2, rho))
1051                     goto err;
1052             }
1053             count++;
1054         } while (BN_is_zero(w) && (count < MAX_ITERATIONS));
1055         if (BN_is_zero(w)) {
1056             BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD_ARR, BN_R_TOO_MANY_ITERATIONS);
1057             goto err;
1058         }
1059     }
1060
1061     if (!BN_GF2m_mod_sqr_arr(w, z, p, ctx))
1062         goto err;
1063     if (!BN_GF2m_add(w, z, w))
1064         goto err;
1065     if (BN_GF2m_cmp(w, a)) {
1066         BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD_ARR, BN_R_NO_SOLUTION);
1067         goto err;
1068     }
1069
1070     if (!BN_copy(r, z))
1071         goto err;
1072     bn_check_top(r);
1073
1074     ret = 1;
1075
1076  err:
1077     BN_CTX_end(ctx);
1078     return ret;
1079 }
1080
1081 /*
1082  * Find r such that r^2 + r = a mod p.  r could be a. If no r exists returns
1083  * 0. This function calls down to the BN_GF2m_mod_solve_quad_arr
1084  * implementation; this wrapper function is only provided for convenience;
1085  * for best performance, use the BN_GF2m_mod_solve_quad_arr function.
1086  */
1087 int BN_GF2m_mod_solve_quad(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
1088                            BN_CTX *ctx)
1089 {
1090     int ret = 0;
1091     const int max = BN_num_bits(p) + 1;
1092     int *arr = NULL;
1093     bn_check_top(a);
1094     bn_check_top(p);
1095     if ((arr = OPENSSL_malloc(sizeof(*arr) * max)) == NULL)
1096         goto err;
1097     ret = BN_GF2m_poly2arr(p, arr, max);
1098     if (!ret || ret > max) {
1099         BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD, BN_R_INVALID_LENGTH);
1100         goto err;
1101     }
1102     ret = BN_GF2m_mod_solve_quad_arr(r, a, arr, ctx);
1103     bn_check_top(r);
1104  err:
1105     OPENSSL_free(arr);
1106     return ret;
1107 }
1108
1109 /*
1110  * Convert the bit-string representation of a polynomial ( \sum_{i=0}^n a_i *
1111  * x^i) into an array of integers corresponding to the bits with non-zero
1112  * coefficient.  Array is terminated with -1. Up to max elements of the array
1113  * will be filled.  Return value is total number of array elements that would
1114  * be filled if array was large enough.
1115  */
1116 int BN_GF2m_poly2arr(const BIGNUM *a, int p[], int max)
1117 {
1118     int i, j, k = 0;
1119     BN_ULONG mask;
1120
1121     if (BN_is_zero(a))
1122         return 0;
1123
1124     for (i = a->top - 1; i >= 0; i--) {
1125         if (!a->d[i])
1126             /* skip word if a->d[i] == 0 */
1127             continue;
1128         mask = BN_TBIT;
1129         for (j = BN_BITS2 - 1; j >= 0; j--) {
1130             if (a->d[i] & mask) {
1131                 if (k < max)
1132                     p[k] = BN_BITS2 * i + j;
1133                 k++;
1134             }
1135             mask >>= 1;
1136         }
1137     }
1138
1139     if (k < max) {
1140         p[k] = -1;
1141         k++;
1142     }
1143
1144     return k;
1145 }
1146
1147 /*
1148  * Convert the coefficient array representation of a polynomial to a
1149  * bit-string.  The array must be terminated by -1.
1150  */
1151 int BN_GF2m_arr2poly(const int p[], BIGNUM *a)
1152 {
1153     int i;
1154
1155     bn_check_top(a);
1156     BN_zero(a);
1157     for (i = 0; p[i] != -1; i++) {
1158         if (BN_set_bit(a, p[i]) == 0)
1159             return 0;
1160     }
1161     bn_check_top(a);
1162
1163     return 1;
1164 }
1165
1166 #endif