1/* crypto/bn/bn_gf2m.c */
2/* ====================================================================
3 * Copyright 2002 Sun Microsystems, Inc. ALL RIGHTS RESERVED.
4 *
5 * The Elliptic Curve Public-Key Crypto Library (ECC Code) included
6 * herein is developed by SUN MICROSYSTEMS, INC., and is contributed
7 * to the OpenSSL project.
8 *
9 * The ECC Code is licensed pursuant to the OpenSSL open source
10 * license provided below.
11 *
12 * In addition, Sun covenants to all licensees who provide a reciprocal
13 * covenant with respect to their own patents if any, not to sue under
14 * current and future patent claims necessarily infringed by the making,
15 * using, practicing, selling, offering for sale and/or otherwise
16 * disposing of the ECC Code as delivered hereunder (or portions thereof),
17 * provided that such covenant shall not apply:
18 *  1) for code that a licensee deletes from the ECC Code;
19 *  2) separates from the ECC Code; or
20 *  3) for infringements caused by:
21 *       i) the modification of the ECC Code or
22 *      ii) the combination of the ECC Code with other software or
23 *          devices where such combination causes the infringement.
24 *
25 * The software is originally written by Sheueling Chang Shantz and
26 * Douglas Stebila of Sun Microsystems Laboratories.
27 *
28 */
29
30/*
31 * NOTE: This file is licensed pursuant to the OpenSSL license below and may
32 * be modified; but after modifications, the above covenant may no longer
33 * apply! In such cases, the corresponding paragraph ["In addition, Sun
34 * covenants ... causes the infringement."] and this note can be edited out;
35 * but please keep the Sun copyright notice and attribution.
36 */
37
38/* ====================================================================
39 * Copyright (c) 1998-2002 The OpenSSL Project.  All rights reserved.
40 *
41 * Redistribution and use in source and binary forms, with or without
42 * modification, are permitted provided that the following conditions
43 * are met:
44 *
45 * 1. Redistributions of source code must retain the above copyright
46 *    notice, this list of conditions and the following disclaimer.
47 *
48 * 2. Redistributions in binary form must reproduce the above copyright
49 *    notice, this list of conditions and the following disclaimer in
50 *    the documentation and/or other materials provided with the
51 *    distribution.
52 *
53 * 3. All advertising materials mentioning features or use of this
54 *    software must display the following acknowledgment:
55 *    "This product includes software developed by the OpenSSL Project
56 *    for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
57 *
58 * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
59 *    endorse or promote products derived from this software without
60 *    prior written permission. For written permission, please contact
61 *    openssl-core@openssl.org.
62 *
63 * 5. Products derived from this software may not be called "OpenSSL"
64 *    nor may "OpenSSL" appear in their names without prior written
65 *    permission of the OpenSSL Project.
66 *
67 * 6. Redistributions of any form whatsoever must retain the following
68 *    acknowledgment:
69 *    "This product includes software developed by the OpenSSL Project
70 *    for use in the OpenSSL Toolkit (http://www.openssl.org/)"
71 *
72 * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
73 * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
74 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
75 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE OpenSSL PROJECT OR
76 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
77 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
78 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
79 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
80 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
81 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
82 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
83 * OF THE POSSIBILITY OF SUCH DAMAGE.
84 * ====================================================================
85 *
86 * This product includes cryptographic software written by Eric Young
87 * (eay@cryptsoft.com).  This product includes software written by Tim
88 * Hudson (tjh@cryptsoft.com).
89 *
90 */
91
92#include <assert.h>
93#include <limits.h>
94#include <stdio.h>
95#include "cryptlib.h"
96#include "bn_lcl.h"
97
98#ifndef OPENSSL_NO_EC2M
99
100/*
101 * Maximum number of iterations before BN_GF2m_mod_solve_quad_arr should
102 * fail.
103 */
104# define MAX_ITERATIONS 50
105
106static const BN_ULONG SQR_tb[16] = { 0, 1, 4, 5, 16, 17, 20, 21,
107    64, 65, 68, 69, 80, 81, 84, 85
108};
109
110/* Platform-specific macros to accelerate squaring. */
111# if defined(SIXTY_FOUR_BIT) || defined(SIXTY_FOUR_BIT_LONG)
112#  define SQR1(w) \
113    SQR_tb[(w) >> 60 & 0xF] << 56 | SQR_tb[(w) >> 56 & 0xF] << 48 | \
114    SQR_tb[(w) >> 52 & 0xF] << 40 | SQR_tb[(w) >> 48 & 0xF] << 32 | \
115    SQR_tb[(w) >> 44 & 0xF] << 24 | SQR_tb[(w) >> 40 & 0xF] << 16 | \
116    SQR_tb[(w) >> 36 & 0xF] <<  8 | SQR_tb[(w) >> 32 & 0xF]
117#  define SQR0(w) \
118    SQR_tb[(w) >> 28 & 0xF] << 56 | SQR_tb[(w) >> 24 & 0xF] << 48 | \
119    SQR_tb[(w) >> 20 & 0xF] << 40 | SQR_tb[(w) >> 16 & 0xF] << 32 | \
120    SQR_tb[(w) >> 12 & 0xF] << 24 | SQR_tb[(w) >>  8 & 0xF] << 16 | \
121    SQR_tb[(w) >>  4 & 0xF] <<  8 | SQR_tb[(w)       & 0xF]
122# endif
123# ifdef THIRTY_TWO_BIT
124#  define SQR1(w) \
125    SQR_tb[(w) >> 28 & 0xF] << 24 | SQR_tb[(w) >> 24 & 0xF] << 16 | \
126    SQR_tb[(w) >> 20 & 0xF] <<  8 | SQR_tb[(w) >> 16 & 0xF]
127#  define SQR0(w) \
128    SQR_tb[(w) >> 12 & 0xF] << 24 | SQR_tb[(w) >>  8 & 0xF] << 16 | \
129    SQR_tb[(w) >>  4 & 0xF] <<  8 | SQR_tb[(w)       & 0xF]
130# endif
131
132# if !defined(OPENSSL_BN_ASM_GF2m)
133/*
134 * Product of two polynomials a, b each with degree < BN_BITS2 - 1, result is
135 * a polynomial r with degree < 2 * BN_BITS - 1 The caller MUST ensure that
136 * the variables have the right amount of space allocated.
137 */
138#  ifdef THIRTY_TWO_BIT
139static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
140                            const BN_ULONG b)
141{
142    register BN_ULONG h, l, s;
143    BN_ULONG tab[8], top2b = a >> 30;
144    register BN_ULONG a1, a2, a4;
145
146    a1 = a & (0x3FFFFFFF);
147    a2 = a1 << 1;
148    a4 = a2 << 1;
149
150    tab[0] = 0;
151    tab[1] = a1;
152    tab[2] = a2;
153    tab[3] = a1 ^ a2;
154    tab[4] = a4;
155    tab[5] = a1 ^ a4;
156    tab[6] = a2 ^ a4;
157    tab[7] = a1 ^ a2 ^ a4;
158
159    s = tab[b & 0x7];
160    l = s;
161    s = tab[b >> 3 & 0x7];
162    l ^= s << 3;
163    h = s >> 29;
164    s = tab[b >> 6 & 0x7];
165    l ^= s << 6;
166    h ^= s >> 26;
167    s = tab[b >> 9 & 0x7];
168    l ^= s << 9;
169    h ^= s >> 23;
170    s = tab[b >> 12 & 0x7];
171    l ^= s << 12;
172    h ^= s >> 20;
173    s = tab[b >> 15 & 0x7];
174    l ^= s << 15;
175    h ^= s >> 17;
176    s = tab[b >> 18 & 0x7];
177    l ^= s << 18;
178    h ^= s >> 14;
179    s = tab[b >> 21 & 0x7];
180    l ^= s << 21;
181    h ^= s >> 11;
182    s = tab[b >> 24 & 0x7];
183    l ^= s << 24;
184    h ^= s >> 8;
185    s = tab[b >> 27 & 0x7];
186    l ^= s << 27;
187    h ^= s >> 5;
188    s = tab[b >> 30];
189    l ^= s << 30;
190    h ^= s >> 2;
191
192    /* compensate for the top two bits of a */
193
194    if (top2b & 01) {
195        l ^= b << 30;
196        h ^= b >> 2;
197    }
198    if (top2b & 02) {
199        l ^= b << 31;
200        h ^= b >> 1;
201    }
202
203    *r1 = h;
204    *r0 = l;
205}
206#  endif
207#  if defined(SIXTY_FOUR_BIT) || defined(SIXTY_FOUR_BIT_LONG)
208static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
209                            const BN_ULONG b)
210{
211    register BN_ULONG h, l, s;
212    BN_ULONG tab[16], top3b = a >> 61;
213    register BN_ULONG a1, a2, a4, a8;
214
215    a1 = a & (0x1FFFFFFFFFFFFFFFULL);
216    a2 = a1 << 1;
217    a4 = a2 << 1;
218    a8 = a4 << 1;
219
220    tab[0] = 0;
221    tab[1] = a1;
222    tab[2] = a2;
223    tab[3] = a1 ^ a2;
224    tab[4] = a4;
225    tab[5] = a1 ^ a4;
226    tab[6] = a2 ^ a4;
227    tab[7] = a1 ^ a2 ^ a4;
228    tab[8] = a8;
229    tab[9] = a1 ^ a8;
230    tab[10] = a2 ^ a8;
231    tab[11] = a1 ^ a2 ^ a8;
232    tab[12] = a4 ^ a8;
233    tab[13] = a1 ^ a4 ^ a8;
234    tab[14] = a2 ^ a4 ^ a8;
235    tab[15] = a1 ^ a2 ^ a4 ^ a8;
236
237    s = tab[b & 0xF];
238    l = s;
239    s = tab[b >> 4 & 0xF];
240    l ^= s << 4;
241    h = s >> 60;
242    s = tab[b >> 8 & 0xF];
243    l ^= s << 8;
244    h ^= s >> 56;
245    s = tab[b >> 12 & 0xF];
246    l ^= s << 12;
247    h ^= s >> 52;
248    s = tab[b >> 16 & 0xF];
249    l ^= s << 16;
250    h ^= s >> 48;
251    s = tab[b >> 20 & 0xF];
252    l ^= s << 20;
253    h ^= s >> 44;
254    s = tab[b >> 24 & 0xF];
255    l ^= s << 24;
256    h ^= s >> 40;
257    s = tab[b >> 28 & 0xF];
258    l ^= s << 28;
259    h ^= s >> 36;
260    s = tab[b >> 32 & 0xF];
261    l ^= s << 32;
262    h ^= s >> 32;
263    s = tab[b >> 36 & 0xF];
264    l ^= s << 36;
265    h ^= s >> 28;
266    s = tab[b >> 40 & 0xF];
267    l ^= s << 40;
268    h ^= s >> 24;
269    s = tab[b >> 44 & 0xF];
270    l ^= s << 44;
271    h ^= s >> 20;
272    s = tab[b >> 48 & 0xF];
273    l ^= s << 48;
274    h ^= s >> 16;
275    s = tab[b >> 52 & 0xF];
276    l ^= s << 52;
277    h ^= s >> 12;
278    s = tab[b >> 56 & 0xF];
279    l ^= s << 56;
280    h ^= s >> 8;
281    s = tab[b >> 60];
282    l ^= s << 60;
283    h ^= s >> 4;
284
285    /* compensate for the top three bits of a */
286
287    if (top3b & 01) {
288        l ^= b << 61;
289        h ^= b >> 3;
290    }
291    if (top3b & 02) {
292        l ^= b << 62;
293        h ^= b >> 2;
294    }
295    if (top3b & 04) {
296        l ^= b << 63;
297        h ^= b >> 1;
298    }
299
300    *r1 = h;
301    *r0 = l;
302}
303#  endif
304
305/*
306 * Product of two polynomials a, b each with degree < 2 * BN_BITS2 - 1,
307 * result is a polynomial r with degree < 4 * BN_BITS2 - 1 The caller MUST
308 * ensure that the variables have the right amount of space allocated.
309 */
310static void bn_GF2m_mul_2x2(BN_ULONG *r, const BN_ULONG a1, const BN_ULONG a0,
311                            const BN_ULONG b1, const BN_ULONG b0)
312{
313    BN_ULONG m1, m0;
314    /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */
315    bn_GF2m_mul_1x1(r + 3, r + 2, a1, b1);
316    bn_GF2m_mul_1x1(r + 1, r, a0, b0);
317    bn_GF2m_mul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1);
318    /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */
319    r[2] ^= m1 ^ r[1] ^ r[3];   /* h0 ^= m1 ^ l1 ^ h1; */
320    r[1] = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0; /* l1 ^= l0 ^ h0 ^ m0; */
321}
322# else
323void bn_GF2m_mul_2x2(BN_ULONG *r, BN_ULONG a1, BN_ULONG a0, BN_ULONG b1,
324                     BN_ULONG b0);
325# endif
326
327/*
328 * Add polynomials a and b and store result in r; r could be a or b, a and b
329 * could be equal; r is the bitwise XOR of a and b.
330 */
331int BN_GF2m_add(BIGNUM *r, const BIGNUM *a, const BIGNUM *b)
332{
333    int i;
334    const BIGNUM *at, *bt;
335
336    bn_check_top(a);
337    bn_check_top(b);
338
339    if (a->top < b->top) {
340        at = b;
341        bt = a;
342    } else {
343        at = a;
344        bt = b;
345    }
346
347    if (bn_wexpand(r, at->top) == NULL)
348        return 0;
349
350    for (i = 0; i < bt->top; i++) {
351        r->d[i] = at->d[i] ^ bt->d[i];
352    }
353    for (; i < at->top; i++) {
354        r->d[i] = at->d[i];
355    }
356
357    r->top = at->top;
358    bn_correct_top(r);
359
360    return 1;
361}
362
363/*-
364 * Some functions allow for representation of the irreducible polynomials
365 * as an int[], say p.  The irreducible f(t) is then of the form:
366 *     t^p[0] + t^p[1] + ... + t^p[k]
367 * where m = p[0] > p[1] > ... > p[k] = 0.
368 */
369
370/* Performs modular reduction of a and store result in r.  r could be a. */
371int BN_GF2m_mod_arr(BIGNUM *r, const BIGNUM *a, const int p[])
372{
373    int j, k;
374    int n, dN, d0, d1;
375    BN_ULONG zz, *z;
376
377    bn_check_top(a);
378
379    if (!p[0]) {
380        /* reduction mod 1 => return 0 */
381        BN_zero(r);
382        return 1;
383    }
384
385    /*
386     * Since the algorithm does reduction in the r value, if a != r, copy the
387     * contents of a into r so we can do reduction in r.
388     */
389    if (a != r) {
390        if (!bn_wexpand(r, a->top))
391            return 0;
392        for (j = 0; j < a->top; j++) {
393            r->d[j] = a->d[j];
394        }
395        r->top = a->top;
396    }
397    z = r->d;
398
399    /* start reduction */
400    dN = p[0] / BN_BITS2;
401    for (j = r->top - 1; j > dN;) {
402        zz = z[j];
403        if (z[j] == 0) {
404            j--;
405            continue;
406        }
407        z[j] = 0;
408
409        for (k = 1; p[k] != 0; k++) {
410            /* reducing component t^p[k] */
411            n = p[0] - p[k];
412            d0 = n % BN_BITS2;
413            d1 = BN_BITS2 - d0;
414            n /= BN_BITS2;
415            z[j - n] ^= (zz >> d0);
416            if (d0)
417                z[j - n - 1] ^= (zz << d1);
418        }
419
420        /* reducing component t^0 */
421        n = dN;
422        d0 = p[0] % BN_BITS2;
423        d1 = BN_BITS2 - d0;
424        z[j - n] ^= (zz >> d0);
425        if (d0)
426            z[j - n - 1] ^= (zz << d1);
427    }
428
429    /* final round of reduction */
430    while (j == dN) {
431
432        d0 = p[0] % BN_BITS2;
433        zz = z[dN] >> d0;
434        if (zz == 0)
435            break;
436        d1 = BN_BITS2 - d0;
437
438        /* clear up the top d1 bits */
439        if (d0)
440            z[dN] = (z[dN] << d1) >> d1;
441        else
442            z[dN] = 0;
443        z[0] ^= zz;             /* reduction t^0 component */
444
445        for (k = 1; p[k] != 0; k++) {
446            BN_ULONG tmp_ulong;
447
448            /* reducing component t^p[k] */
449            n = p[k] / BN_BITS2;
450            d0 = p[k] % BN_BITS2;
451            d1 = BN_BITS2 - d0;
452            z[n] ^= (zz << d0);
453            if (d0 && (tmp_ulong = zz >> d1))
454                z[n + 1] ^= tmp_ulong;
455        }
456
457    }
458
459    bn_correct_top(r);
460    return 1;
461}
462
463/*
464 * Performs modular reduction of a by p and store result in r.  r could be a.
465 * This function calls down to the BN_GF2m_mod_arr implementation; this wrapper
466 * function is only provided for convenience; for best performance, use the
467 * BN_GF2m_mod_arr function.
468 */
469int BN_GF2m_mod(BIGNUM *r, const BIGNUM *a, const BIGNUM *p)
470{
471    int ret = 0;
472    int arr[6];
473    bn_check_top(a);
474    bn_check_top(p);
475    ret = BN_GF2m_poly2arr(p, arr, sizeof(arr) / sizeof(arr[0]));
476    if (!ret || ret > (int)(sizeof(arr) / sizeof(arr[0]))) {
477        BNerr(BN_F_BN_GF2M_MOD, BN_R_INVALID_LENGTH);
478        return 0;
479    }
480    ret = BN_GF2m_mod_arr(r, a, arr);
481    bn_check_top(r);
482    return ret;
483}
484
485/*
486 * Compute the product of two polynomials a and b, reduce modulo p, and store
487 * the result in r.  r could be a or b; a could be b.
488 */
489int BN_GF2m_mod_mul_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
490                        const int p[], BN_CTX *ctx)
491{
492    int zlen, i, j, k, ret = 0;
493    BIGNUM *s;
494    BN_ULONG x1, x0, y1, y0, zz[4];
495
496    bn_check_top(a);
497    bn_check_top(b);
498
499    if (a == b) {
500        return BN_GF2m_mod_sqr_arr(r, a, p, ctx);
501    }
502
503    BN_CTX_start(ctx);
504    if ((s = BN_CTX_get(ctx)) == NULL)
505        goto err;
506
507    zlen = a->top + b->top + 4;
508    if (!bn_wexpand(s, zlen))
509        goto err;
510    s->top = zlen;
511
512    for (i = 0; i < zlen; i++)
513        s->d[i] = 0;
514
515    for (j = 0; j < b->top; j += 2) {
516        y0 = b->d[j];
517        y1 = ((j + 1) == b->top) ? 0 : b->d[j + 1];
518        for (i = 0; i < a->top; i += 2) {
519            x0 = a->d[i];
520            x1 = ((i + 1) == a->top) ? 0 : a->d[i + 1];
521            bn_GF2m_mul_2x2(zz, x1, x0, y1, y0);
522            for (k = 0; k < 4; k++)
523                s->d[i + j + k] ^= zz[k];
524        }
525    }
526
527    bn_correct_top(s);
528    if (BN_GF2m_mod_arr(r, s, p))
529        ret = 1;
530    bn_check_top(r);
531
532 err:
533    BN_CTX_end(ctx);
534    return ret;
535}
536
537/*
538 * Compute the product of two polynomials a and b, reduce modulo p, and store
539 * the result in r.  r could be a or b; a could equal b. This function calls
540 * down to the BN_GF2m_mod_mul_arr implementation; this wrapper function is
541 * only provided for convenience; for best performance, use the
542 * BN_GF2m_mod_mul_arr function.
543 */
544int BN_GF2m_mod_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
545                    const BIGNUM *p, BN_CTX *ctx)
546{
547    int ret = 0;
548    const int max = BN_num_bits(p) + 1;
549    int *arr = NULL;
550    bn_check_top(a);
551    bn_check_top(b);
552    bn_check_top(p);
553    if ((arr = (int *)OPENSSL_malloc(sizeof(int) * max)) == NULL)
554        goto err;
555    ret = BN_GF2m_poly2arr(p, arr, max);
556    if (!ret || ret > max) {
557        BNerr(BN_F_BN_GF2M_MOD_MUL, BN_R_INVALID_LENGTH);
558        goto err;
559    }
560    ret = BN_GF2m_mod_mul_arr(r, a, b, arr, ctx);
561    bn_check_top(r);
562 err:
563    if (arr)
564        OPENSSL_free(arr);
565    return ret;
566}
567
568/* Square a, reduce the result mod p, and store it in a.  r could be a. */
569int BN_GF2m_mod_sqr_arr(BIGNUM *r, const BIGNUM *a, const int p[],
570                        BN_CTX *ctx)
571{
572    int i, ret = 0;
573    BIGNUM *s;
574
575    bn_check_top(a);
576    BN_CTX_start(ctx);
577    if ((s = BN_CTX_get(ctx)) == NULL)
578        goto err;
579    if (!bn_wexpand(s, 2 * a->top))
580        goto err;
581
582    for (i = a->top - 1; i >= 0; i--) {
583        s->d[2 * i + 1] = SQR1(a->d[i]);
584        s->d[2 * i] = SQR0(a->d[i]);
585    }
586
587    s->top = 2 * a->top;
588    bn_correct_top(s);
589    if (!BN_GF2m_mod_arr(r, s, p))
590        goto err;
591    bn_check_top(r);
592    ret = 1;
593 err:
594    BN_CTX_end(ctx);
595    return ret;
596}
597
598/*
599 * Square a, reduce the result mod p, and store it in a.  r could be a. This
600 * function calls down to the BN_GF2m_mod_sqr_arr implementation; this
601 * wrapper function is only provided for convenience; for best performance,
602 * use the BN_GF2m_mod_sqr_arr function.
603 */
604int BN_GF2m_mod_sqr(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
605{
606    int ret = 0;
607    const int max = BN_num_bits(p) + 1;
608    int *arr = NULL;
609
610    bn_check_top(a);
611    bn_check_top(p);
612    if ((arr = (int *)OPENSSL_malloc(sizeof(int) * max)) == NULL)
613        goto err;
614    ret = BN_GF2m_poly2arr(p, arr, max);
615    if (!ret || ret > max) {
616        BNerr(BN_F_BN_GF2M_MOD_SQR, BN_R_INVALID_LENGTH);
617        goto err;
618    }
619    ret = BN_GF2m_mod_sqr_arr(r, a, arr, ctx);
620    bn_check_top(r);
621 err:
622    if (arr)
623        OPENSSL_free(arr);
624    return ret;
625}
626
627/*
628 * Invert a, reduce modulo p, and store the result in r. r could be a. Uses
629 * Modified Almost Inverse Algorithm (Algorithm 10) from Hankerson, D.,
630 * Hernandez, J.L., and Menezes, A.  "Software Implementation of Elliptic
631 * Curve Cryptography Over Binary Fields".
632 */
633int BN_GF2m_mod_inv(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
634{
635    BIGNUM *b, *c = NULL, *u = NULL, *v = NULL, *tmp;
636    int ret = 0;
637
638    bn_check_top(a);
639    bn_check_top(p);
640
641    BN_CTX_start(ctx);
642
643    if ((b = BN_CTX_get(ctx)) == NULL)
644        goto err;
645    if ((c = BN_CTX_get(ctx)) == NULL)
646        goto err;
647    if ((u = BN_CTX_get(ctx)) == NULL)
648        goto err;
649    if ((v = BN_CTX_get(ctx)) == NULL)
650        goto err;
651
652    if (!BN_GF2m_mod(u, a, p))
653        goto err;
654    if (BN_is_zero(u))
655        goto err;
656
657    if (!BN_copy(v, p))
658        goto err;
659# if 0
660    if (!BN_one(b))
661        goto err;
662
663    while (1) {
664        while (!BN_is_odd(u)) {
665            if (BN_is_zero(u))
666                goto err;
667            if (!BN_rshift1(u, u))
668                goto err;
669            if (BN_is_odd(b)) {
670                if (!BN_GF2m_add(b, b, p))
671                    goto err;
672            }
673            if (!BN_rshift1(b, b))
674                goto err;
675        }
676
677        if (BN_abs_is_word(u, 1))
678            break;
679
680        if (BN_num_bits(u) < BN_num_bits(v)) {
681            tmp = u;
682            u = v;
683            v = tmp;
684            tmp = b;
685            b = c;
686            c = tmp;
687        }
688
689        if (!BN_GF2m_add(u, u, v))
690            goto err;
691        if (!BN_GF2m_add(b, b, c))
692            goto err;
693    }
694# else
695    {
696        int i;
697        int ubits = BN_num_bits(u);
698        int vbits = BN_num_bits(v); /* v is copy of p */
699        int top = p->top;
700        BN_ULONG *udp, *bdp, *vdp, *cdp;
701
702        if (!bn_wexpand(u, top))
703            goto err;
704        udp = u->d;
705        for (i = u->top; i < top; i++)
706            udp[i] = 0;
707        u->top = top;
708        if (!bn_wexpand(b, top))
709          goto err;
710        bdp = b->d;
711        bdp[0] = 1;
712        for (i = 1; i < top; i++)
713            bdp[i] = 0;
714        b->top = top;
715        if (!bn_wexpand(c, top))
716          goto err;
717        cdp = c->d;
718        for (i = 0; i < top; i++)
719            cdp[i] = 0;
720        c->top = top;
721        vdp = v->d;             /* It pays off to "cache" *->d pointers,
722                                 * because it allows optimizer to be more
723                                 * aggressive. But we don't have to "cache"
724                                 * p->d, because *p is declared 'const'... */
725        while (1) {
726            while (ubits && !(udp[0] & 1)) {
727                BN_ULONG u0, u1, b0, b1, mask;
728
729                u0 = udp[0];
730                b0 = bdp[0];
731                mask = (BN_ULONG)0 - (b0 & 1);
732                b0 ^= p->d[0] & mask;
733                for (i = 0; i < top - 1; i++) {
734                    u1 = udp[i + 1];
735                    udp[i] = ((u0 >> 1) | (u1 << (BN_BITS2 - 1))) & BN_MASK2;
736                    u0 = u1;
737                    b1 = bdp[i + 1] ^ (p->d[i + 1] & mask);
738                    bdp[i] = ((b0 >> 1) | (b1 << (BN_BITS2 - 1))) & BN_MASK2;
739                    b0 = b1;
740                }
741                udp[i] = u0 >> 1;
742                bdp[i] = b0 >> 1;
743                ubits--;
744            }
745
746            if (ubits <= BN_BITS2) {
747                if (udp[0] == 0) /* poly was reducible */
748                    goto err;
749                if (udp[0] == 1)
750                    break;
751            }
752
753            if (ubits < vbits) {
754                i = ubits;
755                ubits = vbits;
756                vbits = i;
757                tmp = u;
758                u = v;
759                v = tmp;
760                tmp = b;
761                b = c;
762                c = tmp;
763                udp = vdp;
764                vdp = v->d;
765                bdp = cdp;
766                cdp = c->d;
767            }
768            for (i = 0; i < top; i++) {
769                udp[i] ^= vdp[i];
770                bdp[i] ^= cdp[i];
771            }
772            if (ubits == vbits) {
773                BN_ULONG ul;
774                int utop = (ubits - 1) / BN_BITS2;
775
776                while ((ul = udp[utop]) == 0 && utop)
777                    utop--;
778                ubits = utop * BN_BITS2 + BN_num_bits_word(ul);
779            }
780        }
781        bn_correct_top(b);
782    }
783# endif
784
785    if (!BN_copy(r, b))
786        goto err;
787    bn_check_top(r);
788    ret = 1;
789
790 err:
791# ifdef BN_DEBUG                /* BN_CTX_end would complain about the
792                                 * expanded form */
793    bn_correct_top(c);
794    bn_correct_top(u);
795    bn_correct_top(v);
796# endif
797    BN_CTX_end(ctx);
798    return ret;
799}
800
801/*
802 * Invert xx, reduce modulo p, and store the result in r. r could be xx.
803 * This function calls down to the BN_GF2m_mod_inv implementation; this
804 * wrapper function is only provided for convenience; for best performance,
805 * use the BN_GF2m_mod_inv function.
806 */
807int BN_GF2m_mod_inv_arr(BIGNUM *r, const BIGNUM *xx, const int p[],
808                        BN_CTX *ctx)
809{
810    BIGNUM *field;
811    int ret = 0;
812
813    bn_check_top(xx);
814    BN_CTX_start(ctx);
815    if ((field = BN_CTX_get(ctx)) == NULL)
816        goto err;
817    if (!BN_GF2m_arr2poly(p, field))
818        goto err;
819
820    ret = BN_GF2m_mod_inv(r, xx, field, ctx);
821    bn_check_top(r);
822
823 err:
824    BN_CTX_end(ctx);
825    return ret;
826}
827
828# ifndef OPENSSL_SUN_GF2M_DIV
829/*
830 * Divide y by x, reduce modulo p, and store the result in r. r could be x
831 * or y, x could equal y.
832 */
833int BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x,
834                    const BIGNUM *p, BN_CTX *ctx)
835{
836    BIGNUM *xinv = NULL;
837    int ret = 0;
838
839    bn_check_top(y);
840    bn_check_top(x);
841    bn_check_top(p);
842
843    BN_CTX_start(ctx);
844    xinv = BN_CTX_get(ctx);
845    if (xinv == NULL)
846        goto err;
847
848    if (!BN_GF2m_mod_inv(xinv, x, p, ctx))
849        goto err;
850    if (!BN_GF2m_mod_mul(r, y, xinv, p, ctx))
851        goto err;
852    bn_check_top(r);
853    ret = 1;
854
855 err:
856    BN_CTX_end(ctx);
857    return ret;
858}
859# else
860/*
861 * Divide y by x, reduce modulo p, and store the result in r. r could be x
862 * or y, x could equal y. Uses algorithm Modular_Division_GF(2^m) from
863 * Chang-Shantz, S.  "From Euclid's GCD to Montgomery Multiplication to the
864 * Great Divide".
865 */
866int BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x,
867                    const BIGNUM *p, BN_CTX *ctx)
868{
869    BIGNUM *a, *b, *u, *v;
870    int ret = 0;
871
872    bn_check_top(y);
873    bn_check_top(x);
874    bn_check_top(p);
875
876    BN_CTX_start(ctx);
877
878    a = BN_CTX_get(ctx);
879    b = BN_CTX_get(ctx);
880    u = BN_CTX_get(ctx);
881    v = BN_CTX_get(ctx);
882    if (v == NULL)
883        goto err;
884
885    /* reduce x and y mod p */
886    if (!BN_GF2m_mod(u, y, p))
887        goto err;
888    if (!BN_GF2m_mod(a, x, p))
889        goto err;
890    if (!BN_copy(b, p))
891        goto err;
892
893    while (!BN_is_odd(a)) {
894        if (!BN_rshift1(a, a))
895            goto err;
896        if (BN_is_odd(u))
897            if (!BN_GF2m_add(u, u, p))
898                goto err;
899        if (!BN_rshift1(u, u))
900            goto err;
901    }
902
903    do {
904        if (BN_GF2m_cmp(b, a) > 0) {
905            if (!BN_GF2m_add(b, b, a))
906                goto err;
907            if (!BN_GF2m_add(v, v, u))
908                goto err;
909            do {
910                if (!BN_rshift1(b, b))
911                    goto err;
912                if (BN_is_odd(v))
913                    if (!BN_GF2m_add(v, v, p))
914                        goto err;
915                if (!BN_rshift1(v, v))
916                    goto err;
917            } while (!BN_is_odd(b));
918        } else if (BN_abs_is_word(a, 1))
919            break;
920        else {
921            if (!BN_GF2m_add(a, a, b))
922                goto err;
923            if (!BN_GF2m_add(u, u, v))
924                goto err;
925            do {
926                if (!BN_rshift1(a, a))
927                    goto err;
928                if (BN_is_odd(u))
929                    if (!BN_GF2m_add(u, u, p))
930                        goto err;
931                if (!BN_rshift1(u, u))
932                    goto err;
933            } while (!BN_is_odd(a));
934        }
935    } while (1);
936
937    if (!BN_copy(r, u))
938        goto err;
939    bn_check_top(r);
940    ret = 1;
941
942 err:
943    BN_CTX_end(ctx);
944    return ret;
945}
946# endif
947
948/*
949 * Divide yy by xx, reduce modulo p, and store the result in r. r could be xx
950 * * or yy, xx could equal yy. This function calls down to the
951 * BN_GF2m_mod_div implementation; this wrapper function is only provided for
952 * convenience; for best performance, use the BN_GF2m_mod_div function.
953 */
954int BN_GF2m_mod_div_arr(BIGNUM *r, const BIGNUM *yy, const BIGNUM *xx,
955                        const int p[], BN_CTX *ctx)
956{
957    BIGNUM *field;
958    int ret = 0;
959
960    bn_check_top(yy);
961    bn_check_top(xx);
962
963    BN_CTX_start(ctx);
964    if ((field = BN_CTX_get(ctx)) == NULL)
965        goto err;
966    if (!BN_GF2m_arr2poly(p, field))
967        goto err;
968
969    ret = BN_GF2m_mod_div(r, yy, xx, field, ctx);
970    bn_check_top(r);
971
972 err:
973    BN_CTX_end(ctx);
974    return ret;
975}
976
977/*
978 * Compute the bth power of a, reduce modulo p, and store the result in r.  r
979 * could be a. Uses simple square-and-multiply algorithm A.5.1 from IEEE
980 * P1363.
981 */
982int BN_GF2m_mod_exp_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
983                        const int p[], BN_CTX *ctx)
984{
985    int ret = 0, i, n;
986    BIGNUM *u;
987
988    bn_check_top(a);
989    bn_check_top(b);
990
991    if (BN_is_zero(b))
992        return (BN_one(r));
993
994    if (BN_abs_is_word(b, 1))
995        return (BN_copy(r, a) != NULL);
996
997    BN_CTX_start(ctx);
998    if ((u = BN_CTX_get(ctx)) == NULL)
999        goto err;
1000
1001    if (!BN_GF2m_mod_arr(u, a, p))
1002        goto err;
1003
1004    n = BN_num_bits(b) - 1;
1005    for (i = n - 1; i >= 0; i--) {
1006        if (!BN_GF2m_mod_sqr_arr(u, u, p, ctx))
1007            goto err;
1008        if (BN_is_bit_set(b, i)) {
1009            if (!BN_GF2m_mod_mul_arr(u, u, a, p, ctx))
1010                goto err;
1011        }
1012    }
1013    if (!BN_copy(r, u))
1014        goto err;
1015    bn_check_top(r);
1016    ret = 1;
1017 err:
1018    BN_CTX_end(ctx);
1019    return ret;
1020}
1021
1022/*
1023 * Compute the bth power of a, reduce modulo p, and store the result in r.  r
1024 * could be a. This function calls down to the BN_GF2m_mod_exp_arr
1025 * implementation; this wrapper function is only provided for convenience;
1026 * for best performance, use the BN_GF2m_mod_exp_arr function.
1027 */
1028int BN_GF2m_mod_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
1029                    const BIGNUM *p, BN_CTX *ctx)
1030{
1031    int ret = 0;
1032    const int max = BN_num_bits(p) + 1;
1033    int *arr = NULL;
1034    bn_check_top(a);
1035    bn_check_top(b);
1036    bn_check_top(p);
1037    if ((arr = (int *)OPENSSL_malloc(sizeof(int) * max)) == NULL)
1038        goto err;
1039    ret = BN_GF2m_poly2arr(p, arr, max);
1040    if (!ret || ret > max) {
1041        BNerr(BN_F_BN_GF2M_MOD_EXP, BN_R_INVALID_LENGTH);
1042        goto err;
1043    }
1044    ret = BN_GF2m_mod_exp_arr(r, a, b, arr, ctx);
1045    bn_check_top(r);
1046 err:
1047    if (arr)
1048        OPENSSL_free(arr);
1049    return ret;
1050}
1051
1052/*
1053 * Compute the square root of a, reduce modulo p, and store the result in r.
1054 * r could be a. Uses exponentiation as in algorithm A.4.1 from IEEE P1363.
1055 */
1056int BN_GF2m_mod_sqrt_arr(BIGNUM *r, const BIGNUM *a, const int p[],
1057                         BN_CTX *ctx)
1058{
1059    int ret = 0;
1060    BIGNUM *u;
1061
1062    bn_check_top(a);
1063
1064    if (!p[0]) {
1065        /* reduction mod 1 => return 0 */
1066        BN_zero(r);
1067        return 1;
1068    }
1069
1070    BN_CTX_start(ctx);
1071    if ((u = BN_CTX_get(ctx)) == NULL)
1072        goto err;
1073
1074    if (!BN_set_bit(u, p[0] - 1))
1075        goto err;
1076    ret = BN_GF2m_mod_exp_arr(r, a, u, p, ctx);
1077    bn_check_top(r);
1078
1079 err:
1080    BN_CTX_end(ctx);
1081    return ret;
1082}
1083
1084/*
1085 * Compute the square root of a, reduce modulo p, and store the result in r.
1086 * r could be a. This function calls down to the BN_GF2m_mod_sqrt_arr
1087 * implementation; this wrapper function is only provided for convenience;
1088 * for best performance, use the BN_GF2m_mod_sqrt_arr function.
1089 */
1090int BN_GF2m_mod_sqrt(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
1091{
1092    int ret = 0;
1093    const int max = BN_num_bits(p) + 1;
1094    int *arr = NULL;
1095    bn_check_top(a);
1096    bn_check_top(p);
1097    if ((arr = (int *)OPENSSL_malloc(sizeof(int) * max)) == NULL)
1098        goto err;
1099    ret = BN_GF2m_poly2arr(p, arr, max);
1100    if (!ret || ret > max) {
1101        BNerr(BN_F_BN_GF2M_MOD_SQRT, BN_R_INVALID_LENGTH);
1102        goto err;
1103    }
1104    ret = BN_GF2m_mod_sqrt_arr(r, a, arr, ctx);
1105    bn_check_top(r);
1106 err:
1107    if (arr)
1108        OPENSSL_free(arr);
1109    return ret;
1110}
1111
1112/*
1113 * Find r such that r^2 + r = a mod p.  r could be a. If no r exists returns
1114 * 0. Uses algorithms A.4.7 and A.4.6 from IEEE P1363.
1115 */
1116int BN_GF2m_mod_solve_quad_arr(BIGNUM *r, const BIGNUM *a_, const int p[],
1117                               BN_CTX *ctx)
1118{
1119    int ret = 0, count = 0, j;
1120    BIGNUM *a, *z, *rho, *w, *w2, *tmp;
1121
1122    bn_check_top(a_);
1123
1124    if (!p[0]) {
1125        /* reduction mod 1 => return 0 */
1126        BN_zero(r);
1127        return 1;
1128    }
1129
1130    BN_CTX_start(ctx);
1131    a = BN_CTX_get(ctx);
1132    z = BN_CTX_get(ctx);
1133    w = BN_CTX_get(ctx);
1134    if (w == NULL)
1135        goto err;
1136
1137    if (!BN_GF2m_mod_arr(a, a_, p))
1138        goto err;
1139
1140    if (BN_is_zero(a)) {
1141        BN_zero(r);
1142        ret = 1;
1143        goto err;
1144    }
1145
1146    if (p[0] & 0x1) {           /* m is odd */
1147        /* compute half-trace of a */
1148        if (!BN_copy(z, a))
1149            goto err;
1150        for (j = 1; j <= (p[0] - 1) / 2; j++) {
1151            if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1152                goto err;
1153            if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1154                goto err;
1155            if (!BN_GF2m_add(z, z, a))
1156                goto err;
1157        }
1158
1159    } else {                    /* m is even */
1160
1161        rho = BN_CTX_get(ctx);
1162        w2 = BN_CTX_get(ctx);
1163        tmp = BN_CTX_get(ctx);
1164        if (tmp == NULL)
1165            goto err;
1166        do {
1167            if (!BN_rand(rho, p[0], 0, 0))
1168                goto err;
1169            if (!BN_GF2m_mod_arr(rho, rho, p))
1170                goto err;
1171            BN_zero(z);
1172            if (!BN_copy(w, rho))
1173                goto err;
1174            for (j = 1; j <= p[0] - 1; j++) {
1175                if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1176                    goto err;
1177                if (!BN_GF2m_mod_sqr_arr(w2, w, p, ctx))
1178                    goto err;
1179                if (!BN_GF2m_mod_mul_arr(tmp, w2, a, p, ctx))
1180                    goto err;
1181                if (!BN_GF2m_add(z, z, tmp))
1182                    goto err;
1183                if (!BN_GF2m_add(w, w2, rho))
1184                    goto err;
1185            }
1186            count++;
1187        } while (BN_is_zero(w) && (count < MAX_ITERATIONS));
1188        if (BN_is_zero(w)) {
1189            BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD_ARR, BN_R_TOO_MANY_ITERATIONS);
1190            goto err;
1191        }
1192    }
1193
1194    if (!BN_GF2m_mod_sqr_arr(w, z, p, ctx))
1195        goto err;
1196    if (!BN_GF2m_add(w, z, w))
1197        goto err;
1198    if (BN_GF2m_cmp(w, a)) {
1199        BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD_ARR, BN_R_NO_SOLUTION);
1200        goto err;
1201    }
1202
1203    if (!BN_copy(r, z))
1204        goto err;
1205    bn_check_top(r);
1206
1207    ret = 1;
1208
1209 err:
1210    BN_CTX_end(ctx);
1211    return ret;
1212}
1213
1214/*
1215 * Find r such that r^2 + r = a mod p.  r could be a. If no r exists returns
1216 * 0. This function calls down to the BN_GF2m_mod_solve_quad_arr
1217 * implementation; this wrapper function is only provided for convenience;
1218 * for best performance, use the BN_GF2m_mod_solve_quad_arr function.
1219 */
1220int BN_GF2m_mod_solve_quad(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
1221                           BN_CTX *ctx)
1222{
1223    int ret = 0;
1224    const int max = BN_num_bits(p) + 1;
1225    int *arr = NULL;
1226    bn_check_top(a);
1227    bn_check_top(p);
1228    if ((arr = (int *)OPENSSL_malloc(sizeof(int) * max)) == NULL)
1229        goto err;
1230    ret = BN_GF2m_poly2arr(p, arr, max);
1231    if (!ret || ret > max) {
1232        BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD, BN_R_INVALID_LENGTH);
1233        goto err;
1234    }
1235    ret = BN_GF2m_mod_solve_quad_arr(r, a, arr, ctx);
1236    bn_check_top(r);
1237 err:
1238    if (arr)
1239        OPENSSL_free(arr);
1240    return ret;
1241}
1242
1243/*
1244 * Convert the bit-string representation of a polynomial ( \sum_{i=0}^n a_i *
1245 * x^i) into an array of integers corresponding to the bits with non-zero
1246 * coefficient.  Array is terminated with -1. Up to max elements of the array
1247 * will be filled.  Return value is total number of array elements that would
1248 * be filled if array was large enough.
1249 */
1250int BN_GF2m_poly2arr(const BIGNUM *a, int p[], int max)
1251{
1252    int i, j, k = 0;
1253    BN_ULONG mask;
1254
1255    if (BN_is_zero(a))
1256        return 0;
1257
1258    for (i = a->top - 1; i >= 0; i--) {
1259        if (!a->d[i])
1260            /* skip word if a->d[i] == 0 */
1261            continue;
1262        mask = BN_TBIT;
1263        for (j = BN_BITS2 - 1; j >= 0; j--) {
1264            if (a->d[i] & mask) {
1265                if (k < max)
1266                    p[k] = BN_BITS2 * i + j;
1267                k++;
1268            }
1269            mask >>= 1;
1270        }
1271    }
1272
1273    if (k < max) {
1274        p[k] = -1;
1275        k++;
1276    }
1277
1278    return k;
1279}
1280
1281/*
1282 * Convert the coefficient array representation of a polynomial to a
1283 * bit-string.  The array must be terminated by -1.
1284 */
1285int BN_GF2m_arr2poly(const int p[], BIGNUM *a)
1286{
1287    int i;
1288
1289    bn_check_top(a);
1290    BN_zero(a);
1291    for (i = 0; p[i] != -1; i++) {
1292        if (BN_set_bit(a, p[i]) == 0)
1293            return 0;
1294    }
1295    bn_check_top(a);
1296
1297    return 1;
1298}
1299
1300#endif
1301