crypto/bn/bn_sqrt.c

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