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