crypto/bn/bn_sqrt.c

   1 /*
   2  * Written by Lenka Fibikova <fibikova@exp-math.uni-essen.de> and Bodo
   3  * Moeller for the OpenSSL project.
   4  */
   5 /* ====================================================================
   6  * Copyright (c) 1998-2000 The OpenSSL Project.  All rights reserved.
   7  *
   8  * Redistribution and use in source and binary forms, with or without
   9  * modification, are permitted provided that the following conditions
  10  * are met:
  11  *
  12  * 1. Redistributions of source code must retain the above copyright
  13  *    notice, this list of conditions and the following disclaimer.
  14  *
  15  * 2. Redistributions in binary form must reproduce the above copyright
  16  *    notice, this list of conditions and the following disclaimer in
  17  *    the documentation and/or other materials provided with the
  18  *    distribution.
  19  *
  20  * 3. All advertising materials mentioning features or use of this
  21  *    software must display the following acknowledgment:
  22  *    "This product includes software developed by the OpenSSL Project
  23  *    for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
  24  *
  25  * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
  26  *    endorse or promote products derived from this software without
  27  *    prior written permission. For written permission, please contact
  28  *    openssl-core@openssl.org.
  29  *
  30  * 5. Products derived from this software may not be called "OpenSSL"
  31  *    nor may "OpenSSL" appear in their names without prior written
  32  *    permission of the OpenSSL Project.
  33  *
  34  * 6. Redistributions of any form whatsoever must retain the following
  35  *    acknowledgment:
  36  *    "This product includes software developed by the OpenSSL Project
  37  *    for use in the OpenSSL Toolkit (http://www.openssl.org/)"
  38  *
  39  * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
  40  * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  41  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  42  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE OpenSSL PROJECT OR
  43  * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  44  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  45  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  46  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  47  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
  48  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  49  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
  50  * OF THE POSSIBILITY OF SUCH DAMAGE.
  51  * ====================================================================
  52  *
  53  * This product includes cryptographic software written by Eric Young
  54  * (eay@cryptsoft.com).  This product includes software written by Tim
  55  * Hudson (tjh@cryptsoft.com).
  56  *
  57  */
  58
  59 #include "internal/cryptlib.h"
  60 #include "bn_lcl.h"
  61
  62 BIGNUM *BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
  63 /*
  64  * Returns 'ret' such that ret^2 == a (mod p), using the Tonelli/Shanks
  65  * algorithm (cf. Henri Cohen, "A Course in Algebraic Computational Number
  66  * Theory", algorithm 1.5.1). 'p' must be prime!
  67  */
  68 {
  69     BIGNUM *ret = in;
  70     int err = 1;
  71     int r;
  72     BIGNUM *A, *b, *q, *t, *x, *y;
  73     int e, i, j;
  74
  75     if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) {
  76         if (BN_abs_is_word(p, 2)) {
  77             if (ret == NULL)
  78                 ret = BN_new();
  79             if (ret == NULL)
  80                 goto end;
  81             if (!BN_set_word(ret, BN_is_bit_set(a, 0))) {
  82                 if (ret != in)
  83                     BN_free(ret);
  84                 return NULL;
  85             }
  86             bn_check_top(ret);
  87             return ret;
  88         }
  89
  90         BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME);
  91         return (NULL);
  92     }
  93
  94     if (BN_is_zero(a) || BN_is_one(a)) {
  95         if (ret == NULL)
  96             ret = BN_new();
  97         if (ret == NULL)
  98             goto end;
  99         if (!BN_set_word(ret, BN_is_one(a))) {
 100             if (ret != in)
 101                 BN_free(ret);
 102             return NULL;
 103         }
 104         bn_check_top(ret);
 105         return ret;
 106     }
 107
 108     BN_CTX_start(ctx);
 109     A = BN_CTX_get(ctx);
 110     b = BN_CTX_get(ctx);
 111     q = BN_CTX_get(ctx);
 112     t = BN_CTX_get(ctx);
 113     x = BN_CTX_get(ctx);
 114     y = BN_CTX_get(ctx);
 115     if (y == NULL)
 116         goto end;
 117
 118     if (ret == NULL)
 119         ret = BN_new();
 120     if (ret == NULL)
 121         goto end;
 122
 123     /* A = a mod p */
 124     if (!BN_nnmod(A, a, p, ctx))
 125         goto end;
 126
 127     /* now write  |p| - 1  as  2^e*q  where  q  is odd */
 128     e = 1;
 129     while (!BN_is_bit_set(p, e))
 130         e++;
 131     /* we'll set  q  later (if needed) */
 132
 133     if (e == 1) {
 134         /*-
 135          * The easy case:  (|p|-1)/2  is odd, so 2 has an inverse
 136          * modulo  (|p|-1)/2,  and square roots can be computed
 137          * directly by modular exponentiation.
 138          * We have
 139          *     2 * (|p|+1)/4 == 1   (mod (|p|-1)/2),
 140          * so we can use exponent  (|p|+1)/4,  i.e.  (|p|-3)/4 + 1.
 141          */
 142         if (!BN_rshift(q, p, 2))
 143             goto end;
 144         q->neg = 0;
 145         if (!BN_add_word(q, 1))
 146             goto end;
 147         if (!BN_mod_exp(ret, A, q, p, ctx))
 148             goto end;
 149         err = 0;
 150         goto vrfy;
 151     }
 152
 153     if (e == 2) {
 154         /*-
 155          * |p| == 5  (mod 8)
 156          *
 157          * In this case  2  is always a non-square since
 158          * Legendre(2,p) = (-1)^((p^2-1)/8)  for any odd prime.
 159          * So if  a  really is a square, then  2*a  is a non-square.
 160          * Thus for
 161          *      b := (2*a)^((|p|-5)/8),
 162          *      i := (2*a)*b^2
 163          * we have
 164          *     i^2 = (2*a)^((1 + (|p|-5)/4)*2)
 165          *         = (2*a)^((p-1)/2)
 166          *         = -1;
 167          * so if we set
 168          *      x := a*b*(i-1),
 169          * then
 170          *     x^2 = a^2 * b^2 * (i^2 - 2*i + 1)
 171          *         = a^2 * b^2 * (-2*i)
 172          *         = a*(-i)*(2*a*b^2)
 173          *         = a*(-i)*i
 174          *         = a.
 175          *
 176          * (This is due to A.O.L. Atkin,
 177          * <URL: http://listserv.nodak.edu/scripts/wa.exe?A2=ind9211&L=nmbrthry&O=T&P=562>,
 178          * November 1992.)
 179          */
 180
 181         /* t := 2*a */
 182         if (!BN_mod_lshift1_quick(t, A, p))
 183             goto end;
 184
 185         /* b := (2*a)^((|p|-5)/8) */
 186         if (!BN_rshift(q, p, 3))
 187             goto end;
 188         q->neg = 0;
 189         if (!BN_mod_exp(b, t, q, p, ctx))
 190             goto end;
 191
 192         /* y := b^2 */
 193         if (!BN_mod_sqr(y, b, p, ctx))
 194             goto end;
 195
 196         /* t := (2*a)*b^2 - 1 */
 197         if (!BN_mod_mul(t, t, y, p, ctx))
 198             goto end;
 199         if (!BN_sub_word(t, 1))
 200             goto end;
 201
 202         /* x = a*b*t */
 203         if (!BN_mod_mul(x, A, b, p, ctx))
 204             goto end;
 205         if (!BN_mod_mul(x, x, t, p, ctx))
 206             goto end;
 207
 208         if (!BN_copy(ret, x))
 209             goto end;
 210         err = 0;
 211         goto vrfy;
 212     }
 213
 214     /*
 215      * e > 2, so we really have to use the Tonelli/Shanks algorithm. First,
 216      * find some y that is not a square.
 217      */
 218     if (!BN_copy(q, p))
 219         goto end;               /* use 'q' as temp */
 220     q->neg = 0;
 221     i = 2;
 222     do {
 223         /*
 224          * For efficiency, try small numbers first; if this fails, try random
 225          * numbers.
 226          */
 227         if (i < 22) {
 228             if (!BN_set_word(y, i))
 229                 goto end;
 230         } else {
 231             if (!BN_pseudo_rand(y, BN_num_bits(p), 0, 0))
 232                 goto end;
 233             if (BN_ucmp(y, p) >= 0) {
 234                 if (!(p->neg ? BN_add : BN_sub) (y, y, p))
 235                     goto end;
 236             }
 237             /* now 0 <= y < |p| */
 238             if (BN_is_zero(y))
 239                 if (!BN_set_word(y, i))
 240                     goto end;
 241         }
 242
 243         r = BN_kronecker(y, q, ctx); /* here 'q' is |p| */
 244         if (r < -1)
 245             goto end;
 246         if (r == 0) {
 247             /* m divides p */
 248             BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME);
 249             goto end;
 250         }
 251     }
 252     while (r == 1 && ++i < 82);
 253
 254     if (r != -1) {
 255         /*
 256          * Many rounds and still no non-square -- this is more likely a bug
 257          * than just bad luck. Even if p is not prime, we should have found
 258          * some y such that r == -1.
 259          */
 260         BNerr(BN_F_BN_MOD_SQRT, BN_R_TOO_MANY_ITERATIONS);
 261         goto end;
 262     }
 263
 264     /* Here's our actual 'q': */
 265     if (!BN_rshift(q, q, e))
 266         goto end;
 267
 268     /*
 269      * Now that we have some non-square, we can find an element of order 2^e
 270      * by computing its q'th power.
 271      */
 272     if (!BN_mod_exp(y, y, q, p, ctx))
 273         goto end;
 274     if (BN_is_one(y)) {
 275         BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME);
 276         goto end;
 277     }
 278
 279     /*-
 280      * Now we know that (if  p  is indeed prime) there is an integer
 281      * k,  0 <= k < 2^e,  such that
 282      *
 283      *      a^q * y^k == 1   (mod p).
 284      *
 285      * As  a^q  is a square and  y  is not,  k  must be even.
 286      * q+1  is even, too, so there is an element
 287      *
 288      *     X := a^((q+1)/2) * y^(k/2),
 289      *
 290      * and it satisfies
 291      *
 292      *     X^2 = a^q * a     * y^k
 293      *         = a,
 294      *
 295      * so it is the square root that we are looking for.
 296      */
 297
 298     /* t := (q-1)/2  (note that  q  is odd) */
 299     if (!BN_rshift1(t, q))
 300         goto end;
 301
 302     /* x := a^((q-1)/2) */
 303     if (BN_is_zero(t)) {        /* special case: p = 2^e + 1 */
 304         if (!BN_nnmod(t, A, p, ctx))
 305             goto end;
 306         if (BN_is_zero(t)) {
 307             /* special case: a == 0  (mod p) */
 308             BN_zero(ret);
 309             err = 0;
 310             goto end;
 311         } else if (!BN_one(x))
 312             goto end;
 313     } else {
 314         if (!BN_mod_exp(x, A, t, p, ctx))
 315             goto end;
 316         if (BN_is_zero(x)) {
 317             /* special case: a == 0  (mod p) */
 318             BN_zero(ret);
 319             err = 0;
 320             goto end;
 321         }
 322     }
 323
 324     /* b := a*x^2  (= a^q) */
 325     if (!BN_mod_sqr(b, x, p, ctx))
 326         goto end;
 327     if (!BN_mod_mul(b, b, A, p, ctx))
 328         goto end;
 329
 330     /* x := a*x    (= a^((q+1)/2)) */
 331     if (!BN_mod_mul(x, x, A, p, ctx))
 332         goto end;
 333
 334     while (1) {
 335         /*-
 336          * Now  b  is  a^q * y^k  for some even  k  (0 <= k < 2^E
 337          * where  E  refers to the original value of  e,  which we
 338          * don't keep in a variable),  and  x  is  a^((q+1)/2) * y^(k/2).
 339          *
 340          * We have  a*b = x^2,
 341          *    y^2^(e-1) = -1,
 342          *    b^2^(e-1) = 1.
 343          */
 344
 345         if (BN_is_one(b)) {
 346             if (!BN_copy(ret, x))
 347                 goto end;
 348             err = 0;
 349             goto vrfy;
 350         }
 351
 352         /* find smallest  i  such that  b^(2^i) = 1 */
 353         i = 1;
 354         if (!BN_mod_sqr(t, b, p, ctx))
 355             goto end;
 356         while (!BN_is_one(t)) {
 357             i++;
 358             if (i == e) {
 359                 BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE);
 360                 goto end;
 361             }
 362             if (!BN_mod_mul(t, t, t, p, ctx))
 363                 goto end;
 364         }
 365
 366         /* t := y^2^(e - i - 1) */
 367         if (!BN_copy(t, y))
 368             goto end;
 369         for (j = e - i - 1; j > 0; j--) {
 370             if (!BN_mod_sqr(t, t, p, ctx))
 371                 goto end;
 372         }
 373         if (!BN_mod_mul(y, t, t, p, ctx))
 374             goto end;
 375         if (!BN_mod_mul(x, x, t, p, ctx))
 376             goto end;
 377         if (!BN_mod_mul(b, b, y, p, ctx))
 378             goto end;
 379         e = i;
 380     }
 381
 382  vrfy:
 383     if (!err) {
 384         /*
 385          * verify the result -- the input might have been not a square (test
 386          * added in 0.9.8)
 387          */
 388
 389         if (!BN_mod_sqr(x, ret, p, ctx))
 390             err = 1;
 391
 392         if (!err && 0 != BN_cmp(x, A)) {
 393             BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE);
 394             err = 1;
 395         }
 396     }
 397
 398  end:
 399     if (err) {
 400         if (ret != in)
 401             BN_clear_free(ret);
 402         ret = NULL;
 403     }
 404     BN_CTX_end(ctx);
 405     bn_check_top(ret);
 406     return ret;
 407 }