1 /* $OpenBSD: bn_mod_sqrt.c,v 1.3 2023/08/03 18:53:55 tb Exp $ */
2
3 /*
4 * Copyright (c) 2022 Theo Buehler <tb@openbsd.org>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19 #include <openssl/err.h>
20
21 #include "bn_local.h"
22
23 /*
24 * Tonelli-Shanks according to H. Cohen "A Course in Computational Algebraic
25 * Number Theory", Section 1.5.1, Springer GTM volume 138, Berlin, 1996.
26 *
27 * Under the assumption that p is prime and a is a quadratic residue, we know:
28 *
29 * a^[(p-1)/2] = 1 (mod p). (*)
30 *
31 * To find a square root of a (mod p), we handle three cases of increasing
32 * complexity. In the first two cases, we can compute a square root using an
33 * explicit formula, thus avoiding the probabilistic nature of Tonelli-Shanks.
34 *
35 * 1. p = 3 (mod 4).
36 *
37 * Set n = (p+1)/4. Then 2n = 1 + (p-1)/2 and (*) shows that x = a^n (mod p)
38 * is a square root of a: x^2 = a^(2n) = a * a^[(p-1)/2] = a (mod p).
39 *
40 * 2. p = 5 (mod 8).
41 *
42 * This uses a simplification due to Atkin. By Theorem 1.4.7 and 1.4.9, the
43 * Kronecker symbol (2/p) evaluates to (-1)^[(p^2-1)/8]. From p = 5 (mod 8)
44 * we get (p^2-1)/8 = 1 (mod 2), so (2/p) = -1, and thus
45 *
46 * 2^[(p-1)/2] = -1 (mod p). (**)
47 *
48 * Set b = (2a)^[(p-5)/8]. With (p-1)/2 = 2 + (p-5)/2, (*) and (**) show
49 *
50 * i = 2 a b^2 is a square root of -1 (mod p).
51 *
52 * Indeed, i^2 = 2^2 a^2 b^4 = 2^[(p-1)/2] a^[(p-1)/2] = -1 (mod p). Because
53 * of (i-1)^2 = -2i (mod p) and i (-i) = 1 (mod p), a square root of a is
54 *
55 * x = a b (i-1)
56 *
57 * as x^2 = a^2 b^2 (-2i) = a (2 a b^2) (-i) = a (mod p).
58 *
59 * 3. p = 1 (mod 8).
60 *
61 * This is the Tonelli-Shanks algorithm. For a prime p, the multiplicative
62 * group of GF(p) is cyclic of order p - 1 = 2^s q, with odd q. Denote its
63 * 2-Sylow subgroup by S. It is cyclic of order 2^s. The squares in S have
64 * order dividing 2^(s-1). They are the even powers of any generator z of S.
65 * If a is a quadratic residue, 1 = a^[(p-1)/2] = (a^q)^[2^(s-1)], so b = a^q
66 * is a square in S. Therefore there is an integer k such that b z^(2k) = 1.
67 * Set x = a^[(q+1)/2] z^k, and find x^2 = a (mod p).
68 *
69 * The problem is thus reduced to finding a generator z of the 2-Sylow
70 * subgroup S of GF(p)* and finding k. An iterative constructions avoids
71 * the need for an explicit k, a generator is found by a randomized search.
72 *
73 * While we do not actually know that p is a prime number, we can still apply
74 * the formulas in cases 1 and 2 and verify that we have indeed found a square
75 * root of p. Similarly, in case 3, we can try to find a quadratic non-residue,
76 * which will fail for example if p is a square. The iterative construction
77 * may or may not find a candidate square root which we can then validate.
78 */
79
80 /*
81 * Handle the cases where p is 2, p isn't odd or p is one. Since BN_mod_sqrt()
82 * can run on untrusted data, a primality check is too expensive. Also treat
83 * the obvious cases where a is 0 or 1.
84 */
85
86 static int
bn_mod_sqrt_trivial_cases(int * done,BIGNUM * out_sqrt,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)87 bn_mod_sqrt_trivial_cases(int *done, BIGNUM *out_sqrt, const BIGNUM *a,
88 const BIGNUM *p, BN_CTX *ctx)
89 {
90 *done = 1;
91
92 if (BN_abs_is_word(p, 2))
93 return BN_set_word(out_sqrt, BN_is_odd(a));
94
95 if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) {
96 BNerror(BN_R_P_IS_NOT_PRIME);
97 return 0;
98 }
99
100 if (BN_is_zero(a) || BN_is_one(a))
101 return BN_set_word(out_sqrt, BN_is_one(a));
102
103 *done = 0;
104
105 return 1;
106 }
107
108 /*
109 * Case 1. We know that (a/p) = 1 and that p = 3 (mod 4).
110 */
111
112 static int
bn_mod_sqrt_p_is_3_mod_4(BIGNUM * out_sqrt,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)113 bn_mod_sqrt_p_is_3_mod_4(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
114 BN_CTX *ctx)
115 {
116 BIGNUM *n;
117 int ret = 0;
118
119 BN_CTX_start(ctx);
120
121 if ((n = BN_CTX_get(ctx)) == NULL)
122 goto err;
123
124 /* Calculate n = (|p| + 1) / 4. */
125 if (!BN_uadd(n, p, BN_value_one()))
126 goto err;
127 if (!BN_rshift(n, n, 2))
128 goto err;
129
130 /* By case 1 above, out_sqrt = a^n is a square root of a (mod p). */
131 if (!BN_mod_exp_ct(out_sqrt, a, n, p, ctx))
132 goto err;
133
134 ret = 1;
135
136 err:
137 BN_CTX_end(ctx);
138
139 return ret;
140 }
141
142 /*
143 * Case 2. We know that (a/p) = 1 and that p = 5 (mod 8).
144 */
145
146 static int
bn_mod_sqrt_p_is_5_mod_8(BIGNUM * out_sqrt,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)147 bn_mod_sqrt_p_is_5_mod_8(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
148 BN_CTX *ctx)
149 {
150 BIGNUM *b, *i, *n, *tmp;
151 int ret = 0;
152
153 BN_CTX_start(ctx);
154
155 if ((b = BN_CTX_get(ctx)) == NULL)
156 goto err;
157 if ((i = BN_CTX_get(ctx)) == NULL)
158 goto err;
159 if ((n = BN_CTX_get(ctx)) == NULL)
160 goto err;
161 if ((tmp = BN_CTX_get(ctx)) == NULL)
162 goto err;
163
164 /* Calculate n = (|p| - 5) / 8. Since p = 5 (mod 8), simply shift. */
165 if (!BN_rshift(n, p, 3))
166 goto err;
167 BN_set_negative(n, 0);
168
169 /* Compute tmp = 2a (mod p) for later use. */
170 if (!BN_mod_lshift1(tmp, a, p, ctx))
171 goto err;
172
173 /* Calculate b = (2a)^n (mod p). */
174 if (!BN_mod_exp_ct(b, tmp, n, p, ctx))
175 goto err;
176
177 /* Calculate i = 2 a b^2 (mod p). */
178 if (!BN_mod_sqr(i, b, p, ctx))
179 goto err;
180 if (!BN_mod_mul(i, tmp, i, p, ctx))
181 goto err;
182
183 /* A square root is out_sqrt = a b (i-1) (mod p). */
184 if (!BN_sub_word(i, 1))
185 goto err;
186 if (!BN_mod_mul(out_sqrt, a, b, p, ctx))
187 goto err;
188 if (!BN_mod_mul(out_sqrt, out_sqrt, i, p, ctx))
189 goto err;
190
191 ret = 1;
192
193 err:
194 BN_CTX_end(ctx);
195
196 return ret;
197 }
198
199 /*
200 * Case 3. We know that (a/p) = 1 and that p = 1 (mod 8).
201 */
202
203 /*
204 * Simple helper. To find a generator of the 2-Sylow subgroup of GF(p)*, we
205 * need to find a quadratic non-residue of p, i.e., n such that (n/p) = -1.
206 */
207
208 static int
bn_mod_sqrt_n_is_non_residue(int * is_non_residue,const BIGNUM * n,const BIGNUM * p,BN_CTX * ctx)209 bn_mod_sqrt_n_is_non_residue(int *is_non_residue, const BIGNUM *n,
210 const BIGNUM *p, BN_CTX *ctx)
211 {
212 switch (BN_kronecker(n, p, ctx)) {
213 case -1:
214 *is_non_residue = 1;
215 return 1;
216 case 1:
217 *is_non_residue = 0;
218 return 1;
219 case 0:
220 /* n divides p, so ... */
221 BNerror(BN_R_P_IS_NOT_PRIME);
222 return 0;
223 default:
224 return 0;
225 }
226 }
227
228 /*
229 * The following is the only non-deterministic part preparing Tonelli-Shanks.
230 *
231 * If we find n such that (n/p) = -1, then n^q (mod p) is a generator of the
232 * 2-Sylow subgroup of GF(p)*. To find such n, first try some small numbers,
233 * then random ones.
234 */
235
236 static int
bn_mod_sqrt_find_sylow_generator(BIGNUM * out_generator,const BIGNUM * p,const BIGNUM * q,BN_CTX * ctx)237 bn_mod_sqrt_find_sylow_generator(BIGNUM *out_generator, const BIGNUM *p,
238 const BIGNUM *q, BN_CTX *ctx)
239 {
240 BIGNUM *n, *p_abs;
241 int i, is_non_residue;
242 int ret = 0;
243
244 BN_CTX_start(ctx);
245
246 if ((n = BN_CTX_get(ctx)) == NULL)
247 goto err;
248 if ((p_abs = BN_CTX_get(ctx)) == NULL)
249 goto err;
250
251 for (i = 2; i < 32; i++) {
252 if (!BN_set_word(n, i))
253 goto err;
254 if (!bn_mod_sqrt_n_is_non_residue(&is_non_residue, n, p, ctx))
255 goto err;
256 if (is_non_residue)
257 goto found;
258 }
259
260 if (!bn_copy(p_abs, p))
261 goto err;
262 BN_set_negative(p_abs, 0);
263
264 for (i = 0; i < 128; i++) {
265 if (!bn_rand_interval(n, 32, p_abs))
266 goto err;
267 if (!bn_mod_sqrt_n_is_non_residue(&is_non_residue, n, p, ctx))
268 goto err;
269 if (is_non_residue)
270 goto found;
271 }
272
273 /*
274 * The probability to get here is < 2^(-128) for prime p. For squares
275 * it is easy: for p = 1369 = 37^2 this happens in ~3% of runs.
276 */
277
278 BNerror(BN_R_TOO_MANY_ITERATIONS);
279 goto err;
280
281 found:
282 /*
283 * If p is prime, n^q generates the 2-Sylow subgroup S of GF(p)*.
284 */
285
286 if (!BN_mod_exp_ct(out_generator, n, q, p, ctx))
287 goto err;
288
289 /* Sanity: p is not necessarily prime, so we could have found 0 or 1. */
290 if (BN_is_zero(out_generator) || BN_is_one(out_generator)) {
291 BNerror(BN_R_P_IS_NOT_PRIME);
292 goto err;
293 }
294
295 ret = 1;
296
297 err:
298 BN_CTX_end(ctx);
299
300 return ret;
301 }
302
303 /*
304 * Initialization step for Tonelli-Shanks.
305 *
306 * In the end, b = a^q (mod p) and x = a^[(q+1)/2] (mod p). Cohen optimizes this
307 * to minimize taking powers of a. This is a bit confusing and distracting, so
308 * factor this into a separate function.
309 */
310
311 static int
bn_mod_sqrt_tonelli_shanks_initialize(BIGNUM * b,BIGNUM * x,const BIGNUM * a,const BIGNUM * p,const BIGNUM * q,BN_CTX * ctx)312 bn_mod_sqrt_tonelli_shanks_initialize(BIGNUM *b, BIGNUM *x, const BIGNUM *a,
313 const BIGNUM *p, const BIGNUM *q, BN_CTX *ctx)
314 {
315 BIGNUM *k;
316 int ret = 0;
317
318 BN_CTX_start(ctx);
319
320 if ((k = BN_CTX_get(ctx)) == NULL)
321 goto err;
322
323 /* k = (q-1)/2. Since q is odd, we can shift. */
324 if (!BN_rshift1(k, q))
325 goto err;
326
327 /* x = a^[(q-1)/2] (mod p). */
328 if (!BN_mod_exp_ct(x, a, k, p, ctx))
329 goto err;
330
331 /* b = ax^2 = a^q (mod p). */
332 if (!BN_mod_sqr(b, x, p, ctx))
333 goto err;
334 if (!BN_mod_mul(b, a, b, p, ctx))
335 goto err;
336
337 /* x = ax = a^[(q+1)/2] (mod p). */
338 if (!BN_mod_mul(x, a, x, p, ctx))
339 goto err;
340
341 ret = 1;
342
343 err:
344 BN_CTX_end(ctx);
345
346 return ret;
347 }
348
349 /*
350 * Find smallest exponent m such that b^(2^m) = 1 (mod p). Assuming that a
351 * is a quadratic residue and p is a prime, we know that 1 <= m < r.
352 */
353
354 static int
bn_mod_sqrt_tonelli_shanks_find_exponent(int * out_exponent,const BIGNUM * b,const BIGNUM * p,int r,BN_CTX * ctx)355 bn_mod_sqrt_tonelli_shanks_find_exponent(int *out_exponent, const BIGNUM *b,
356 const BIGNUM *p, int r, BN_CTX *ctx)
357 {
358 BIGNUM *x;
359 int m;
360 int ret = 0;
361
362 BN_CTX_start(ctx);
363
364 if ((x = BN_CTX_get(ctx)) == NULL)
365 goto err;
366
367 /*
368 * If r <= 1, the Tonelli-Shanks iteration should have terminated as
369 * r == 1 implies b == 1.
370 */
371 if (r <= 1) {
372 BNerror(BN_R_P_IS_NOT_PRIME);
373 goto err;
374 }
375
376 /*
377 * Sanity check to ensure taking squares actually does something:
378 * If b is 1, the Tonelli-Shanks iteration should have terminated.
379 * If b is 0, something's very wrong, in particular p can't be prime.
380 */
381 if (BN_is_zero(b) || BN_is_one(b)) {
382 BNerror(BN_R_P_IS_NOT_PRIME);
383 goto err;
384 }
385
386 if (!bn_copy(x, b))
387 goto err;
388
389 for (m = 1; m < r; m++) {
390 if (!BN_mod_sqr(x, x, p, ctx))
391 goto err;
392 if (BN_is_one(x))
393 break;
394 }
395
396 if (m >= r) {
397 /* This means a is not a quadratic residue. As (a/p) = 1, ... */
398 BNerror(BN_R_P_IS_NOT_PRIME);
399 goto err;
400 }
401
402 *out_exponent = m;
403
404 ret = 1;
405
406 err:
407 BN_CTX_end(ctx);
408
409 return ret;
410 }
411
412 /*
413 * The update step. With the minimal m such that b^(2^m) = 1 (mod m),
414 * set t = y^[2^(r-m-1)] (mod p) and update x = xt, y = t^2, b = by.
415 * This preserves the loop invariants a b = x^2, y^[2^(r-1)] = -1 and
416 * b^[2^(r-1)] = 1.
417 */
418
419 static int
bn_mod_sqrt_tonelli_shanks_update(BIGNUM * b,BIGNUM * x,BIGNUM * y,const BIGNUM * p,int m,int r,BN_CTX * ctx)420 bn_mod_sqrt_tonelli_shanks_update(BIGNUM *b, BIGNUM *x, BIGNUM *y,
421 const BIGNUM *p, int m, int r, BN_CTX *ctx)
422 {
423 BIGNUM *t;
424 int ret = 0;
425
426 BN_CTX_start(ctx);
427
428 if ((t = BN_CTX_get(ctx)) == NULL)
429 goto err;
430
431 /* t = y^[2^(r-m-1)] (mod p). */
432 if (!BN_set_bit(t, r - m - 1))
433 goto err;
434 if (!BN_mod_exp_ct(t, y, t, p, ctx))
435 goto err;
436
437 /* x = xt (mod p). */
438 if (!BN_mod_mul(x, x, t, p, ctx))
439 goto err;
440
441 /* y = t^2 = y^[2^(r-m)] (mod p). */
442 if (!BN_mod_sqr(y, t, p, ctx))
443 goto err;
444
445 /* b = by (mod p). */
446 if (!BN_mod_mul(b, b, y, p, ctx))
447 goto err;
448
449 ret = 1;
450
451 err:
452 BN_CTX_end(ctx);
453
454 return ret;
455 }
456
457 static int
bn_mod_sqrt_p_is_1_mod_8(BIGNUM * out_sqrt,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)458 bn_mod_sqrt_p_is_1_mod_8(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
459 BN_CTX *ctx)
460 {
461 BIGNUM *b, *q, *x, *y;
462 int e, m, r;
463 int ret = 0;
464
465 BN_CTX_start(ctx);
466
467 if ((b = BN_CTX_get(ctx)) == NULL)
468 goto err;
469 if ((q = BN_CTX_get(ctx)) == NULL)
470 goto err;
471 if ((x = BN_CTX_get(ctx)) == NULL)
472 goto err;
473 if ((y = BN_CTX_get(ctx)) == NULL)
474 goto err;
475
476 /*
477 * Factor p - 1 = 2^e q with odd q. Since p = 1 (mod 8), we know e >= 3.
478 */
479
480 e = 1;
481 while (!BN_is_bit_set(p, e))
482 e++;
483 if (!BN_rshift(q, p, e))
484 goto err;
485
486 if (!bn_mod_sqrt_find_sylow_generator(y, p, q, ctx))
487 goto err;
488
489 /*
490 * Set b = a^q (mod p) and x = a^[(q+1)/2] (mod p).
491 */
492 if (!bn_mod_sqrt_tonelli_shanks_initialize(b, x, a, p, q, ctx))
493 goto err;
494
495 /*
496 * The Tonelli-Shanks iteration. Starting with r = e, the following loop
497 * invariants hold at the start of the loop.
498 *
499 * a b = x^2 (mod p)
500 * y^[2^(r-1)] = -1 (mod p)
501 * b^[2^(r-1)] = 1 (mod p)
502 *
503 * In particular, if b = 1 (mod p), x is a square root of a.
504 *
505 * Since p - 1 = 2^e q, we have 2^(e-1) q = (p - 1) / 2, so in the first
506 * iteration this follows from (a/p) = 1, (n/p) = -1, y = n^q, b = a^q.
507 *
508 * In subsequent iterations, t = y^[2^(r-m-1)], where m is the smallest
509 * m such that b^(2^m) = 1. With x = xt (mod p) and b = bt^2 (mod p) the
510 * first invariant is preserved, the second and third follow from
511 * y = t^2 (mod p) and r = m as well as the choice of m.
512 *
513 * Finally, r is strictly decreasing in each iteration. If p is prime,
514 * let S be the 2-Sylow subgroup of GF(p)*. We can prove the algorithm
515 * stops: Let S_r be the subgroup of S consisting of elements of order
516 * dividing 2^r. Then S_r = <y> and b is in S_(r-1). The S_r form a
517 * descending filtration of S and when r = 1, then b = 1.
518 */
519
520 for (r = e; r >= 1; r = m) {
521 /*
522 * Termination condition. If b == 1 then x is a square root.
523 */
524 if (BN_is_one(b))
525 goto done;
526
527 /* Find smallest exponent 1 <= m < r such that b^(2^m) == 1. */
528 if (!bn_mod_sqrt_tonelli_shanks_find_exponent(&m, b, p, r, ctx))
529 goto err;
530
531 /*
532 * With t = y^[2^(r-m-1)], update x = xt, y = t^2, b = by.
533 */
534 if (!bn_mod_sqrt_tonelli_shanks_update(b, x, y, p, m, r, ctx))
535 goto err;
536
537 /*
538 * Sanity check to make sure we don't loop indefinitely.
539 * bn_mod_sqrt_tonelli_shanks_find_exponent() ensures m < r.
540 */
541 if (r <= m)
542 goto err;
543 }
544
545 /*
546 * If p is prime, we should not get here.
547 */
548
549 BNerror(BN_R_NOT_A_SQUARE);
550 goto err;
551
552 done:
553 if (!bn_copy(out_sqrt, x))
554 goto err;
555
556 ret = 1;
557
558 err:
559 BN_CTX_end(ctx);
560
561 return ret;
562 }
563
564 /*
565 * Choose the smaller of sqrt and |p| - sqrt.
566 */
567
568 static int
bn_mod_sqrt_normalize(BIGNUM * sqrt,const BIGNUM * p,BN_CTX * ctx)569 bn_mod_sqrt_normalize(BIGNUM *sqrt, const BIGNUM *p, BN_CTX *ctx)
570 {
571 BIGNUM *x;
572 int ret = 0;
573
574 BN_CTX_start(ctx);
575
576 if ((x = BN_CTX_get(ctx)) == NULL)
577 goto err;
578
579 if (!BN_lshift1(x, sqrt))
580 goto err;
581
582 if (BN_ucmp(x, p) > 0) {
583 if (!BN_usub(sqrt, p, sqrt))
584 goto err;
585 }
586
587 ret = 1;
588
589 err:
590 BN_CTX_end(ctx);
591
592 return ret;
593 }
594
595 /*
596 * Verify that a = (sqrt_a)^2 (mod p). Requires that a is reduced (mod p).
597 */
598
599 static int
bn_mod_sqrt_verify(const BIGNUM * a,const BIGNUM * sqrt_a,const BIGNUM * p,BN_CTX * ctx)600 bn_mod_sqrt_verify(const BIGNUM *a, const BIGNUM *sqrt_a, const BIGNUM *p,
601 BN_CTX *ctx)
602 {
603 BIGNUM *x;
604 int ret = 0;
605
606 BN_CTX_start(ctx);
607
608 if ((x = BN_CTX_get(ctx)) == NULL)
609 goto err;
610
611 if (!BN_mod_sqr(x, sqrt_a, p, ctx))
612 goto err;
613
614 if (BN_cmp(x, a) != 0) {
615 BNerror(BN_R_NOT_A_SQUARE);
616 goto err;
617 }
618
619 ret = 1;
620
621 err:
622 BN_CTX_end(ctx);
623
624 return ret;
625 }
626
627 static int
bn_mod_sqrt_internal(BIGNUM * out_sqrt,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)628 bn_mod_sqrt_internal(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
629 BN_CTX *ctx)
630 {
631 BIGNUM *a_mod_p, *sqrt;
632 BN_ULONG lsw;
633 int done;
634 int kronecker;
635 int ret = 0;
636
637 BN_CTX_start(ctx);
638
639 if ((a_mod_p = BN_CTX_get(ctx)) == NULL)
640 goto err;
641 if ((sqrt = BN_CTX_get(ctx)) == NULL)
642 goto err;
643
644 if (!BN_nnmod(a_mod_p, a, p, ctx))
645 goto err;
646
647 if (!bn_mod_sqrt_trivial_cases(&done, sqrt, a_mod_p, p, ctx))
648 goto err;
649 if (done)
650 goto verify;
651
652 /*
653 * Make sure that the Kronecker symbol (a/p) == 1. In case p is prime
654 * this is equivalent to a having a square root (mod p). The cost of
655 * BN_kronecker() is O(log^2(n)). This is small compared to the cost
656 * O(log^4(n)) of Tonelli-Shanks.
657 */
658
659 if ((kronecker = BN_kronecker(a_mod_p, p, ctx)) == -2)
660 goto err;
661 if (kronecker <= 0) {
662 /* This error is only accurate if p is known to be a prime. */
663 BNerror(BN_R_NOT_A_SQUARE);
664 goto err;
665 }
666
667 lsw = BN_lsw(p);
668
669 if (lsw % 4 == 3) {
670 if (!bn_mod_sqrt_p_is_3_mod_4(sqrt, a_mod_p, p, ctx))
671 goto err;
672 } else if (lsw % 8 == 5) {
673 if (!bn_mod_sqrt_p_is_5_mod_8(sqrt, a_mod_p, p, ctx))
674 goto err;
675 } else if (lsw % 8 == 1) {
676 if (!bn_mod_sqrt_p_is_1_mod_8(sqrt, a_mod_p, p, ctx))
677 goto err;
678 } else {
679 /* Impossible to hit since the trivial cases ensure p is odd. */
680 BNerror(BN_R_P_IS_NOT_PRIME);
681 goto err;
682 }
683
684 if (!bn_mod_sqrt_normalize(sqrt, p, ctx))
685 goto err;
686
687 verify:
688 if (!bn_mod_sqrt_verify(a_mod_p, sqrt, p, ctx))
689 goto err;
690
691 if (!bn_copy(out_sqrt, sqrt))
692 goto err;
693
694 ret = 1;
695
696 err:
697 BN_CTX_end(ctx);
698
699 return ret;
700 }
701
702 BIGNUM *
BN_mod_sqrt(BIGNUM * in,const BIGNUM * a,const BIGNUM * p,BN_CTX * ctx)703 BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
704 {
705 BIGNUM *out_sqrt;
706
707 if ((out_sqrt = in) == NULL)
708 out_sqrt = BN_new();
709 if (out_sqrt == NULL)
710 goto err;
711
712 if (!bn_mod_sqrt_internal(out_sqrt, a, p, ctx))
713 goto err;
714
715 return out_sqrt;
716
717 err:
718 if (out_sqrt != in)
719 BN_free(out_sqrt);
720
721 return NULL;
722 }
723 LCRYPTO_ALIAS(BN_mod_sqrt);
724