crypto/openssl/crypto/bn/bn_gf2m.c

   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 /*
  99  * Maximum number of iterations before BN_GF2m_mod_solve_quad_arr should
 100  * fail.
 101  */
 102 #define MAX_ITERATIONS 50
 103
 104 static const BN_ULONG SQR_tb[16] = { 0, 1, 4, 5, 16, 17, 20, 21,
 105     64, 65, 68, 69, 80, 81, 84, 85
 106 };
 107
 108 /* Platform-specific macros to accelerate squaring. */
 109 #if defined(SIXTY_FOUR_BIT) || defined(SIXTY_FOUR_BIT_LONG)
 110 # define SQR1(w) \
 111     SQR_tb[(w) >> 60 & 0xF] << 56 | SQR_tb[(w) >> 56 & 0xF] << 48 | \
 112     SQR_tb[(w) >> 52 & 0xF] << 40 | SQR_tb[(w) >> 48 & 0xF] << 32 | \
 113     SQR_tb[(w) >> 44 & 0xF] << 24 | SQR_tb[(w) >> 40 & 0xF] << 16 | \
 114     SQR_tb[(w) >> 36 & 0xF] <<  8 | SQR_tb[(w) >> 32 & 0xF]
 115 # define SQR0(w) \
 116     SQR_tb[(w) >> 28 & 0xF] << 56 | SQR_tb[(w) >> 24 & 0xF] << 48 | \
 117     SQR_tb[(w) >> 20 & 0xF] << 40 | SQR_tb[(w) >> 16 & 0xF] << 32 | \
 118     SQR_tb[(w) >> 12 & 0xF] << 24 | SQR_tb[(w) >>  8 & 0xF] << 16 | \
 119     SQR_tb[(w) >>  4 & 0xF] <<  8 | SQR_tb[(w)       & 0xF]
 120 #endif
 121 #ifdef THIRTY_TWO_BIT
 122 # define SQR1(w) \
 123     SQR_tb[(w) >> 28 & 0xF] << 24 | SQR_tb[(w) >> 24 & 0xF] << 16 | \
 124     SQR_tb[(w) >> 20 & 0xF] <<  8 | SQR_tb[(w) >> 16 & 0xF]
 125 # define SQR0(w) \
 126     SQR_tb[(w) >> 12 & 0xF] << 24 | SQR_tb[(w) >>  8 & 0xF] << 16 | \
 127     SQR_tb[(w) >>  4 & 0xF] <<  8 | SQR_tb[(w)       & 0xF]
 128 #endif
 129 #ifdef SIXTEEN_BIT
 130 # define SQR1(w) \
 131     SQR_tb[(w) >> 12 & 0xF] <<  8 | SQR_tb[(w) >>  8 & 0xF]
 132 # define SQR0(w) \
 133     SQR_tb[(w) >>  4 & 0xF] <<  8 | SQR_tb[(w)       & 0xF]
 134 #endif
 135 #ifdef EIGHT_BIT
 136 # define SQR1(w) \
 137     SQR_tb[(w) >>  4 & 0xF]
 138 # define SQR0(w) \
 139     SQR_tb[(w)       & 15]
 140 #endif
 141
 142 /*
 143  * Product of two polynomials a, b each with degree < BN_BITS2 - 1, result is
 144  * a polynomial r with degree < 2 * BN_BITS - 1 The caller MUST ensure that
 145  * the variables have the right amount of space allocated.
 146  */
 147 #ifdef EIGHT_BIT
 148 static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
 149                             const BN_ULONG b)
 150 {
 151     register BN_ULONG h, l, s;
 152     BN_ULONG tab[4], top1b = a >> 7;
 153     register BN_ULONG a1, a2;
 154
 155     a1 = a & (0x7F);
 156     a2 = a1 << 1;
 157
 158     tab[0] = 0;
 159     tab[1] = a1;
 160     tab[2] = a2;
 161     tab[3] = a1 ^ a2;
 162
 163     s = tab[b & 0x3];
 164     l = s;
 165     s = tab[b >> 2 & 0x3];
 166     l ^= s << 2;
 167     h = s >> 6;
 168     s = tab[b >> 4 & 0x3];
 169     l ^= s << 4;
 170     h ^= s >> 4;
 171     s = tab[b >> 6];
 172     l ^= s << 6;
 173     h ^= s >> 2;
 174
 175     /* compensate for the top bit of a */
 176
 177     if (top1b & 01) {
 178         l ^= b << 7;
 179         h ^= b >> 1;
 180     }
 181
 182     *r1 = h;
 183     *r0 = l;
 184 }
 185 #endif
 186 #ifdef SIXTEEN_BIT
 187 static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
 188                             const BN_ULONG b)
 189 {
 190     register BN_ULONG h, l, s;
 191     BN_ULONG tab[4], top1b = a >> 15;
 192     register BN_ULONG a1, a2;
 193
 194     a1 = a & (0x7FFF);
 195     a2 = a1 << 1;
 196
 197     tab[0] = 0;
 198     tab[1] = a1;
 199     tab[2] = a2;
 200     tab[3] = a1 ^ a2;
 201
 202     s = tab[b & 0x3];
 203     l = s;
 204     s = tab[b >> 2 & 0x3];
 205     l ^= s << 2;
 206     h = s >> 14;
 207     s = tab[b >> 4 & 0x3];
 208     l ^= s << 4;
 209     h ^= s >> 12;
 210     s = tab[b >> 6 & 0x3];
 211     l ^= s << 6;
 212     h ^= s >> 10;
 213     s = tab[b >> 8 & 0x3];
 214     l ^= s << 8;
 215     h ^= s >> 8;
 216     s = tab[b >> 10 & 0x3];
 217     l ^= s << 10;
 218     h ^= s >> 6;
 219     s = tab[b >> 12 & 0x3];
 220     l ^= s << 12;
 221     h ^= s >> 4;
 222     s = tab[b >> 14];
 223     l ^= s << 14;
 224     h ^= s >> 2;
 225
 226     /* compensate for the top bit of a */
 227
 228     if (top1b & 01) {
 229         l ^= b << 15;
 230         h ^= b >> 1;
 231     }
 232
 233     *r1 = h;
 234     *r0 = l;
 235 }
 236 #endif
 237 #ifdef THIRTY_TWO_BIT
 238 static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
 239                             const BN_ULONG b)
 240 {
 241     register BN_ULONG h, l, s;
 242     BN_ULONG tab[8], top2b = a >> 30;
 243     register BN_ULONG a1, a2, a4;
 244
 245     a1 = a & (0x3FFFFFFF);
 246     a2 = a1 << 1;
 247     a4 = a2 << 1;
 248
 249     tab[0] = 0;
 250     tab[1] = a1;
 251     tab[2] = a2;
 252     tab[3] = a1 ^ a2;
 253     tab[4] = a4;
 254     tab[5] = a1 ^ a4;
 255     tab[6] = a2 ^ a4;
 256     tab[7] = a1 ^ a2 ^ a4;
 257
 258     s = tab[b & 0x7];
 259     l = s;
 260     s = tab[b >> 3 & 0x7];
 261     l ^= s << 3;
 262     h = s >> 29;
 263     s = tab[b >> 6 & 0x7];
 264     l ^= s << 6;
 265     h ^= s >> 26;
 266     s = tab[b >> 9 & 0x7];
 267     l ^= s << 9;
 268     h ^= s >> 23;
 269     s = tab[b >> 12 & 0x7];
 270     l ^= s << 12;
 271     h ^= s >> 20;
 272     s = tab[b >> 15 & 0x7];
 273     l ^= s << 15;
 274     h ^= s >> 17;
 275     s = tab[b >> 18 & 0x7];
 276     l ^= s << 18;
 277     h ^= s >> 14;
 278     s = tab[b >> 21 & 0x7];
 279     l ^= s << 21;
 280     h ^= s >> 11;
 281     s = tab[b >> 24 & 0x7];
 282     l ^= s << 24;
 283     h ^= s >> 8;
 284     s = tab[b >> 27 & 0x7];
 285     l ^= s << 27;
 286     h ^= s >> 5;
 287     s = tab[b >> 30];
 288     l ^= s << 30;
 289     h ^= s >> 2;
 290
 291     /* compensate for the top two bits of a */
 292
 293     if (top2b & 01) {
 294         l ^= b << 30;
 295         h ^= b >> 2;
 296     }
 297     if (top2b & 02) {
 298         l ^= b << 31;
 299         h ^= b >> 1;
 300     }
 301
 302     *r1 = h;
 303     *r0 = l;
 304 }
 305 #endif
 306 #if defined(SIXTY_FOUR_BIT) || defined(SIXTY_FOUR_BIT_LONG)
 307 static void bn_GF2m_mul_1x1(BN_ULONG *r1, BN_ULONG *r0, const BN_ULONG a,
 308                             const BN_ULONG b)
 309 {
 310     register BN_ULONG h, l, s;
 311     BN_ULONG tab[16], top3b = a >> 61;
 312     register BN_ULONG a1, a2, a4, a8;
 313
 314     a1 = a & (0x1FFFFFFFFFFFFFFFULL);
 315     a2 = a1 << 1;
 316     a4 = a2 << 1;
 317     a8 = a4 << 1;
 318
 319     tab[0] = 0;
 320     tab[1] = a1;
 321     tab[2] = a2;
 322     tab[3] = a1 ^ a2;
 323     tab[4] = a4;
 324     tab[5] = a1 ^ a4;
 325     tab[6] = a2 ^ a4;
 326     tab[7] = a1 ^ a2 ^ a4;
 327     tab[8] = a8;
 328     tab[9] = a1 ^ a8;
 329     tab[10] = a2 ^ a8;
 330     tab[11] = a1 ^ a2 ^ a8;
 331     tab[12] = a4 ^ a8;
 332     tab[13] = a1 ^ a4 ^ a8;
 333     tab[14] = a2 ^ a4 ^ a8;
 334     tab[15] = a1 ^ a2 ^ a4 ^ a8;
 335
 336     s = tab[b & 0xF];
 337     l = s;
 338     s = tab[b >> 4 & 0xF];
 339     l ^= s << 4;
 340     h = s >> 60;
 341     s = tab[b >> 8 & 0xF];
 342     l ^= s << 8;
 343     h ^= s >> 56;
 344     s = tab[b >> 12 & 0xF];
 345     l ^= s << 12;
 346     h ^= s >> 52;
 347     s = tab[b >> 16 & 0xF];
 348     l ^= s << 16;
 349     h ^= s >> 48;
 350     s = tab[b >> 20 & 0xF];
 351     l ^= s << 20;
 352     h ^= s >> 44;
 353     s = tab[b >> 24 & 0xF];
 354     l ^= s << 24;
 355     h ^= s >> 40;
 356     s = tab[b >> 28 & 0xF];
 357     l ^= s << 28;
 358     h ^= s >> 36;
 359     s = tab[b >> 32 & 0xF];
 360     l ^= s << 32;
 361     h ^= s >> 32;
 362     s = tab[b >> 36 & 0xF];
 363     l ^= s << 36;
 364     h ^= s >> 28;
 365     s = tab[b >> 40 & 0xF];
 366     l ^= s << 40;
 367     h ^= s >> 24;
 368     s = tab[b >> 44 & 0xF];
 369     l ^= s << 44;
 370     h ^= s >> 20;
 371     s = tab[b >> 48 & 0xF];
 372     l ^= s << 48;
 373     h ^= s >> 16;
 374     s = tab[b >> 52 & 0xF];
 375     l ^= s << 52;
 376     h ^= s >> 12;
 377     s = tab[b >> 56 & 0xF];
 378     l ^= s << 56;
 379     h ^= s >> 8;
 380     s = tab[b >> 60];
 381     l ^= s << 60;
 382     h ^= s >> 4;
 383
 384     /* compensate for the top three bits of a */
 385
 386     if (top3b & 01) {
 387         l ^= b << 61;
 388         h ^= b >> 3;
 389     }
 390     if (top3b & 02) {
 391         l ^= b << 62;
 392         h ^= b >> 2;
 393     }
 394     if (top3b & 04) {
 395         l ^= b << 63;
 396         h ^= b >> 1;
 397     }
 398
 399     *r1 = h;
 400     *r0 = l;
 401 }
 402 #endif
 403
 404 /*
 405  * Product of two polynomials a, b each with degree < 2 * BN_BITS2 - 1,
 406  * result is a polynomial r with degree < 4 * BN_BITS2 - 1 The caller MUST
 407  * ensure that the variables have the right amount of space allocated.
 408  */
 409 static void bn_GF2m_mul_2x2(BN_ULONG *r, const BN_ULONG a1, const BN_ULONG a0,
 410                             const BN_ULONG b1, const BN_ULONG b0)
 411 {
 412     BN_ULONG m1, m0;
 413     /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */
 414     bn_GF2m_mul_1x1(r + 3, r + 2, a1, b1);
 415     bn_GF2m_mul_1x1(r + 1, r, a0, b0);
 416     bn_GF2m_mul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1);
 417     /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */
 418     r[2] ^= m1 ^ r[1] ^ r[3];   /* h0 ^= m1 ^ l1 ^ h1; */
 419     r[1] = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0; /* l1 ^= l0 ^ h0 ^ m0; */
 420 }
 421
 422 /*
 423  * Add polynomials a and b and store result in r; r could be a or b, a and b
 424  * could be equal; r is the bitwise XOR of a and b.
 425  */
 426 int BN_GF2m_add(BIGNUM *r, const BIGNUM *a, const BIGNUM *b)
 427 {
 428     int i;
 429     const BIGNUM *at, *bt;
 430
 431     bn_check_top(a);
 432     bn_check_top(b);
 433
 434     if (a->top < b->top) {
 435         at = b;
 436         bt = a;
 437     } else {
 438         at = a;
 439         bt = b;
 440     }
 441
 442     if (bn_wexpand(r, at->top) == NULL)
 443         return 0;
 444
 445     for (i = 0; i < bt->top; i++) {
 446         r->d[i] = at->d[i] ^ bt->d[i];
 447     }
 448     for (; i < at->top; i++) {
 449         r->d[i] = at->d[i];
 450     }
 451
 452     r->top = at->top;
 453     bn_correct_top(r);
 454
 455     return 1;
 456 }
 457
 458 /*-
 459  * Some functions allow for representation of the irreducible polynomials
 460  * as an int[], say p.  The irreducible f(t) is then of the form:
 461  *     t^p[0] + t^p[1] + ... + t^p[k]
 462  * where m = p[0] > p[1] > ... > p[k] = 0.
 463  */
 464
 465 /* Performs modular reduction of a and store result in r.  r could be a. */
 466 int BN_GF2m_mod_arr(BIGNUM *r, const BIGNUM *a, const unsigned int p[])
 467 {
 468     int j, k;
 469     int n, dN, d0, d1;
 470     BN_ULONG zz, *z;
 471
 472     bn_check_top(a);
 473
 474     if (!p[0]) {
 475         /* reduction mod 1 => return 0 */
 476         BN_zero(r);
 477         return 1;
 478     }
 479
 480     /*
 481      * Since the algorithm does reduction in the r value, if a != r, copy the
 482      * contents of a into r so we can do reduction in r.
 483      */
 484     if (a != r) {
 485         if (!bn_wexpand(r, a->top))
 486             return 0;
 487         for (j = 0; j < a->top; j++) {
 488             r->d[j] = a->d[j];
 489         }
 490         r->top = a->top;
 491     }
 492     z = r->d;
 493
 494     /* start reduction */
 495     dN = p[0] / BN_BITS2;
 496     for (j = r->top - 1; j > dN;) {
 497         zz = z[j];
 498         if (z[j] == 0) {
 499             j--;
 500             continue;
 501         }
 502         z[j] = 0;
 503
 504         for (k = 1; p[k] != 0; k++) {
 505             /* reducing component t^p[k] */
 506             n = p[0] - p[k];
 507             d0 = n % BN_BITS2;
 508             d1 = BN_BITS2 - d0;
 509             n /= BN_BITS2;
 510             z[j - n] ^= (zz >> d0);
 511             if (d0)
 512                 z[j - n - 1] ^= (zz << d1);
 513         }
 514
 515         /* reducing component t^0 */
 516         n = dN;
 517         d0 = p[0] % BN_BITS2;
 518         d1 = BN_BITS2 - d0;
 519         z[j - n] ^= (zz >> d0);
 520         if (d0)
 521             z[j - n - 1] ^= (zz << d1);
 522     }
 523
 524     /* final round of reduction */
 525     while (j == dN) {
 526
 527         d0 = p[0] % BN_BITS2;
 528         zz = z[dN] >> d0;
 529         if (zz == 0)
 530             break;
 531         d1 = BN_BITS2 - d0;
 532
 533         /* clear up the top d1 bits */
 534         if (d0)
 535             z[dN] = (z[dN] << d1) >> d1;
 536         else
 537             z[dN] = 0;
 538         z[0] ^= zz;             /* reduction t^0 component */
 539
 540         for (k = 1; p[k] != 0; k++) {
 541             BN_ULONG tmp_ulong;
 542
 543             /* reducing component t^p[k] */
 544             n = p[k] / BN_BITS2;
 545             d0 = p[k] % BN_BITS2;
 546             d1 = BN_BITS2 - d0;
 547             z[n] ^= (zz << d0);
 548             tmp_ulong = zz >> d1;
 549             if (d0 && tmp_ulong)
 550                 z[n + 1] ^= tmp_ulong;
 551         }
 552
 553     }
 554
 555     bn_correct_top(r);
 556     return 1;
 557 }
 558
 559 /*
 560  * Performs modular reduction of a by p and store result in r.  r could be a.
 561  * This function calls down to the BN_GF2m_mod_arr implementation; this wrapper
 562  * function is only provided for convenience; for best performance, use the
 563  * BN_GF2m_mod_arr function.
 564  */
 565 int BN_GF2m_mod(BIGNUM *r, const BIGNUM *a, const BIGNUM *p)
 566 {
 567     int ret = 0;
 568     const int max = BN_num_bits(p);
 569     unsigned int *arr = NULL;
 570     bn_check_top(a);
 571     bn_check_top(p);
 572     if ((arr =
 573          (unsigned int *)OPENSSL_malloc(sizeof(unsigned int) * max)) == NULL)
 574         goto err;
 575     ret = BN_GF2m_poly2arr(p, arr, max);
 576     if (!ret || ret > max) {
 577         BNerr(BN_F_BN_GF2M_MOD, BN_R_INVALID_LENGTH);
 578         goto err;
 579     }
 580     ret = BN_GF2m_mod_arr(r, a, arr);
 581     bn_check_top(r);
 582  err:
 583     if (arr)
 584         OPENSSL_free(arr);
 585     return ret;
 586 }
 587
 588 /*
 589  * Compute the product of two polynomials a and b, reduce modulo p, and store
 590  * the result in r.  r could be a or b; a could be b.
 591  */
 592 int BN_GF2m_mod_mul_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
 593                         const unsigned int p[], BN_CTX *ctx)
 594 {
 595     int zlen, i, j, k, ret = 0;
 596     BIGNUM *s;
 597     BN_ULONG x1, x0, y1, y0, zz[4];
 598
 599     bn_check_top(a);
 600     bn_check_top(b);
 601
 602     if (a == b) {
 603         return BN_GF2m_mod_sqr_arr(r, a, p, ctx);
 604     }
 605
 606     BN_CTX_start(ctx);
 607     if ((s = BN_CTX_get(ctx)) == NULL)
 608         goto err;
 609
 610     zlen = a->top + b->top + 4;
 611     if (!bn_wexpand(s, zlen))
 612         goto err;
 613     s->top = zlen;
 614
 615     for (i = 0; i < zlen; i++)
 616         s->d[i] = 0;
 617
 618     for (j = 0; j < b->top; j += 2) {
 619         y0 = b->d[j];
 620         y1 = ((j + 1) == b->top) ? 0 : b->d[j + 1];
 621         for (i = 0; i < a->top; i += 2) {
 622             x0 = a->d[i];
 623             x1 = ((i + 1) == a->top) ? 0 : a->d[i + 1];
 624             bn_GF2m_mul_2x2(zz, x1, x0, y1, y0);
 625             for (k = 0; k < 4; k++)
 626                 s->d[i + j + k] ^= zz[k];
 627         }
 628     }
 629
 630     bn_correct_top(s);
 631     if (BN_GF2m_mod_arr(r, s, p))
 632         ret = 1;
 633     bn_check_top(r);
 634
 635  err:
 636     BN_CTX_end(ctx);
 637     return ret;
 638 }
 639
 640 /*
 641  * Compute the product of two polynomials a and b, reduce modulo p, and store
 642  * the result in r.  r could be a or b; a could equal b. This function calls
 643  * down to the BN_GF2m_mod_mul_arr implementation; this wrapper function is
 644  * only provided for convenience; for best performance, use the
 645  * BN_GF2m_mod_mul_arr function.
 646  */
 647 int BN_GF2m_mod_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
 648                     const BIGNUM *p, BN_CTX *ctx)
 649 {
 650     int ret = 0;
 651     const int max = BN_num_bits(p);
 652     unsigned int *arr = NULL;
 653     bn_check_top(a);
 654     bn_check_top(b);
 655     bn_check_top(p);
 656     if ((arr =
 657          (unsigned int *)OPENSSL_malloc(sizeof(unsigned int) * max)) == NULL)
 658         goto err;
 659     ret = BN_GF2m_poly2arr(p, arr, max);
 660     if (!ret || ret > max) {
 661         BNerr(BN_F_BN_GF2M_MOD_MUL, BN_R_INVALID_LENGTH);
 662         goto err;
 663     }
 664     ret = BN_GF2m_mod_mul_arr(r, a, b, arr, ctx);
 665     bn_check_top(r);
 666  err:
 667     if (arr)
 668         OPENSSL_free(arr);
 669     return ret;
 670 }
 671
 672 /* Square a, reduce the result mod p, and store it in a.  r could be a. */
 673 int BN_GF2m_mod_sqr_arr(BIGNUM *r, const BIGNUM *a, const unsigned int p[],
 674                         BN_CTX *ctx)
 675 {
 676     int i, ret = 0;
 677     BIGNUM *s;
 678
 679     bn_check_top(a);
 680     BN_CTX_start(ctx);
 681     if ((s = BN_CTX_get(ctx)) == NULL)
 682         return 0;
 683     if (!bn_wexpand(s, 2 * a->top))
 684         goto err;
 685
 686     for (i = a->top - 1; i >= 0; i--) {
 687         s->d[2 * i + 1] = SQR1(a->d[i]);
 688         s->d[2 * i] = SQR0(a->d[i]);
 689     }
 690
 691     s->top = 2 * a->top;
 692     bn_correct_top(s);
 693     if (!BN_GF2m_mod_arr(r, s, p))
 694         goto err;
 695     bn_check_top(r);
 696     ret = 1;
 697  err:
 698     BN_CTX_end(ctx);
 699     return ret;
 700 }
 701
 702 /*
 703  * Square a, reduce the result mod p, and store it in a.  r could be a. This
 704  * function calls down to the BN_GF2m_mod_sqr_arr implementation; this
 705  * wrapper function is only provided for convenience; for best performance,
 706  * use the BN_GF2m_mod_sqr_arr function.
 707  */
 708 int BN_GF2m_mod_sqr(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
 709 {
 710     int ret = 0;
 711     const int max = BN_num_bits(p);
 712     unsigned int *arr = NULL;
 713
 714     bn_check_top(a);
 715     bn_check_top(p);
 716     if ((arr =
 717          (unsigned int *)OPENSSL_malloc(sizeof(unsigned int) * max)) == NULL)
 718         goto err;
 719     ret = BN_GF2m_poly2arr(p, arr, max);
 720     if (!ret || ret > max) {
 721         BNerr(BN_F_BN_GF2M_MOD_SQR, BN_R_INVALID_LENGTH);
 722         goto err;
 723     }
 724     ret = BN_GF2m_mod_sqr_arr(r, a, arr, ctx);
 725     bn_check_top(r);
 726  err:
 727     if (arr)
 728         OPENSSL_free(arr);
 729     return ret;
 730 }
 731
 732 /*
 733  * Invert a, reduce modulo p, and store the result in r. r could be a. Uses
 734  * Modified Almost Inverse Algorithm (Algorithm 10) from Hankerson, D.,
 735  * Hernandez, J.L., and Menezes, A.  "Software Implementation of Elliptic
 736  * Curve Cryptography Over Binary Fields".
 737  */
 738 int BN_GF2m_mod_inv(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
 739 {
 740     BIGNUM *b, *c, *u, *v, *tmp;
 741     int ret = 0;
 742
 743     bn_check_top(a);
 744     bn_check_top(p);
 745
 746     BN_CTX_start(ctx);
 747
 748     b = BN_CTX_get(ctx);
 749     c = BN_CTX_get(ctx);
 750     u = BN_CTX_get(ctx);
 751     v = BN_CTX_get(ctx);
 752     if (v == NULL)
 753         goto err;
 754
 755     if (!BN_one(b))
 756         goto err;
 757     if (!BN_GF2m_mod(u, a, p))
 758         goto err;
 759     if (!BN_copy(v, p))
 760         goto err;
 761
 762     if (BN_is_zero(u))
 763         goto err;
 764
 765     while (1) {
 766         while (!BN_is_odd(u)) {
 767             if (BN_is_zero(u))
 768                 goto err;
 769             if (!BN_rshift1(u, u))
 770                 goto err;
 771             if (BN_is_odd(b)) {
 772                 if (!BN_GF2m_add(b, b, p))
 773                     goto err;
 774             }
 775             if (!BN_rshift1(b, b))
 776                 goto err;
 777         }
 778
 779         if (BN_abs_is_word(u, 1))
 780             break;
 781
 782         if (BN_num_bits(u) < BN_num_bits(v)) {
 783             tmp = u;
 784             u = v;
 785             v = tmp;
 786             tmp = b;
 787             b = c;
 788             c = tmp;
 789         }
 790
 791         if (!BN_GF2m_add(u, u, v))
 792             goto err;
 793         if (!BN_GF2m_add(b, b, c))
 794             goto err;
 795     }
 796
 797     if (!BN_copy(r, b))
 798         goto err;
 799     bn_check_top(r);
 800     ret = 1;
 801
 802  err:
 803     BN_CTX_end(ctx);
 804     return ret;
 805 }
 806
 807 /*
 808  * Invert xx, reduce modulo p, and store the result in r. r could be xx.
 809  * This function calls down to the BN_GF2m_mod_inv implementation; this
 810  * wrapper function is only provided for convenience; for best performance,
 811  * use the BN_GF2m_mod_inv function.
 812  */
 813 int BN_GF2m_mod_inv_arr(BIGNUM *r, const BIGNUM *xx, const unsigned int p[],
 814                         BN_CTX *ctx)
 815 {
 816     BIGNUM *field;
 817     int ret = 0;
 818
 819     bn_check_top(xx);
 820     BN_CTX_start(ctx);
 821     if ((field = BN_CTX_get(ctx)) == NULL)
 822         goto err;
 823     if (!BN_GF2m_arr2poly(p, field))
 824         goto err;
 825
 826     ret = BN_GF2m_mod_inv(r, xx, field, ctx);
 827     bn_check_top(r);
 828
 829  err:
 830     BN_CTX_end(ctx);
 831     return ret;
 832 }
 833
 834 #ifndef OPENSSL_SUN_GF2M_DIV
 835 /*
 836  * Divide y by x, reduce modulo p, and store the result in r. r could be x
 837  * or y, x could equal y.
 838  */
 839 int BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x,
 840                     const BIGNUM *p, BN_CTX *ctx)
 841 {
 842     BIGNUM *xinv = NULL;
 843     int ret = 0;
 844
 845     bn_check_top(y);
 846     bn_check_top(x);
 847     bn_check_top(p);
 848
 849     BN_CTX_start(ctx);
 850     xinv = BN_CTX_get(ctx);
 851     if (xinv == NULL)
 852         goto err;
 853
 854     if (!BN_GF2m_mod_inv(xinv, x, p, ctx))
 855         goto err;
 856     if (!BN_GF2m_mod_mul(r, y, xinv, p, ctx))
 857         goto err;
 858     bn_check_top(r);
 859     ret = 1;
 860
 861  err:
 862     BN_CTX_end(ctx);
 863     return ret;
 864 }
 865 #else
 866 /*
 867  * Divide y by x, reduce modulo p, and store the result in r. r could be x
 868  * or y, x could equal y. Uses algorithm Modular_Division_GF(2^m) from
 869  * Chang-Shantz, S.  "From Euclid's GCD to Montgomery Multiplication to the
 870  * Great Divide".
 871  */
 872 int BN_GF2m_mod_div(BIGNUM *r, const BIGNUM *y, const BIGNUM *x,
 873                     const BIGNUM *p, BN_CTX *ctx)
 874 {
 875     BIGNUM *a, *b, *u, *v;
 876     int ret = 0;
 877
 878     bn_check_top(y);
 879     bn_check_top(x);
 880     bn_check_top(p);
 881
 882     BN_CTX_start(ctx);
 883
 884     a = BN_CTX_get(ctx);
 885     b = BN_CTX_get(ctx);
 886     u = BN_CTX_get(ctx);
 887     v = BN_CTX_get(ctx);
 888     if (v == NULL)
 889         goto err;
 890
 891     /* reduce x and y mod p */
 892     if (!BN_GF2m_mod(u, y, p))
 893         goto err;
 894     if (!BN_GF2m_mod(a, x, p))
 895         goto err;
 896     if (!BN_copy(b, p))
 897         goto err;
 898
 899     while (!BN_is_odd(a)) {
 900         if (!BN_rshift1(a, a))
 901             goto err;
 902         if (BN_is_odd(u))
 903             if (!BN_GF2m_add(u, u, p))
 904                 goto err;
 905         if (!BN_rshift1(u, u))
 906             goto err;
 907     }
 908
 909     do {
 910         if (BN_GF2m_cmp(b, a) > 0) {
 911             if (!BN_GF2m_add(b, b, a))
 912                 goto err;
 913             if (!BN_GF2m_add(v, v, u))
 914                 goto err;
 915             do {
 916                 if (!BN_rshift1(b, b))
 917                     goto err;
 918                 if (BN_is_odd(v))
 919                     if (!BN_GF2m_add(v, v, p))
 920                         goto err;
 921                 if (!BN_rshift1(v, v))
 922                     goto err;
 923             } while (!BN_is_odd(b));
 924         } else if (BN_abs_is_word(a, 1))
 925             break;
 926         else {
 927             if (!BN_GF2m_add(a, a, b))
 928                 goto err;
 929             if (!BN_GF2m_add(u, u, v))
 930                 goto err;
 931             do {
 932                 if (!BN_rshift1(a, a))
 933                     goto err;
 934                 if (BN_is_odd(u))
 935                     if (!BN_GF2m_add(u, u, p))
 936                         goto err;
 937                 if (!BN_rshift1(u, u))
 938                     goto err;
 939             } while (!BN_is_odd(a));
 940         }
 941     } while (1);
 942
 943     if (!BN_copy(r, u))
 944         goto err;
 945     bn_check_top(r);
 946     ret = 1;
 947
 948  err:
 949     BN_CTX_end(ctx);
 950     return ret;
 951 }
 952 #endif
 953
 954 /*
 955  * Divide yy by xx, reduce modulo p, and store the result in r. r could be xx
 956  * * or yy, xx could equal yy. This function calls down to the
 957  * BN_GF2m_mod_div implementation; this wrapper function is only provided for
 958  * convenience; for best performance, use the BN_GF2m_mod_div function.
 959  */
 960 int BN_GF2m_mod_div_arr(BIGNUM *r, const BIGNUM *yy, const BIGNUM *xx,
 961                         const unsigned int p[], BN_CTX *ctx)
 962 {
 963     BIGNUM *field;
 964     int ret = 0;
 965
 966     bn_check_top(yy);
 967     bn_check_top(xx);
 968
 969     BN_CTX_start(ctx);
 970     if ((field = BN_CTX_get(ctx)) == NULL)
 971         goto err;
 972     if (!BN_GF2m_arr2poly(p, field))
 973         goto err;
 974
 975     ret = BN_GF2m_mod_div(r, yy, xx, field, ctx);
 976     bn_check_top(r);
 977
 978  err:
 979     BN_CTX_end(ctx);
 980     return ret;
 981 }
 982
 983 /*
 984  * Compute the bth power of a, reduce modulo p, and store the result in r.  r
 985  * could be a. Uses simple square-and-multiply algorithm A.5.1 from IEEE
 986  * P1363.
 987  */
 988 int BN_GF2m_mod_exp_arr(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
 989                         const unsigned int p[], BN_CTX *ctx)
 990 {
 991     int ret = 0, i, n;
 992     BIGNUM *u;
 993
 994     bn_check_top(a);
 995     bn_check_top(b);
 996
 997     if (BN_is_zero(b))
 998         return (BN_one(r));
 999
1000     if (BN_abs_is_word(b, 1))
1001         return (BN_copy(r, a) != NULL);
1002
1003     BN_CTX_start(ctx);
1004     if ((u = BN_CTX_get(ctx)) == NULL)
1005         goto err;
1006
1007     if (!BN_GF2m_mod_arr(u, a, p))
1008         goto err;
1009
1010     n = BN_num_bits(b) - 1;
1011     for (i = n - 1; i >= 0; i--) {
1012         if (!BN_GF2m_mod_sqr_arr(u, u, p, ctx))
1013             goto err;
1014         if (BN_is_bit_set(b, i)) {
1015             if (!BN_GF2m_mod_mul_arr(u, u, a, p, ctx))
1016                 goto err;
1017         }
1018     }
1019     if (!BN_copy(r, u))
1020         goto err;
1021     bn_check_top(r);
1022     ret = 1;
1023  err:
1024     BN_CTX_end(ctx);
1025     return ret;
1026 }
1027
1028 /*
1029  * Compute the bth power of a, reduce modulo p, and store the result in r.  r
1030  * could be a. This function calls down to the BN_GF2m_mod_exp_arr
1031  * implementation; this wrapper function is only provided for convenience;
1032  * for best performance, use the BN_GF2m_mod_exp_arr function.
1033  */
1034 int BN_GF2m_mod_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *b,
1035                     const BIGNUM *p, BN_CTX *ctx)
1036 {
1037     int ret = 0;
1038     const int max = BN_num_bits(p);
1039     unsigned int *arr = NULL;
1040     bn_check_top(a);
1041     bn_check_top(b);
1042     bn_check_top(p);
1043     if ((arr =
1044          (unsigned int *)OPENSSL_malloc(sizeof(unsigned int) * max)) == NULL)
1045         goto err;
1046     ret = BN_GF2m_poly2arr(p, arr, max);
1047     if (!ret || ret > max) {
1048         BNerr(BN_F_BN_GF2M_MOD_EXP, BN_R_INVALID_LENGTH);
1049         goto err;
1050     }
1051     ret = BN_GF2m_mod_exp_arr(r, a, b, arr, ctx);
1052     bn_check_top(r);
1053  err:
1054     if (arr)
1055         OPENSSL_free(arr);
1056     return ret;
1057 }
1058
1059 /*
1060  * Compute the square root of a, reduce modulo p, and store the result in r.
1061  * r could be a. Uses exponentiation as in algorithm A.4.1 from IEEE P1363.
1062  */
1063 int BN_GF2m_mod_sqrt_arr(BIGNUM *r, const BIGNUM *a, const unsigned int p[],
1064                          BN_CTX *ctx)
1065 {
1066     int ret = 0;
1067     BIGNUM *u;
1068
1069     bn_check_top(a);
1070
1071     if (!p[0]) {
1072         /* reduction mod 1 => return 0 */
1073         BN_zero(r);
1074         return 1;
1075     }
1076
1077     BN_CTX_start(ctx);
1078     if ((u = BN_CTX_get(ctx)) == NULL)
1079         goto err;
1080
1081     if (!BN_set_bit(u, p[0] - 1))
1082         goto err;
1083     ret = BN_GF2m_mod_exp_arr(r, a, u, p, ctx);
1084     bn_check_top(r);
1085
1086  err:
1087     BN_CTX_end(ctx);
1088     return ret;
1089 }
1090
1091 /*
1092  * Compute the square root of a, reduce modulo p, and store the result in r.
1093  * r could be a. This function calls down to the BN_GF2m_mod_sqrt_arr
1094  * implementation; this wrapper function is only provided for convenience;
1095  * for best performance, use the BN_GF2m_mod_sqrt_arr function.
1096  */
1097 int BN_GF2m_mod_sqrt(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
1098 {
1099     int ret = 0;
1100     const int max = BN_num_bits(p);
1101     unsigned int *arr = NULL;
1102     bn_check_top(a);
1103     bn_check_top(p);
1104     if ((arr =
1105          (unsigned int *)OPENSSL_malloc(sizeof(unsigned int) * max)) == NULL)
1106         goto err;
1107     ret = BN_GF2m_poly2arr(p, arr, max);
1108     if (!ret || ret > max) {
1109         BNerr(BN_F_BN_GF2M_MOD_SQRT, BN_R_INVALID_LENGTH);
1110         goto err;
1111     }
1112     ret = BN_GF2m_mod_sqrt_arr(r, a, arr, ctx);
1113     bn_check_top(r);
1114  err:
1115     if (arr)
1116         OPENSSL_free(arr);
1117     return ret;
1118 }
1119
1120 /*
1121  * Find r such that r^2 + r = a mod p.  r could be a. If no r exists returns
1122  * 0. Uses algorithms A.4.7 and A.4.6 from IEEE P1363.
1123  */
1124 int BN_GF2m_mod_solve_quad_arr(BIGNUM *r, const BIGNUM *a_,
1125                                const unsigned int p[], BN_CTX *ctx)
1126 {
1127     int ret = 0, count = 0;
1128     unsigned int j;
1129     BIGNUM *a, *z, *rho, *w, *w2, *tmp;
1130
1131     bn_check_top(a_);
1132
1133     if (!p[0]) {
1134         /* reduction mod 1 => return 0 */
1135         BN_zero(r);
1136         return 1;
1137     }
1138
1139     BN_CTX_start(ctx);
1140     a = BN_CTX_get(ctx);
1141     z = BN_CTX_get(ctx);
1142     w = BN_CTX_get(ctx);
1143     if (w == NULL)
1144         goto err;
1145
1146     if (!BN_GF2m_mod_arr(a, a_, p))
1147         goto err;
1148
1149     if (BN_is_zero(a)) {
1150         BN_zero(r);
1151         ret = 1;
1152         goto err;
1153     }
1154
1155     if (p[0] & 0x1) {           /* m is odd */
1156         /* compute half-trace of a */
1157         if (!BN_copy(z, a))
1158             goto err;
1159         for (j = 1; j <= (p[0] - 1) / 2; j++) {
1160             if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1161                 goto err;
1162             if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1163                 goto err;
1164             if (!BN_GF2m_add(z, z, a))
1165                 goto err;
1166         }
1167
1168     } else {                    /* m is even */
1169
1170         rho = BN_CTX_get(ctx);
1171         w2 = BN_CTX_get(ctx);
1172         tmp = BN_CTX_get(ctx);
1173         if (tmp == NULL)
1174             goto err;
1175         do {
1176             if (!BN_rand(rho, p[0], 0, 0))
1177                 goto err;
1178             if (!BN_GF2m_mod_arr(rho, rho, p))
1179                 goto err;
1180             BN_zero(z);
1181             if (!BN_copy(w, rho))
1182                 goto err;
1183             for (j = 1; j <= p[0] - 1; j++) {
1184                 if (!BN_GF2m_mod_sqr_arr(z, z, p, ctx))
1185                     goto err;
1186                 if (!BN_GF2m_mod_sqr_arr(w2, w, p, ctx))
1187                     goto err;
1188                 if (!BN_GF2m_mod_mul_arr(tmp, w2, a, p, ctx))
1189                     goto err;
1190                 if (!BN_GF2m_add(z, z, tmp))
1191                     goto err;
1192                 if (!BN_GF2m_add(w, w2, rho))
1193                     goto err;
1194             }
1195             count++;
1196         } while (BN_is_zero(w) && (count < MAX_ITERATIONS));
1197         if (BN_is_zero(w)) {
1198             BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD_ARR, BN_R_TOO_MANY_ITERATIONS);
1199             goto err;
1200         }
1201     }
1202
1203     if (!BN_GF2m_mod_sqr_arr(w, z, p, ctx))
1204         goto err;
1205     if (!BN_GF2m_add(w, z, w))
1206         goto err;
1207     if (BN_GF2m_cmp(w, a)) {
1208         BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD_ARR, BN_R_NO_SOLUTION);
1209         goto err;
1210     }
1211
1212     if (!BN_copy(r, z))
1213         goto err;
1214     bn_check_top(r);
1215
1216     ret = 1;
1217
1218  err:
1219     BN_CTX_end(ctx);
1220     return ret;
1221 }
1222
1223 /*
1224  * Find r such that r^2 + r = a mod p.  r could be a. If no r exists returns
1225  * 0. This function calls down to the BN_GF2m_mod_solve_quad_arr
1226  * implementation; this wrapper function is only provided for convenience;
1227  * for best performance, use the BN_GF2m_mod_solve_quad_arr function.
1228  */
1229 int BN_GF2m_mod_solve_quad(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
1230                            BN_CTX *ctx)
1231 {
1232     int ret = 0;
1233     const int max = BN_num_bits(p);
1234     unsigned int *arr = NULL;
1235     bn_check_top(a);
1236     bn_check_top(p);
1237     if ((arr = (unsigned int *)OPENSSL_malloc(sizeof(unsigned int) *
1238                                               max)) == NULL)
1239         goto err;
1240     ret = BN_GF2m_poly2arr(p, arr, max);
1241     if (!ret || ret > max) {
1242         BNerr(BN_F_BN_GF2M_MOD_SOLVE_QUAD, BN_R_INVALID_LENGTH);
1243         goto err;
1244     }
1245     ret = BN_GF2m_mod_solve_quad_arr(r, a, arr, ctx);
1246     bn_check_top(r);
1247  err:
1248     if (arr)
1249         OPENSSL_free(arr);
1250     return ret;
1251 }
1252
1253 /*
1254  * Convert the bit-string representation of a polynomial ( \sum_{i=0}^n a_i *
1255  * x^i , where a_0 is *not* zero) into an array of integers corresponding to
1256  * the bits with non-zero coefficient. Up to max elements of the array will
1257  * be filled.  Return value is total number of coefficients that would be
1258  * extracted if array was large enough.
1259  */
1260 int BN_GF2m_poly2arr(const BIGNUM *a, unsigned int p[], int max)
1261 {
1262     int i, j, k = 0;
1263     BN_ULONG mask;
1264
1265     if (BN_is_zero(a) || !BN_is_bit_set(a, 0))
1266         /*
1267          * a_0 == 0 => return error (the unsigned int array must be
1268          * terminated by 0)
1269          */
1270         return 0;
1271
1272     for (i = a->top - 1; i >= 0; i--) {
1273         if (!a->d[i])
1274             /* skip word if a->d[i] == 0 */
1275             continue;
1276         mask = BN_TBIT;
1277         for (j = BN_BITS2 - 1; j >= 0; j--) {
1278             if (a->d[i] & mask) {
1279                 if (k < max)
1280                     p[k] = BN_BITS2 * i + j;
1281                 k++;
1282             }
1283             mask >>= 1;
1284         }
1285     }
1286
1287     return k;
1288 }
1289
1290 /*
1291  * Convert the coefficient array representation of a polynomial to a
1292  * bit-string.  The array must be terminated by 0.
1293  */
1294 int BN_GF2m_arr2poly(const unsigned int p[], BIGNUM *a)
1295 {
1296     int i;
1297
1298     bn_check_top(a);
1299     BN_zero(a);
1300     for (i = 0; p[i] != 0; i++) {
1301         if (BN_set_bit(a, p[i]) == 0)
1302             return 0;
1303     }
1304     BN_set_bit(a, 0);
1305     bn_check_top(a);
1306
1307     return 1;
1308 }
1309
1310 /*
1311  * Constant-time conditional swap of a and b.
1312  * a and b are swapped if condition is not 0.  The code assumes that at most one bit of condition is set.
1313  * nwords is the number of words to swap.  The code assumes that at least nwords are allocated in both a and b,
1314  * and that no more than nwords are used by either a or b.
1315  * a and b cannot be the same number
1316  */
1317 void BN_consttime_swap(BN_ULONG condition, BIGNUM *a, BIGNUM *b, int nwords)
1318 {
1319     BN_ULONG t;
1320     int i;
1321
1322     bn_wcheck_size(a, nwords);
1323     bn_wcheck_size(b, nwords);
1324
1325     assert(a != b);
1326     assert((condition & (condition - 1)) == 0);
1327     assert(sizeof(BN_ULONG) >= sizeof(int));
1328
1329     condition = ((condition - 1) >> (BN_BITS2 - 1)) - 1;
1330
1331     t = (a->top ^ b->top) & condition;
1332     a->top ^= t;
1333     b->top ^= t;
1334
1335 #define BN_CONSTTIME_SWAP(ind) \
1336         do { \
1337                 t = (a->d[ind] ^ b->d[ind]) & condition; \
1338                 a->d[ind] ^= t; \
1339                 b->d[ind] ^= t; \
1340         } while (0)
1341
1342     switch (nwords) {
1343     default:
1344         for (i = 10; i < nwords; i++)
1345             BN_CONSTTIME_SWAP(i);
1346         /* Fallthrough */
1347     case 10:
1348         BN_CONSTTIME_SWAP(9);   /* Fallthrough */
1349     case 9:
1350         BN_CONSTTIME_SWAP(8);   /* Fallthrough */
1351     case 8:
1352         BN_CONSTTIME_SWAP(7);   /* Fallthrough */
1353     case 7:
1354         BN_CONSTTIME_SWAP(6);   /* Fallthrough */
1355     case 6:
1356         BN_CONSTTIME_SWAP(5);   /* Fallthrough */
1357     case 5:
1358         BN_CONSTTIME_SWAP(4);   /* Fallthrough */
1359     case 4:
1360         BN_CONSTTIME_SWAP(3);   /* Fallthrough */
1361     case 3:
1362         BN_CONSTTIME_SWAP(2);   /* Fallthrough */
1363     case 2:
1364         BN_CONSTTIME_SWAP(1);   /* Fallthrough */
1365     case 1:
1366         BN_CONSTTIME_SWAP(0);
1367     }
1368 #undef BN_CONSTTIME_SWAP
1369 }