1160814Ssimon/* crypto/bn/bn_sqrt.c */ 2280304Sjkim/* 3280304Sjkim * Written by Lenka Fibikova <fibikova@exp-math.uni-essen.de> and Bodo 4280304Sjkim * Moeller for the OpenSSL project. 5280304Sjkim */ 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 14280304Sjkim * 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 63280304SjkimBIGNUM *BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) 64280304Sjkim/* 65280304Sjkim * Returns 'ret' such that ret^2 == a (mod p), using the Tonelli/Shanks 66280304Sjkim * algorithm (cf. Henri Cohen, "A Course in Algebraic Computational Number 67280304Sjkim * Theory", algorithm 1.5.1). 'p' must be prime! 68109998Smarkm */ 69280304Sjkim{ 70280304Sjkim BIGNUM *ret = in; 71280304Sjkim int err = 1; 72280304Sjkim int r; 73280304Sjkim BIGNUM *A, *b, *q, *t, *x, *y; 74280304Sjkim int e, i, j; 75109998Smarkm 76280304Sjkim if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) { 77280304Sjkim if (BN_abs_is_word(p, 2)) { 78280304Sjkim if (ret == NULL) 79280304Sjkim ret = BN_new(); 80280304Sjkim if (ret == NULL) 81280304Sjkim goto end; 82280304Sjkim if (!BN_set_word(ret, BN_is_bit_set(a, 0))) { 83280304Sjkim if (ret != in) 84280304Sjkim BN_free(ret); 85280304Sjkim return NULL; 86280304Sjkim } 87280304Sjkim bn_check_top(ret); 88280304Sjkim return ret; 89280304Sjkim } 90109998Smarkm 91280304Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 92280304Sjkim return (NULL); 93280304Sjkim } 94109998Smarkm 95280304Sjkim if (BN_is_zero(a) || BN_is_one(a)) { 96280304Sjkim if (ret == NULL) 97280304Sjkim ret = BN_new(); 98280304Sjkim if (ret == NULL) 99280304Sjkim goto end; 100280304Sjkim if (!BN_set_word(ret, BN_is_one(a))) { 101280304Sjkim if (ret != in) 102280304Sjkim BN_free(ret); 103280304Sjkim return NULL; 104280304Sjkim } 105280304Sjkim bn_check_top(ret); 106280304Sjkim return ret; 107280304Sjkim } 108109998Smarkm 109280304Sjkim BN_CTX_start(ctx); 110280304Sjkim A = BN_CTX_get(ctx); 111280304Sjkim b = BN_CTX_get(ctx); 112280304Sjkim q = BN_CTX_get(ctx); 113280304Sjkim t = BN_CTX_get(ctx); 114280304Sjkim x = BN_CTX_get(ctx); 115280304Sjkim y = BN_CTX_get(ctx); 116280304Sjkim if (y == NULL) 117280304Sjkim goto end; 118160814Ssimon 119280304Sjkim if (ret == NULL) 120280304Sjkim ret = BN_new(); 121280304Sjkim if (ret == NULL) 122280304Sjkim goto end; 123109998Smarkm 124280304Sjkim /* A = a mod p */ 125280304Sjkim if (!BN_nnmod(A, a, p, ctx)) 126280304Sjkim goto end; 127109998Smarkm 128280304Sjkim /* now write |p| - 1 as 2^e*q where q is odd */ 129280304Sjkim e = 1; 130280304Sjkim while (!BN_is_bit_set(p, e)) 131280304Sjkim e++; 132280304Sjkim /* we'll set q later (if needed) */ 133109998Smarkm 134280304Sjkim if (e == 1) { 135280304Sjkim /*- 136280304Sjkim * The easy case: (|p|-1)/2 is odd, so 2 has an inverse 137280304Sjkim * modulo (|p|-1)/2, and square roots can be computed 138280304Sjkim * directly by modular exponentiation. 139280304Sjkim * We have 140280304Sjkim * 2 * (|p|+1)/4 == 1 (mod (|p|-1)/2), 141280304Sjkim * so we can use exponent (|p|+1)/4, i.e. (|p|-3)/4 + 1. 142280304Sjkim */ 143280304Sjkim if (!BN_rshift(q, p, 2)) 144280304Sjkim goto end; 145280304Sjkim q->neg = 0; 146280304Sjkim if (!BN_add_word(q, 1)) 147280304Sjkim goto end; 148280304Sjkim if (!BN_mod_exp(ret, A, q, p, ctx)) 149280304Sjkim goto end; 150280304Sjkim err = 0; 151280304Sjkim goto vrfy; 152280304Sjkim } 153109998Smarkm 154280304Sjkim if (e == 2) { 155280304Sjkim /*- 156280304Sjkim * |p| == 5 (mod 8) 157280304Sjkim * 158280304Sjkim * In this case 2 is always a non-square since 159280304Sjkim * Legendre(2,p) = (-1)^((p^2-1)/8) for any odd prime. 160280304Sjkim * So if a really is a square, then 2*a is a non-square. 161280304Sjkim * Thus for 162280304Sjkim * b := (2*a)^((|p|-5)/8), 163280304Sjkim * i := (2*a)*b^2 164280304Sjkim * we have 165280304Sjkim * i^2 = (2*a)^((1 + (|p|-5)/4)*2) 166280304Sjkim * = (2*a)^((p-1)/2) 167280304Sjkim * = -1; 168280304Sjkim * so if we set 169280304Sjkim * x := a*b*(i-1), 170280304Sjkim * then 171280304Sjkim * x^2 = a^2 * b^2 * (i^2 - 2*i + 1) 172280304Sjkim * = a^2 * b^2 * (-2*i) 173280304Sjkim * = a*(-i)*(2*a*b^2) 174280304Sjkim * = a*(-i)*i 175280304Sjkim * = a. 176280304Sjkim * 177280304Sjkim * (This is due to A.O.L. Atkin, 178280304Sjkim * <URL: http://listserv.nodak.edu/scripts/wa.exe?A2=ind9211&L=nmbrthry&O=T&P=562>, 179280304Sjkim * November 1992.) 180280304Sjkim */ 181109998Smarkm 182280304Sjkim /* t := 2*a */ 183280304Sjkim if (!BN_mod_lshift1_quick(t, A, p)) 184280304Sjkim goto end; 185109998Smarkm 186280304Sjkim /* b := (2*a)^((|p|-5)/8) */ 187280304Sjkim if (!BN_rshift(q, p, 3)) 188280304Sjkim goto end; 189280304Sjkim q->neg = 0; 190280304Sjkim if (!BN_mod_exp(b, t, q, p, ctx)) 191280304Sjkim goto end; 192109998Smarkm 193280304Sjkim /* y := b^2 */ 194280304Sjkim if (!BN_mod_sqr(y, b, p, ctx)) 195280304Sjkim goto end; 196109998Smarkm 197280304Sjkim /* t := (2*a)*b^2 - 1 */ 198280304Sjkim if (!BN_mod_mul(t, t, y, p, ctx)) 199280304Sjkim goto end; 200280304Sjkim if (!BN_sub_word(t, 1)) 201280304Sjkim goto end; 202109998Smarkm 203280304Sjkim /* x = a*b*t */ 204280304Sjkim if (!BN_mod_mul(x, A, b, p, ctx)) 205280304Sjkim goto end; 206280304Sjkim if (!BN_mod_mul(x, x, t, p, ctx)) 207280304Sjkim goto end; 208109998Smarkm 209280304Sjkim if (!BN_copy(ret, x)) 210280304Sjkim goto end; 211280304Sjkim err = 0; 212280304Sjkim goto vrfy; 213280304Sjkim } 214109998Smarkm 215280304Sjkim /* 216280304Sjkim * e > 2, so we really have to use the Tonelli/Shanks algorithm. First, 217280304Sjkim * find some y that is not a square. 218280304Sjkim */ 219280304Sjkim if (!BN_copy(q, p)) 220280304Sjkim goto end; /* use 'q' as temp */ 221280304Sjkim q->neg = 0; 222280304Sjkim i = 2; 223280304Sjkim do { 224280304Sjkim /* 225280304Sjkim * For efficiency, try small numbers first; if this fails, try random 226280304Sjkim * numbers. 227280304Sjkim */ 228280304Sjkim if (i < 22) { 229280304Sjkim if (!BN_set_word(y, i)) 230280304Sjkim goto end; 231280304Sjkim } else { 232280304Sjkim if (!BN_pseudo_rand(y, BN_num_bits(p), 0, 0)) 233280304Sjkim goto end; 234280304Sjkim if (BN_ucmp(y, p) >= 0) { 235280304Sjkim if (!(p->neg ? BN_add : BN_sub) (y, y, p)) 236280304Sjkim goto end; 237280304Sjkim } 238280304Sjkim /* now 0 <= y < |p| */ 239280304Sjkim if (BN_is_zero(y)) 240280304Sjkim if (!BN_set_word(y, i)) 241280304Sjkim goto end; 242280304Sjkim } 243109998Smarkm 244280304Sjkim r = BN_kronecker(y, q, ctx); /* here 'q' is |p| */ 245280304Sjkim if (r < -1) 246280304Sjkim goto end; 247280304Sjkim if (r == 0) { 248280304Sjkim /* m divides p */ 249280304Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 250280304Sjkim goto end; 251280304Sjkim } 252280304Sjkim } 253280304Sjkim while (r == 1 && ++i < 82); 254109998Smarkm 255280304Sjkim if (r != -1) { 256280304Sjkim /* 257280304Sjkim * Many rounds and still no non-square -- this is more likely a bug 258280304Sjkim * than just bad luck. Even if p is not prime, we should have found 259280304Sjkim * some y such that r == -1. 260280304Sjkim */ 261280304Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_TOO_MANY_ITERATIONS); 262280304Sjkim goto end; 263280304Sjkim } 264109998Smarkm 265280304Sjkim /* Here's our actual 'q': */ 266280304Sjkim if (!BN_rshift(q, q, e)) 267280304Sjkim goto end; 268109998Smarkm 269280304Sjkim /* 270280304Sjkim * Now that we have some non-square, we can find an element of order 2^e 271280304Sjkim * by computing its q'th power. 272280304Sjkim */ 273280304Sjkim if (!BN_mod_exp(y, y, q, p, ctx)) 274280304Sjkim goto end; 275280304Sjkim if (BN_is_one(y)) { 276280304Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 277280304Sjkim goto end; 278280304Sjkim } 279109998Smarkm 280280304Sjkim /*- 281280304Sjkim * Now we know that (if p is indeed prime) there is an integer 282280304Sjkim * k, 0 <= k < 2^e, such that 283280304Sjkim * 284280304Sjkim * a^q * y^k == 1 (mod p). 285280304Sjkim * 286280304Sjkim * As a^q is a square and y is not, k must be even. 287280304Sjkim * q+1 is even, too, so there is an element 288280304Sjkim * 289280304Sjkim * X := a^((q+1)/2) * y^(k/2), 290280304Sjkim * 291280304Sjkim * and it satisfies 292280304Sjkim * 293280304Sjkim * X^2 = a^q * a * y^k 294280304Sjkim * = a, 295280304Sjkim * 296280304Sjkim * so it is the square root that we are looking for. 297280304Sjkim */ 298109998Smarkm 299280304Sjkim /* t := (q-1)/2 (note that q is odd) */ 300280304Sjkim if (!BN_rshift1(t, q)) 301280304Sjkim goto end; 302280304Sjkim 303280304Sjkim /* x := a^((q-1)/2) */ 304280304Sjkim if (BN_is_zero(t)) { /* special case: p = 2^e + 1 */ 305280304Sjkim if (!BN_nnmod(t, A, p, ctx)) 306280304Sjkim goto end; 307280304Sjkim if (BN_is_zero(t)) { 308280304Sjkim /* special case: a == 0 (mod p) */ 309280304Sjkim BN_zero(ret); 310280304Sjkim err = 0; 311280304Sjkim goto end; 312280304Sjkim } else if (!BN_one(x)) 313280304Sjkim goto end; 314280304Sjkim } else { 315280304Sjkim if (!BN_mod_exp(x, A, t, p, ctx)) 316280304Sjkim goto end; 317280304Sjkim if (BN_is_zero(x)) { 318280304Sjkim /* special case: a == 0 (mod p) */ 319280304Sjkim BN_zero(ret); 320280304Sjkim err = 0; 321280304Sjkim goto end; 322280304Sjkim } 323280304Sjkim } 324280304Sjkim 325280304Sjkim /* b := a*x^2 (= a^q) */ 326280304Sjkim if (!BN_mod_sqr(b, x, p, ctx)) 327280304Sjkim goto end; 328280304Sjkim if (!BN_mod_mul(b, b, A, p, ctx)) 329280304Sjkim goto end; 330280304Sjkim 331280304Sjkim /* x := a*x (= a^((q+1)/2)) */ 332280304Sjkim if (!BN_mod_mul(x, x, A, p, ctx)) 333280304Sjkim goto end; 334280304Sjkim 335280304Sjkim while (1) { 336280304Sjkim /*- 337280304Sjkim * Now b is a^q * y^k for some even k (0 <= k < 2^E 338280304Sjkim * where E refers to the original value of e, which we 339280304Sjkim * don't keep in a variable), and x is a^((q+1)/2) * y^(k/2). 340280304Sjkim * 341280304Sjkim * We have a*b = x^2, 342280304Sjkim * y^2^(e-1) = -1, 343280304Sjkim * b^2^(e-1) = 1. 344280304Sjkim */ 345280304Sjkim 346280304Sjkim if (BN_is_one(b)) { 347280304Sjkim if (!BN_copy(ret, x)) 348280304Sjkim goto end; 349280304Sjkim err = 0; 350280304Sjkim goto vrfy; 351280304Sjkim } 352280304Sjkim 353280304Sjkim /* find smallest i such that b^(2^i) = 1 */ 354280304Sjkim i = 1; 355280304Sjkim if (!BN_mod_sqr(t, b, p, ctx)) 356280304Sjkim goto end; 357280304Sjkim while (!BN_is_one(t)) { 358280304Sjkim i++; 359280304Sjkim if (i == e) { 360280304Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE); 361280304Sjkim goto end; 362280304Sjkim } 363280304Sjkim if (!BN_mod_mul(t, t, t, p, ctx)) 364280304Sjkim goto end; 365280304Sjkim } 366280304Sjkim 367280304Sjkim /* t := y^2^(e - i - 1) */ 368280304Sjkim if (!BN_copy(t, y)) 369280304Sjkim goto end; 370280304Sjkim for (j = e - i - 1; j > 0; j--) { 371280304Sjkim if (!BN_mod_sqr(t, t, p, ctx)) 372280304Sjkim goto end; 373280304Sjkim } 374280304Sjkim if (!BN_mod_mul(y, t, t, p, ctx)) 375280304Sjkim goto end; 376280304Sjkim if (!BN_mod_mul(x, x, t, p, ctx)) 377280304Sjkim goto end; 378280304Sjkim if (!BN_mod_mul(b, b, y, p, ctx)) 379280304Sjkim goto end; 380280304Sjkim e = i; 381280304Sjkim } 382280304Sjkim 383160814Ssimon vrfy: 384280304Sjkim if (!err) { 385280304Sjkim /* 386280304Sjkim * verify the result -- the input might have been not a square (test 387280304Sjkim * added in 0.9.8) 388280304Sjkim */ 389160814Ssimon 390280304Sjkim if (!BN_mod_sqr(x, ret, p, ctx)) 391280304Sjkim err = 1; 392280304Sjkim 393280304Sjkim if (!err && 0 != BN_cmp(x, A)) { 394280304Sjkim BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE); 395280304Sjkim err = 1; 396280304Sjkim } 397280304Sjkim } 398280304Sjkim 399109998Smarkm end: 400280304Sjkim if (err) { 401280304Sjkim if (ret != NULL && ret != in) { 402280304Sjkim BN_clear_free(ret); 403280304Sjkim } 404280304Sjkim ret = NULL; 405280304Sjkim } 406280304Sjkim BN_CTX_end(ctx); 407280304Sjkim bn_check_top(ret); 408280304Sjkim return ret; 409280304Sjkim} 410