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