1 /* 2 * Copyright 2000-2018 The OpenSSL Project Authors. All Rights Reserved. 3 * 4 * Licensed under the OpenSSL license (the "License"). You may not use 5 * this file except in compliance with the License. You can obtain a copy 6 * in the file LICENSE in the source distribution or at 7 * https://www.openssl.org/source/license.html 8 */ 9 10 #include "internal/cryptlib.h" 11 #include "bn_lcl.h" 12 13 BIGNUM *BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) 14 /* 15 * Returns 'ret' such that ret^2 == a (mod p), using the Tonelli/Shanks 16 * algorithm (cf. Henri Cohen, "A Course in Algebraic Computational Number 17 * Theory", algorithm 1.5.1). 'p' must be prime! 18 */ 19 { 20 BIGNUM *ret = in; 21 int err = 1; 22 int r; 23 BIGNUM *A, *b, *q, *t, *x, *y; 24 int e, i, j; 25 26 if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) { 27 if (BN_abs_is_word(p, 2)) { 28 if (ret == NULL) 29 ret = BN_new(); 30 if (ret == NULL) 31 goto end; 32 if (!BN_set_word(ret, BN_is_bit_set(a, 0))) { 33 if (ret != in) 34 BN_free(ret); 35 return NULL; 36 } 37 bn_check_top(ret); 38 return ret; 39 } 40 41 BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 42 return NULL; 43 } 44 45 if (BN_is_zero(a) || BN_is_one(a)) { 46 if (ret == NULL) 47 ret = BN_new(); 48 if (ret == NULL) 49 goto end; 50 if (!BN_set_word(ret, BN_is_one(a))) { 51 if (ret != in) 52 BN_free(ret); 53 return NULL; 54 } 55 bn_check_top(ret); 56 return ret; 57 } 58 59 BN_CTX_start(ctx); 60 A = BN_CTX_get(ctx); 61 b = BN_CTX_get(ctx); 62 q = BN_CTX_get(ctx); 63 t = BN_CTX_get(ctx); 64 x = BN_CTX_get(ctx); 65 y = BN_CTX_get(ctx); 66 if (y == NULL) 67 goto end; 68 69 if (ret == NULL) 70 ret = BN_new(); 71 if (ret == NULL) 72 goto end; 73 74 /* A = a mod p */ 75 if (!BN_nnmod(A, a, p, ctx)) 76 goto end; 77 78 /* now write |p| - 1 as 2^e*q where q is odd */ 79 e = 1; 80 while (!BN_is_bit_set(p, e)) 81 e++; 82 /* we'll set q later (if needed) */ 83 84 if (e == 1) { 85 /*- 86 * The easy case: (|p|-1)/2 is odd, so 2 has an inverse 87 * modulo (|p|-1)/2, and square roots can be computed 88 * directly by modular exponentiation. 89 * We have 90 * 2 * (|p|+1)/4 == 1 (mod (|p|-1)/2), 91 * so we can use exponent (|p|+1)/4, i.e. (|p|-3)/4 + 1. 92 */ 93 if (!BN_rshift(q, p, 2)) 94 goto end; 95 q->neg = 0; 96 if (!BN_add_word(q, 1)) 97 goto end; 98 if (!BN_mod_exp(ret, A, q, p, ctx)) 99 goto end; 100 err = 0; 101 goto vrfy; 102 } 103 104 if (e == 2) { 105 /*- 106 * |p| == 5 (mod 8) 107 * 108 * In this case 2 is always a non-square since 109 * Legendre(2,p) = (-1)^((p^2-1)/8) for any odd prime. 110 * So if a really is a square, then 2*a is a non-square. 111 * Thus for 112 * b := (2*a)^((|p|-5)/8), 113 * i := (2*a)*b^2 114 * we have 115 * i^2 = (2*a)^((1 + (|p|-5)/4)*2) 116 * = (2*a)^((p-1)/2) 117 * = -1; 118 * so if we set 119 * x := a*b*(i-1), 120 * then 121 * x^2 = a^2 * b^2 * (i^2 - 2*i + 1) 122 * = a^2 * b^2 * (-2*i) 123 * = a*(-i)*(2*a*b^2) 124 * = a*(-i)*i 125 * = a. 126 * 127 * (This is due to A.O.L. Atkin, 128 * <URL: http://listserv.nodak.edu/scripts/wa.exe?A2=ind9211&L=nmbrthry&O=T&P=562>, 129 * November 1992.) 130 */ 131 132 /* t := 2*a */ 133 if (!BN_mod_lshift1_quick(t, A, p)) 134 goto end; 135 136 /* b := (2*a)^((|p|-5)/8) */ 137 if (!BN_rshift(q, p, 3)) 138 goto end; 139 q->neg = 0; 140 if (!BN_mod_exp(b, t, q, p, ctx)) 141 goto end; 142 143 /* y := b^2 */ 144 if (!BN_mod_sqr(y, b, p, ctx)) 145 goto end; 146 147 /* t := (2*a)*b^2 - 1 */ 148 if (!BN_mod_mul(t, t, y, p, ctx)) 149 goto end; 150 if (!BN_sub_word(t, 1)) 151 goto end; 152 153 /* x = a*b*t */ 154 if (!BN_mod_mul(x, A, b, p, ctx)) 155 goto end; 156 if (!BN_mod_mul(x, x, t, p, ctx)) 157 goto end; 158 159 if (!BN_copy(ret, x)) 160 goto end; 161 err = 0; 162 goto vrfy; 163 } 164 165 /* 166 * e > 2, so we really have to use the Tonelli/Shanks algorithm. First, 167 * find some y that is not a square. 168 */ 169 if (!BN_copy(q, p)) 170 goto end; /* use 'q' as temp */ 171 q->neg = 0; 172 i = 2; 173 do { 174 /* 175 * For efficiency, try small numbers first; if this fails, try random 176 * numbers. 177 */ 178 if (i < 22) { 179 if (!BN_set_word(y, i)) 180 goto end; 181 } else { 182 if (!BN_priv_rand(y, BN_num_bits(p), 0, 0)) 183 goto end; 184 if (BN_ucmp(y, p) >= 0) { 185 if (!(p->neg ? BN_add : BN_sub) (y, y, p)) 186 goto end; 187 } 188 /* now 0 <= y < |p| */ 189 if (BN_is_zero(y)) 190 if (!BN_set_word(y, i)) 191 goto end; 192 } 193 194 r = BN_kronecker(y, q, ctx); /* here 'q' is |p| */ 195 if (r < -1) 196 goto end; 197 if (r == 0) { 198 /* m divides p */ 199 BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 200 goto end; 201 } 202 } 203 while (r == 1 && ++i < 82); 204 205 if (r != -1) { 206 /* 207 * Many rounds and still no non-square -- this is more likely a bug 208 * than just bad luck. Even if p is not prime, we should have found 209 * some y such that r == -1. 210 */ 211 BNerr(BN_F_BN_MOD_SQRT, BN_R_TOO_MANY_ITERATIONS); 212 goto end; 213 } 214 215 /* Here's our actual 'q': */ 216 if (!BN_rshift(q, q, e)) 217 goto end; 218 219 /* 220 * Now that we have some non-square, we can find an element of order 2^e 221 * by computing its q'th power. 222 */ 223 if (!BN_mod_exp(y, y, q, p, ctx)) 224 goto end; 225 if (BN_is_one(y)) { 226 BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME); 227 goto end; 228 } 229 230 /*- 231 * Now we know that (if p is indeed prime) there is an integer 232 * k, 0 <= k < 2^e, such that 233 * 234 * a^q * y^k == 1 (mod p). 235 * 236 * As a^q is a square and y is not, k must be even. 237 * q+1 is even, too, so there is an element 238 * 239 * X := a^((q+1)/2) * y^(k/2), 240 * 241 * and it satisfies 242 * 243 * X^2 = a^q * a * y^k 244 * = a, 245 * 246 * so it is the square root that we are looking for. 247 */ 248 249 /* t := (q-1)/2 (note that q is odd) */ 250 if (!BN_rshift1(t, q)) 251 goto end; 252 253 /* x := a^((q-1)/2) */ 254 if (BN_is_zero(t)) { /* special case: p = 2^e + 1 */ 255 if (!BN_nnmod(t, A, p, ctx)) 256 goto end; 257 if (BN_is_zero(t)) { 258 /* special case: a == 0 (mod p) */ 259 BN_zero(ret); 260 err = 0; 261 goto end; 262 } else if (!BN_one(x)) 263 goto end; 264 } else { 265 if (!BN_mod_exp(x, A, t, p, ctx)) 266 goto end; 267 if (BN_is_zero(x)) { 268 /* special case: a == 0 (mod p) */ 269 BN_zero(ret); 270 err = 0; 271 goto end; 272 } 273 } 274 275 /* b := a*x^2 (= a^q) */ 276 if (!BN_mod_sqr(b, x, p, ctx)) 277 goto end; 278 if (!BN_mod_mul(b, b, A, p, ctx)) 279 goto end; 280 281 /* x := a*x (= a^((q+1)/2)) */ 282 if (!BN_mod_mul(x, x, A, p, ctx)) 283 goto end; 284 285 while (1) { 286 /*- 287 * Now b is a^q * y^k for some even k (0 <= k < 2^E 288 * where E refers to the original value of e, which we 289 * don't keep in a variable), and x is a^((q+1)/2) * y^(k/2). 290 * 291 * We have a*b = x^2, 292 * y^2^(e-1) = -1, 293 * b^2^(e-1) = 1. 294 */ 295 296 if (BN_is_one(b)) { 297 if (!BN_copy(ret, x)) 298 goto end; 299 err = 0; 300 goto vrfy; 301 } 302 303 /* find smallest i such that b^(2^i) = 1 */ 304 i = 1; 305 if (!BN_mod_sqr(t, b, p, ctx)) 306 goto end; 307 while (!BN_is_one(t)) { 308 i++; 309 if (i == e) { 310 BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE); 311 goto end; 312 } 313 if (!BN_mod_mul(t, t, t, p, ctx)) 314 goto end; 315 } 316 317 /* t := y^2^(e - i - 1) */ 318 if (!BN_copy(t, y)) 319 goto end; 320 for (j = e - i - 1; j > 0; j--) { 321 if (!BN_mod_sqr(t, t, p, ctx)) 322 goto end; 323 } 324 if (!BN_mod_mul(y, t, t, p, ctx)) 325 goto end; 326 if (!BN_mod_mul(x, x, t, p, ctx)) 327 goto end; 328 if (!BN_mod_mul(b, b, y, p, ctx)) 329 goto end; 330 e = i; 331 } 332 333 vrfy: 334 if (!err) { 335 /* 336 * verify the result -- the input might have been not a square (test 337 * added in 0.9.8) 338 */ 339 340 if (!BN_mod_sqr(x, ret, p, ctx)) 341 err = 1; 342 343 if (!err && 0 != BN_cmp(x, A)) { 344 BNerr(BN_F_BN_MOD_SQRT, BN_R_NOT_A_SQUARE); 345 err = 1; 346 } 347 } 348 349 end: 350 if (err) { 351 if (ret != in) 352 BN_clear_free(ret); 353 ret = NULL; 354 } 355 BN_CTX_end(ctx); 356 bn_check_top(ret); 357 return ret; 358 } 359