1160814Ssimon/* crypto/bn/bn_sqrt.c */ 2280297Sjkim/* 3280297Sjkim * Written by Lenka Fibikova <fibikova@exp-math.uni-essen.de> and Bodo 4280297Sjkim * Moeller for the OpenSSL project. 5280297Sjkim */ 6109998Smarkm/* ==================================================================== 7109998Smarkm * Copyright (c) 1998-2000 The OpenSSL Project. All rights reserved. 8109998Smarkm * 9109998Smarkm * Redistribution and use in source and binary forms, with or without 10109998Smarkm * modification, are permitted provided that the following conditions 11109998Smarkm * are met: 12109998Smarkm * 13109998Smarkm * 1. Redistributions of source code must retain the above copyright 14280297Sjkim * notice, this list of conditions and the following disclaimer. 15109998Smarkm * 16109998Smarkm * 2. Redistributions in binary form must reproduce the above copyright 17109998Smarkm * notice, this list of conditions and the following disclaimer in 18109998Smarkm * the documentation and/or other materials provided with the 19109998Smarkm * distribution. 20109998Smarkm * 21109998Smarkm * 3. All advertising materials mentioning features or use of this 22109998Smarkm * software must display the following acknowledgment: 23109998Smarkm * "This product includes software developed by the OpenSSL Project 24109998Smarkm * for use in the OpenSSL Toolkit. (http://www.openssl.org/)" 25109998Smarkm * 26109998Smarkm * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to 27109998Smarkm * endorse or promote products derived from this software without 28109998Smarkm * prior written permission. For written permission, please contact 29109998Smarkm * openssl-core@openssl.org. 30109998Smarkm * 31109998Smarkm * 5. Products derived from this software may not be called "OpenSSL" 32109998Smarkm * nor may "OpenSSL" appear in their names without prior written 33109998Smarkm * permission of the OpenSSL Project. 34109998Smarkm * 35109998Smarkm * 6. Redistributions of any form whatsoever must retain the following 36109998Smarkm * acknowledgment: 37109998Smarkm * "This product includes software developed by the OpenSSL Project 38109998Smarkm * for use in the OpenSSL Toolkit (http://www.openssl.org/)" 39109998Smarkm * 40109998Smarkm * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY 41109998Smarkm * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 42109998Smarkm * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 43109998Smarkm * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE OpenSSL PROJECT OR 44109998Smarkm * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 45109998Smarkm * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 46109998Smarkm * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 47109998Smarkm * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 48109998Smarkm * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 49109998Smarkm * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 50109998Smarkm * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED 51109998Smarkm * OF THE POSSIBILITY OF SUCH DAMAGE. 52109998Smarkm * ==================================================================== 53109998Smarkm * 54109998Smarkm * This product includes cryptographic software written by Eric Young 55109998Smarkm * (eay@cryptsoft.com). This product includes software written by Tim 56109998Smarkm * Hudson (tjh@cryptsoft.com). 57109998Smarkm * 58109998Smarkm */ 59109998Smarkm 60109998Smarkm#include "cryptlib.h" 61109998Smarkm#include "bn_lcl.h" 62109998Smarkm 63280297SjkimBIGNUM *BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) 64280297Sjkim/* 65280297Sjkim * Returns 'ret' such that ret^2 == a (mod p), using the Tonelli/Shanks 66280297Sjkim * algorithm (cf. Henri Cohen, "A Course in Algebraic Computational Number 67280297Sjkim * Theory", algorithm 1.5.1). 'p' must be prime! 68109998Smarkm */ 69280297Sjkim{ 70280297Sjkim BIGNUM *ret = in; 71280297Sjkim int err = 1; 72280297Sjkim int r; 73280297Sjkim BIGNUM *A, *b, *q, *t, *x, *y; 74280297Sjkim int e, i, j; 75109998Smarkm 76280297Sjkim if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) { 77280297Sjkim if (BN_abs_is_word(p, 2)) { 78280297Sjkim if (ret == NULL) 79280297Sjkim ret = BN_new(); 80280297Sjkim if (ret == NULL) 81280297Sjkim goto end; 82280297Sjkim if (!BN_set_word(ret, BN_is_bit_set(a, 0))) { 83280297Sjkim if (ret != in) 84280297Sjkim BN_free(ret); 85280297Sjkim return NULL; 86280297Sjkim } 87280297Sjkim bn_check_top(ret); 88280297Sjkim return ret; 89280297Sjkim } 90109998Smarkm 91280297Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 92280297Sjkim return (NULL); 93280297Sjkim } 94109998Smarkm 95280297Sjkim if (BN_is_zero(a) || BN_is_one(a)) { 96280297Sjkim if (ret == NULL) 97280297Sjkim ret = BN_new(); 98280297Sjkim if (ret == NULL) 99280297Sjkim goto end; 100280297Sjkim if (!BN_set_word(ret, BN_is_one(a))) { 101280297Sjkim if (ret != in) 102280297Sjkim BN_free(ret); 103280297Sjkim return NULL; 104280297Sjkim } 105280297Sjkim bn_check_top(ret); 106280297Sjkim return ret; 107280297Sjkim } 108109998Smarkm 109280297Sjkim BN_CTX_start(ctx); 110280297Sjkim A = BN_CTX_get(ctx); 111280297Sjkim b = BN_CTX_get(ctx); 112280297Sjkim q = BN_CTX_get(ctx); 113280297Sjkim t = BN_CTX_get(ctx); 114280297Sjkim x = BN_CTX_get(ctx); 115280297Sjkim y = BN_CTX_get(ctx); 116280297Sjkim if (y == NULL) 117280297Sjkim goto end; 118160814Ssimon 119280297Sjkim if (ret == NULL) 120280297Sjkim ret = BN_new(); 121280297Sjkim if (ret == NULL) 122280297Sjkim goto end; 123109998Smarkm 124280297Sjkim /* A = a mod p */ 125280297Sjkim if (!BN_nnmod(A, a, p, ctx)) 126280297Sjkim goto end; 127109998Smarkm 128280297Sjkim /* now write |p| - 1 as 2^e*q where q is odd */ 129280297Sjkim e = 1; 130280297Sjkim while (!BN_is_bit_set(p, e)) 131280297Sjkim e++; 132280297Sjkim /* we'll set q later (if needed) */ 133109998Smarkm 134280297Sjkim if (e == 1) { 135280297Sjkim /*- 136280297Sjkim * The easy case: (|p|-1)/2 is odd, so 2 has an inverse 137280297Sjkim * modulo (|p|-1)/2, and square roots can be computed 138280297Sjkim * directly by modular exponentiation. 139280297Sjkim * We have 140280297Sjkim * 2 * (|p|+1)/4 == 1 (mod (|p|-1)/2), 141280297Sjkim * so we can use exponent (|p|+1)/4, i.e. (|p|-3)/4 + 1. 142280297Sjkim */ 143280297Sjkim if (!BN_rshift(q, p, 2)) 144280297Sjkim goto end; 145280297Sjkim q->neg = 0; 146280297Sjkim if (!BN_add_word(q, 1)) 147280297Sjkim goto end; 148280297Sjkim if (!BN_mod_exp(ret, A, q, p, ctx)) 149280297Sjkim goto end; 150280297Sjkim err = 0; 151280297Sjkim goto vrfy; 152280297Sjkim } 153109998Smarkm 154280297Sjkim if (e == 2) { 155280297Sjkim /*- 156280297Sjkim * |p| == 5 (mod 8) 157280297Sjkim * 158280297Sjkim * In this case 2 is always a non-square since 159280297Sjkim * Legendre(2,p) = (-1)^((p^2-1)/8) for any odd prime. 160280297Sjkim * So if a really is a square, then 2*a is a non-square. 161280297Sjkim * Thus for 162280297Sjkim * b := (2*a)^((|p|-5)/8), 163280297Sjkim * i := (2*a)*b^2 164280297Sjkim * we have 165280297Sjkim * i^2 = (2*a)^((1 + (|p|-5)/4)*2) 166280297Sjkim * = (2*a)^((p-1)/2) 167280297Sjkim * = -1; 168280297Sjkim * so if we set 169280297Sjkim * x := a*b*(i-1), 170280297Sjkim * then 171280297Sjkim * x^2 = a^2 * b^2 * (i^2 - 2*i + 1) 172280297Sjkim * = a^2 * b^2 * (-2*i) 173280297Sjkim * = a*(-i)*(2*a*b^2) 174280297Sjkim * = a*(-i)*i 175280297Sjkim * = a. 176280297Sjkim * 177280297Sjkim * (This is due to A.O.L. Atkin, 178280297Sjkim * <URL: http://listserv.nodak.edu/scripts/wa.exe?A2=ind9211&L=nmbrthry&O=T&P=562>, 179280297Sjkim * November 1992.) 180280297Sjkim */ 181109998Smarkm 182280297Sjkim /* t := 2*a */ 183280297Sjkim if (!BN_mod_lshift1_quick(t, A, p)) 184280297Sjkim goto end; 185109998Smarkm 186280297Sjkim /* b := (2*a)^((|p|-5)/8) */ 187280297Sjkim if (!BN_rshift(q, p, 3)) 188280297Sjkim goto end; 189280297Sjkim q->neg = 0; 190280297Sjkim if (!BN_mod_exp(b, t, q, p, ctx)) 191280297Sjkim goto end; 192109998Smarkm 193280297Sjkim /* y := b^2 */ 194280297Sjkim if (!BN_mod_sqr(y, b, p, ctx)) 195280297Sjkim goto end; 196109998Smarkm 197280297Sjkim /* t := (2*a)*b^2 - 1 */ 198280297Sjkim if (!BN_mod_mul(t, t, y, p, ctx)) 199280297Sjkim goto end; 200280297Sjkim if (!BN_sub_word(t, 1)) 201280297Sjkim goto end; 202109998Smarkm 203280297Sjkim /* x = a*b*t */ 204280297Sjkim if (!BN_mod_mul(x, A, b, p, ctx)) 205280297Sjkim goto end; 206280297Sjkim if (!BN_mod_mul(x, x, t, p, ctx)) 207280297Sjkim goto end; 208109998Smarkm 209280297Sjkim if (!BN_copy(ret, x)) 210280297Sjkim goto end; 211280297Sjkim err = 0; 212280297Sjkim goto vrfy; 213280297Sjkim } 214109998Smarkm 215280297Sjkim /* 216280297Sjkim * e > 2, so we really have to use the Tonelli/Shanks algorithm. First, 217280297Sjkim * find some y that is not a square. 218280297Sjkim */ 219280297Sjkim if (!BN_copy(q, p)) 220280297Sjkim goto end; /* use 'q' as temp */ 221280297Sjkim q->neg = 0; 222280297Sjkim i = 2; 223280297Sjkim do { 224280297Sjkim /* 225280297Sjkim * For efficiency, try small numbers first; if this fails, try random 226280297Sjkim * numbers. 227280297Sjkim */ 228280297Sjkim if (i < 22) { 229280297Sjkim if (!BN_set_word(y, i)) 230280297Sjkim goto end; 231280297Sjkim } else { 232280297Sjkim if (!BN_pseudo_rand(y, BN_num_bits(p), 0, 0)) 233280297Sjkim goto end; 234280297Sjkim if (BN_ucmp(y, p) >= 0) { 235280297Sjkim if (!(p->neg ? BN_add : BN_sub) (y, y, p)) 236280297Sjkim goto end; 237280297Sjkim } 238280297Sjkim /* now 0 <= y < |p| */ 239280297Sjkim if (BN_is_zero(y)) 240280297Sjkim if (!BN_set_word(y, i)) 241280297Sjkim goto end; 242280297Sjkim } 243109998Smarkm 244280297Sjkim r = BN_kronecker(y, q, ctx); /* here 'q' is |p| */ 245280297Sjkim if (r < -1) 246280297Sjkim goto end; 247280297Sjkim if (r == 0) { 248280297Sjkim /* m divides p */ 249280297Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 250280297Sjkim goto end; 251280297Sjkim } 252280297Sjkim } 253280297Sjkim while (r == 1 && ++i < 82); 254109998Smarkm 255280297Sjkim if (r != -1) { 256280297Sjkim /* 257280297Sjkim * Many rounds and still no non-square -- this is more likely a bug 258280297Sjkim * than just bad luck. Even if p is not prime, we should have found 259280297Sjkim * some y such that r == -1. 260280297Sjkim */ 261280297Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_TOO_MANY_ITERATIONS); 262280297Sjkim goto end; 263280297Sjkim } 264109998Smarkm 265280297Sjkim /* Here's our actual 'q': */ 266280297Sjkim if (!BN_rshift(q, q, e)) 267280297Sjkim goto end; 268109998Smarkm 269280297Sjkim /* 270280297Sjkim * Now that we have some non-square, we can find an element of order 2^e 271280297Sjkim * by computing its q'th power. 272280297Sjkim */ 273280297Sjkim if (!BN_mod_exp(y, y, q, p, ctx)) 274280297Sjkim goto end; 275280297Sjkim if (BN_is_one(y)) { 276280297Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 277280297Sjkim goto end; 278280297Sjkim } 279109998Smarkm 280280297Sjkim /*- 281280297Sjkim * Now we know that (if p is indeed prime) there is an integer 282280297Sjkim * k, 0 <= k < 2^e, such that 283280297Sjkim * 284280297Sjkim * a^q * y^k == 1 (mod p). 285280297Sjkim * 286280297Sjkim * As a^q is a square and y is not, k must be even. 287280297Sjkim * q+1 is even, too, so there is an element 288280297Sjkim * 289280297Sjkim * X := a^((q+1)/2) * y^(k/2), 290280297Sjkim * 291280297Sjkim * and it satisfies 292280297Sjkim * 293280297Sjkim * X^2 = a^q * a * y^k 294280297Sjkim * = a, 295280297Sjkim * 296280297Sjkim * so it is the square root that we are looking for. 297280297Sjkim */ 298109998Smarkm 299280297Sjkim /* t := (q-1)/2 (note that q is odd) */ 300280297Sjkim if (!BN_rshift1(t, q)) 301280297Sjkim goto end; 302280297Sjkim 303280297Sjkim /* x := a^((q-1)/2) */ 304280297Sjkim if (BN_is_zero(t)) { /* special case: p = 2^e + 1 */ 305280297Sjkim if (!BN_nnmod(t, A, p, ctx)) 306280297Sjkim goto end; 307280297Sjkim if (BN_is_zero(t)) { 308280297Sjkim /* special case: a == 0 (mod p) */ 309280297Sjkim BN_zero(ret); 310280297Sjkim err = 0; 311280297Sjkim goto end; 312280297Sjkim } else if (!BN_one(x)) 313280297Sjkim goto end; 314280297Sjkim } else { 315280297Sjkim if (!BN_mod_exp(x, A, t, p, ctx)) 316280297Sjkim goto end; 317280297Sjkim if (BN_is_zero(x)) { 318280297Sjkim /* special case: a == 0 (mod p) */ 319280297Sjkim BN_zero(ret); 320280297Sjkim err = 0; 321280297Sjkim goto end; 322280297Sjkim } 323280297Sjkim } 324280297Sjkim 325280297Sjkim /* b := a*x^2 (= a^q) */ 326280297Sjkim if (!BN_mod_sqr(b, x, p, ctx)) 327280297Sjkim goto end; 328280297Sjkim if (!BN_mod_mul(b, b, A, p, ctx)) 329280297Sjkim goto end; 330280297Sjkim 331280297Sjkim /* x := a*x (= a^((q+1)/2)) */ 332280297Sjkim if (!BN_mod_mul(x, x, A, p, ctx)) 333280297Sjkim goto end; 334280297Sjkim 335280297Sjkim while (1) { 336280297Sjkim /*- 337280297Sjkim * Now b is a^q * y^k for some even k (0 <= k < 2^E 338280297Sjkim * where E refers to the original value of e, which we 339280297Sjkim * don't keep in a variable), and x is a^((q+1)/2) * y^(k/2). 340280297Sjkim * 341280297Sjkim * We have a*b = x^2, 342280297Sjkim * y^2^(e-1) = -1, 343280297Sjkim * b^2^(e-1) = 1. 344280297Sjkim */ 345280297Sjkim 346280297Sjkim if (BN_is_one(b)) { 347280297Sjkim if (!BN_copy(ret, x)) 348280297Sjkim goto end; 349280297Sjkim err = 0; 350280297Sjkim goto vrfy; 351280297Sjkim } 352280297Sjkim 353280297Sjkim /* find smallest i such that b^(2^i) = 1 */ 354280297Sjkim i = 1; 355280297Sjkim if (!BN_mod_sqr(t, b, p, ctx)) 356280297Sjkim goto end; 357280297Sjkim while (!BN_is_one(t)) { 358280297Sjkim i++; 359280297Sjkim if (i == e) { 360280297Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE); 361280297Sjkim goto end; 362280297Sjkim } 363280297Sjkim if (!BN_mod_mul(t, t, t, p, ctx)) 364280297Sjkim goto end; 365280297Sjkim } 366280297Sjkim 367280297Sjkim /* t := y^2^(e - i - 1) */ 368280297Sjkim if (!BN_copy(t, y)) 369280297Sjkim goto end; 370280297Sjkim for (j = e - i - 1; j > 0; j--) { 371280297Sjkim if (!BN_mod_sqr(t, t, p, ctx)) 372280297Sjkim goto end; 373280297Sjkim } 374280297Sjkim if (!BN_mod_mul(y, t, t, p, ctx)) 375280297Sjkim goto end; 376280297Sjkim if (!BN_mod_mul(x, x, t, p, ctx)) 377280297Sjkim goto end; 378280297Sjkim if (!BN_mod_mul(b, b, y, p, ctx)) 379280297Sjkim goto end; 380280297Sjkim e = i; 381280297Sjkim } 382280297Sjkim 383160814Ssimon vrfy: 384280297Sjkim if (!err) { 385280297Sjkim /* 386280297Sjkim * verify the result -- the input might have been not a square (test 387280297Sjkim * added in 0.9.8) 388280297Sjkim */ 389160814Ssimon 390280297Sjkim if (!BN_mod_sqr(x, ret, p, ctx)) 391280297Sjkim err = 1; 392280297Sjkim 393280297Sjkim if (!err && 0 != BN_cmp(x, A)) { 394280297Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE); 395280297Sjkim err = 1; 396280297Sjkim } 397280297Sjkim } 398280297Sjkim 399109998Smarkm end: 400280297Sjkim if (err) { 401280297Sjkim if (ret != NULL && ret != in) { 402280297Sjkim BN_clear_free(ret); 403280297Sjkim } 404280297Sjkim ret = NULL; 405280297Sjkim } 406280297Sjkim BN_CTX_end(ctx); 407280297Sjkim bn_check_top(ret); 408280297Sjkim return ret; 409280297Sjkim} 410