bn_sqrt.c revision 109998
1109998Smarkm/* crypto/bn/bn_mod.c */ 2109998Smarkm/* Written by Lenka Fibikova <fibikova@exp-math.uni-essen.de> 3109998Smarkm * and Bodo Moeller for the OpenSSL project. */ 4109998Smarkm/* ==================================================================== 5109998Smarkm * Copyright (c) 1998-2000 The OpenSSL Project. All rights reserved. 6109998Smarkm * 7109998Smarkm * Redistribution and use in source and binary forms, with or without 8109998Smarkm * modification, are permitted provided that the following conditions 9109998Smarkm * are met: 10109998Smarkm * 11109998Smarkm * 1. Redistributions of source code must retain the above copyright 12109998Smarkm * notice, this list of conditions and the following disclaimer. 13109998Smarkm * 14109998Smarkm * 2. Redistributions in binary form must reproduce the above copyright 15109998Smarkm * notice, this list of conditions and the following disclaimer in 16109998Smarkm * the documentation and/or other materials provided with the 17109998Smarkm * distribution. 18109998Smarkm * 19109998Smarkm * 3. All advertising materials mentioning features or use of this 20109998Smarkm * software must display the following acknowledgment: 21109998Smarkm * "This product includes software developed by the OpenSSL Project 22109998Smarkm * for use in the OpenSSL Toolkit. (http://www.openssl.org/)" 23109998Smarkm * 24109998Smarkm * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to 25109998Smarkm * endorse or promote products derived from this software without 26109998Smarkm * prior written permission. For written permission, please contact 27109998Smarkm * openssl-core@openssl.org. 28109998Smarkm * 29109998Smarkm * 5. Products derived from this software may not be called "OpenSSL" 30109998Smarkm * nor may "OpenSSL" appear in their names without prior written 31109998Smarkm * permission of the OpenSSL Project. 32109998Smarkm * 33109998Smarkm * 6. Redistributions of any form whatsoever must retain the following 34109998Smarkm * acknowledgment: 35109998Smarkm * "This product includes software developed by the OpenSSL Project 36109998Smarkm * for use in the OpenSSL Toolkit (http://www.openssl.org/)" 37109998Smarkm * 38109998Smarkm * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY 39109998Smarkm * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 40109998Smarkm * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 41109998Smarkm * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE OpenSSL PROJECT OR 42109998Smarkm * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 43109998Smarkm * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 44109998Smarkm * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 45109998Smarkm * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 46109998Smarkm * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 47109998Smarkm * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 48109998Smarkm * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED 49109998Smarkm * OF THE POSSIBILITY OF SUCH DAMAGE. 50109998Smarkm * ==================================================================== 51109998Smarkm * 52109998Smarkm * This product includes cryptographic software written by Eric Young 53109998Smarkm * (eay@cryptsoft.com). This product includes software written by Tim 54109998Smarkm * Hudson (tjh@cryptsoft.com). 55109998Smarkm * 56109998Smarkm */ 57109998Smarkm 58109998Smarkm#include "cryptlib.h" 59109998Smarkm#include "bn_lcl.h" 60109998Smarkm 61109998Smarkm 62109998SmarkmBIGNUM *BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) 63109998Smarkm/* Returns 'ret' such that 64109998Smarkm * ret^2 == a (mod p), 65109998Smarkm * using the Tonelli/Shanks algorithm (cf. Henri Cohen, "A Course 66109998Smarkm * in Algebraic Computational Number Theory", algorithm 1.5.1). 67109998Smarkm * 'p' must be prime! 68109998Smarkm * If 'a' is not a square, this is not necessarily detected by 69109998Smarkm * the algorithms; a bogus result must be expected in this case. 70109998Smarkm */ 71109998Smarkm { 72109998Smarkm BIGNUM *ret = in; 73109998Smarkm int err = 1; 74109998Smarkm int r; 75109998Smarkm BIGNUM *b, *q, *t, *x, *y; 76109998Smarkm int e, i, j; 77109998Smarkm 78109998Smarkm if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) 79109998Smarkm { 80109998Smarkm if (BN_abs_is_word(p, 2)) 81109998Smarkm { 82109998Smarkm if (ret == NULL) 83109998Smarkm ret = BN_new(); 84109998Smarkm if (ret == NULL) 85109998Smarkm goto end; 86109998Smarkm if (!BN_set_word(ret, BN_is_bit_set(a, 0))) 87109998Smarkm { 88109998Smarkm BN_free(ret); 89109998Smarkm return NULL; 90109998Smarkm } 91109998Smarkm return ret; 92109998Smarkm } 93109998Smarkm 94109998Smarkm BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 95109998Smarkm return(NULL); 96109998Smarkm } 97109998Smarkm 98109998Smarkm if (BN_is_zero(a) || BN_is_one(a)) 99109998Smarkm { 100109998Smarkm if (ret == NULL) 101109998Smarkm ret = BN_new(); 102109998Smarkm if (ret == NULL) 103109998Smarkm goto end; 104109998Smarkm if (!BN_set_word(ret, BN_is_one(a))) 105109998Smarkm { 106109998Smarkm BN_free(ret); 107109998Smarkm return NULL; 108109998Smarkm } 109109998Smarkm return ret; 110109998Smarkm } 111109998Smarkm 112109998Smarkm#if 0 /* if BN_mod_sqrt is used with correct input, this just wastes time */ 113109998Smarkm r = BN_kronecker(a, p, ctx); 114109998Smarkm if (r < -1) return NULL; 115109998Smarkm if (r == -1) 116109998Smarkm { 117109998Smarkm BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE); 118109998Smarkm return(NULL); 119109998Smarkm } 120109998Smarkm#endif 121109998Smarkm 122109998Smarkm BN_CTX_start(ctx); 123109998Smarkm b = BN_CTX_get(ctx); 124109998Smarkm q = BN_CTX_get(ctx); 125109998Smarkm t = BN_CTX_get(ctx); 126109998Smarkm x = BN_CTX_get(ctx); 127109998Smarkm y = BN_CTX_get(ctx); 128109998Smarkm if (y == NULL) goto end; 129109998Smarkm 130109998Smarkm if (ret == NULL) 131109998Smarkm ret = BN_new(); 132109998Smarkm if (ret == NULL) goto end; 133109998Smarkm 134109998Smarkm /* now write |p| - 1 as 2^e*q where q is odd */ 135109998Smarkm e = 1; 136109998Smarkm while (!BN_is_bit_set(p, e)) 137109998Smarkm e++; 138109998Smarkm /* we'll set q later (if needed) */ 139109998Smarkm 140109998Smarkm if (e == 1) 141109998Smarkm { 142109998Smarkm /* The easy case: (|p|-1)/2 is odd, so 2 has an inverse 143109998Smarkm * modulo (|p|-1)/2, and square roots can be computed 144109998Smarkm * directly by modular exponentiation. 145109998Smarkm * We have 146109998Smarkm * 2 * (|p|+1)/4 == 1 (mod (|p|-1)/2), 147109998Smarkm * so we can use exponent (|p|+1)/4, i.e. (|p|-3)/4 + 1. 148109998Smarkm */ 149109998Smarkm if (!BN_rshift(q, p, 2)) goto end; 150109998Smarkm q->neg = 0; 151109998Smarkm if (!BN_add_word(q, 1)) goto end; 152109998Smarkm if (!BN_mod_exp(ret, a, q, p, ctx)) goto end; 153109998Smarkm err = 0; 154109998Smarkm goto end; 155109998Smarkm } 156109998Smarkm 157109998Smarkm if (e == 2) 158109998Smarkm { 159109998Smarkm /* |p| == 5 (mod 8) 160109998Smarkm * 161109998Smarkm * In this case 2 is always a non-square since 162109998Smarkm * Legendre(2,p) = (-1)^((p^2-1)/8) for any odd prime. 163109998Smarkm * So if a really is a square, then 2*a is a non-square. 164109998Smarkm * Thus for 165109998Smarkm * b := (2*a)^((|p|-5)/8), 166109998Smarkm * i := (2*a)*b^2 167109998Smarkm * we have 168109998Smarkm * i^2 = (2*a)^((1 + (|p|-5)/4)*2) 169109998Smarkm * = (2*a)^((p-1)/2) 170109998Smarkm * = -1; 171109998Smarkm * so if we set 172109998Smarkm * x := a*b*(i-1), 173109998Smarkm * then 174109998Smarkm * x^2 = a^2 * b^2 * (i^2 - 2*i + 1) 175109998Smarkm * = a^2 * b^2 * (-2*i) 176109998Smarkm * = a*(-i)*(2*a*b^2) 177109998Smarkm * = a*(-i)*i 178109998Smarkm * = a. 179109998Smarkm * 180109998Smarkm * (This is due to A.O.L. Atkin, 181109998Smarkm * <URL: http://listserv.nodak.edu/scripts/wa.exe?A2=ind9211&L=nmbrthry&O=T&P=562>, 182109998Smarkm * November 1992.) 183109998Smarkm */ 184109998Smarkm 185109998Smarkm /* make sure that a is reduced modulo p */ 186109998Smarkm if (a->neg || BN_ucmp(a, p) >= 0) 187109998Smarkm { 188109998Smarkm if (!BN_nnmod(x, a, p, ctx)) goto end; 189109998Smarkm a = x; /* use x as temporary variable */ 190109998Smarkm } 191109998Smarkm 192109998Smarkm /* t := 2*a */ 193109998Smarkm if (!BN_mod_lshift1_quick(t, a, p)) goto end; 194109998Smarkm 195109998Smarkm /* b := (2*a)^((|p|-5)/8) */ 196109998Smarkm if (!BN_rshift(q, p, 3)) goto end; 197109998Smarkm q->neg = 0; 198109998Smarkm if (!BN_mod_exp(b, t, q, p, ctx)) goto end; 199109998Smarkm 200109998Smarkm /* y := b^2 */ 201109998Smarkm if (!BN_mod_sqr(y, b, p, ctx)) goto end; 202109998Smarkm 203109998Smarkm /* t := (2*a)*b^2 - 1*/ 204109998Smarkm if (!BN_mod_mul(t, t, y, p, ctx)) goto end; 205109998Smarkm if (!BN_sub_word(t, 1)) goto end; 206109998Smarkm 207109998Smarkm /* x = a*b*t */ 208109998Smarkm if (!BN_mod_mul(x, a, b, p, ctx)) goto end; 209109998Smarkm if (!BN_mod_mul(x, x, t, p, ctx)) goto end; 210109998Smarkm 211109998Smarkm if (!BN_copy(ret, x)) goto end; 212109998Smarkm err = 0; 213109998Smarkm goto end; 214109998Smarkm } 215109998Smarkm 216109998Smarkm /* e > 2, so we really have to use the Tonelli/Shanks algorithm. 217109998Smarkm * First, find some y that is not a square. */ 218109998Smarkm if (!BN_copy(q, p)) goto end; /* use 'q' as temp */ 219109998Smarkm q->neg = 0; 220109998Smarkm i = 2; 221109998Smarkm do 222109998Smarkm { 223109998Smarkm /* For efficiency, try small numbers first; 224109998Smarkm * if this fails, try random numbers. 225109998Smarkm */ 226109998Smarkm if (i < 22) 227109998Smarkm { 228109998Smarkm if (!BN_set_word(y, i)) goto end; 229109998Smarkm } 230109998Smarkm else 231109998Smarkm { 232109998Smarkm if (!BN_pseudo_rand(y, BN_num_bits(p), 0, 0)) goto end; 233109998Smarkm if (BN_ucmp(y, p) >= 0) 234109998Smarkm { 235109998Smarkm if (!(p->neg ? BN_add : BN_sub)(y, y, p)) goto end; 236109998Smarkm } 237109998Smarkm /* now 0 <= y < |p| */ 238109998Smarkm if (BN_is_zero(y)) 239109998Smarkm if (!BN_set_word(y, i)) goto end; 240109998Smarkm } 241109998Smarkm 242109998Smarkm r = BN_kronecker(y, q, ctx); /* here 'q' is |p| */ 243109998Smarkm if (r < -1) goto end; 244109998Smarkm if (r == 0) 245109998Smarkm { 246109998Smarkm /* m divides p */ 247109998Smarkm BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 248109998Smarkm goto end; 249109998Smarkm } 250109998Smarkm } 251109998Smarkm while (r == 1 && ++i < 82); 252109998Smarkm 253109998Smarkm if (r != -1) 254109998Smarkm { 255109998Smarkm /* Many rounds and still no non-square -- this is more likely 256109998Smarkm * a bug than just bad luck. 257109998Smarkm * Even if p is not prime, we should have found some y 258109998Smarkm * such that r == -1. 259109998Smarkm */ 260109998Smarkm BNerr(BN_F_BN_MOD_SQRT, BN_R_TOO_MANY_ITERATIONS); 261109998Smarkm goto end; 262109998Smarkm } 263109998Smarkm 264109998Smarkm /* Here's our actual 'q': */ 265109998Smarkm if (!BN_rshift(q, q, e)) goto end; 266109998Smarkm 267109998Smarkm /* Now that we have some non-square, we can find an element 268109998Smarkm * of order 2^e by computing its q'th power. */ 269109998Smarkm if (!BN_mod_exp(y, y, q, p, ctx)) goto end; 270109998Smarkm if (BN_is_one(y)) 271109998Smarkm { 272109998Smarkm BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 273109998Smarkm goto end; 274109998Smarkm } 275109998Smarkm 276109998Smarkm /* Now we know that (if p is indeed prime) there is an integer 277109998Smarkm * k, 0 <= k < 2^e, such that 278109998Smarkm * 279109998Smarkm * a^q * y^k == 1 (mod p). 280109998Smarkm * 281109998Smarkm * As a^q is a square and y is not, k must be even. 282109998Smarkm * q+1 is even, too, so there is an element 283109998Smarkm * 284109998Smarkm * X := a^((q+1)/2) * y^(k/2), 285109998Smarkm * 286109998Smarkm * and it satisfies 287109998Smarkm * 288109998Smarkm * X^2 = a^q * a * y^k 289109998Smarkm * = a, 290109998Smarkm * 291109998Smarkm * so it is the square root that we are looking for. 292109998Smarkm */ 293109998Smarkm 294109998Smarkm /* t := (q-1)/2 (note that q is odd) */ 295109998Smarkm if (!BN_rshift1(t, q)) goto end; 296109998Smarkm 297109998Smarkm /* x := a^((q-1)/2) */ 298109998Smarkm if (BN_is_zero(t)) /* special case: p = 2^e + 1 */ 299109998Smarkm { 300109998Smarkm if (!BN_nnmod(t, a, p, ctx)) goto end; 301109998Smarkm if (BN_is_zero(t)) 302109998Smarkm { 303109998Smarkm /* special case: a == 0 (mod p) */ 304109998Smarkm if (!BN_zero(ret)) goto end; 305109998Smarkm err = 0; 306109998Smarkm goto end; 307109998Smarkm } 308109998Smarkm else 309109998Smarkm if (!BN_one(x)) goto end; 310109998Smarkm } 311109998Smarkm else 312109998Smarkm { 313109998Smarkm if (!BN_mod_exp(x, a, t, p, ctx)) goto end; 314109998Smarkm if (BN_is_zero(x)) 315109998Smarkm { 316109998Smarkm /* special case: a == 0 (mod p) */ 317109998Smarkm if (!BN_zero(ret)) goto end; 318109998Smarkm err = 0; 319109998Smarkm goto end; 320109998Smarkm } 321109998Smarkm } 322109998Smarkm 323109998Smarkm /* b := a*x^2 (= a^q) */ 324109998Smarkm if (!BN_mod_sqr(b, x, p, ctx)) goto end; 325109998Smarkm if (!BN_mod_mul(b, b, a, p, ctx)) goto end; 326109998Smarkm 327109998Smarkm /* x := a*x (= a^((q+1)/2)) */ 328109998Smarkm if (!BN_mod_mul(x, x, a, p, ctx)) goto end; 329109998Smarkm 330109998Smarkm while (1) 331109998Smarkm { 332109998Smarkm /* Now b is a^q * y^k for some even k (0 <= k < 2^E 333109998Smarkm * where E refers to the original value of e, which we 334109998Smarkm * don't keep in a variable), and x is a^((q+1)/2) * y^(k/2). 335109998Smarkm * 336109998Smarkm * We have a*b = x^2, 337109998Smarkm * y^2^(e-1) = -1, 338109998Smarkm * b^2^(e-1) = 1. 339109998Smarkm */ 340109998Smarkm 341109998Smarkm if (BN_is_one(b)) 342109998Smarkm { 343109998Smarkm if (!BN_copy(ret, x)) goto end; 344109998Smarkm err = 0; 345109998Smarkm goto end; 346109998Smarkm } 347109998Smarkm 348109998Smarkm 349109998Smarkm /* find smallest i such that b^(2^i) = 1 */ 350109998Smarkm i = 1; 351109998Smarkm if (!BN_mod_sqr(t, b, p, ctx)) goto end; 352109998Smarkm while (!BN_is_one(t)) 353109998Smarkm { 354109998Smarkm i++; 355109998Smarkm if (i == e) 356109998Smarkm { 357109998Smarkm BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE); 358109998Smarkm goto end; 359109998Smarkm } 360109998Smarkm if (!BN_mod_mul(t, t, t, p, ctx)) goto end; 361109998Smarkm } 362109998Smarkm 363109998Smarkm 364109998Smarkm /* t := y^2^(e - i - 1) */ 365109998Smarkm if (!BN_copy(t, y)) goto end; 366109998Smarkm for (j = e - i - 1; j > 0; j--) 367109998Smarkm { 368109998Smarkm if (!BN_mod_sqr(t, t, p, ctx)) goto end; 369109998Smarkm } 370109998Smarkm if (!BN_mod_mul(y, t, t, p, ctx)) goto end; 371109998Smarkm if (!BN_mod_mul(x, x, t, p, ctx)) goto end; 372109998Smarkm if (!BN_mod_mul(b, b, y, p, ctx)) goto end; 373109998Smarkm e = i; 374109998Smarkm } 375109998Smarkm 376109998Smarkm end: 377109998Smarkm if (err) 378109998Smarkm { 379109998Smarkm if (ret != NULL && ret != in) 380109998Smarkm { 381109998Smarkm BN_clear_free(ret); 382109998Smarkm } 383109998Smarkm ret = NULL; 384109998Smarkm } 385109998Smarkm BN_CTX_end(ctx); 386109998Smarkm return ret; 387109998Smarkm } 388