summaryrefslogtreecommitdiff
path: root/src/lib/libc
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--src/lib/libcrypto/Makefile4
-rw-r--r--src/lib/libcrypto/bn/bn_mod_sqrt.c726
-rw-r--r--src/lib/libcrypto/bn/bn_sqrt.c409
3 files changed, 728 insertions, 411 deletions
diff --git a/src/lib/libcrypto/Makefile b/src/lib/libcrypto/Makefile
index b9a251a1d9..5405d79449 100644
--- a/src/lib/libcrypto/Makefile
+++ b/src/lib/libcrypto/Makefile
@@ -1,4 +1,4 @@
1# $OpenBSD: Makefile,v 1.99 2023/03/01 11:28:30 tb Exp $ 1# $OpenBSD: Makefile,v 1.100 2023/04/11 10:08:44 tb Exp $
2 2
3LIB= crypto 3LIB= crypto
4LIBREBUILD=y 4LIBREBUILD=y
@@ -191,6 +191,7 @@ SRCS+= bn_isqrt.c
191SRCS+= bn_kron.c 191SRCS+= bn_kron.c
192SRCS+= bn_lib.c 192SRCS+= bn_lib.c
193SRCS+= bn_mod.c 193SRCS+= bn_mod.c
194SRCS+= bn_mod_sqrt.c
194SRCS+= bn_mont.c 195SRCS+= bn_mont.c
195SRCS+= bn_mpi.c 196SRCS+= bn_mpi.c
196SRCS+= bn_mul.c 197SRCS+= bn_mul.c
@@ -202,7 +203,6 @@ SRCS+= bn_recp.c
202SRCS+= bn_shift.c 203SRCS+= bn_shift.c
203SRCS+= bn_small_primes.c 204SRCS+= bn_small_primes.c
204SRCS+= bn_sqr.c 205SRCS+= bn_sqr.c
205SRCS+= bn_sqrt.c
206SRCS+= bn_word.c 206SRCS+= bn_word.c
207SRCS+= bn_x931p.c 207SRCS+= bn_x931p.c
208 208
diff --git a/src/lib/libcrypto/bn/bn_mod_sqrt.c b/src/lib/libcrypto/bn/bn_mod_sqrt.c
new file mode 100644
index 0000000000..acca540e25
--- /dev/null
+++ b/src/lib/libcrypto/bn/bn_mod_sqrt.c
@@ -0,0 +1,726 @@
1/* $OpenBSD: bn_mod_sqrt.c,v 1.1 2023/04/11 10:08:44 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
86static int
87bn_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
112static int
113bn_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
146static int
147bn_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
208static int
209bn_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
236static int
237bn_mod_sqrt_find_sylow_generator(BIGNUM *out_generator, const BIGNUM *p,
238 const BIGNUM *q, BN_CTX *ctx)
239{
240 BIGNUM *n, *p_abs, *thirty_two;
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 ((thirty_two = BN_CTX_get(ctx)) == NULL)
249 goto err;
250 if ((p_abs = BN_CTX_get(ctx)) == NULL)
251 goto err;
252
253 for (i = 2; i < 32; i++) {
254 if (!BN_set_word(n, i))
255 goto err;
256 if (!bn_mod_sqrt_n_is_non_residue(&is_non_residue, n, p, ctx))
257 goto err;
258 if (is_non_residue)
259 goto found;
260 }
261
262 if (!BN_set_word(thirty_two, 32))
263 goto err;
264 if (!bn_copy(p_abs, p))
265 goto err;
266 BN_set_negative(p_abs, 0);
267
268 for (i = 0; i < 128; i++) {
269 if (!bn_rand_interval(n, thirty_two, p_abs))
270 goto err;
271 if (!bn_mod_sqrt_n_is_non_residue(&is_non_residue, n, p, ctx))
272 goto err;
273 if (is_non_residue)
274 goto found;
275 }
276
277 /*
278 * The probability to get here is < 2^(-128) for prime p. For squares
279 * it is easy: for p = 1369 = 37^2 this happens in ~3% of runs.
280 */
281
282 BNerror(BN_R_TOO_MANY_ITERATIONS);
283 goto err;
284
285 found:
286 /*
287 * If p is prime, n^q generates the 2-Sylow subgroup S of GF(p)*.
288 */
289
290 if (!BN_mod_exp_ct(out_generator, n, q, p, ctx))
291 goto err;
292
293 /* Sanity: p is not necessarily prime, so we could have found 0 or 1. */
294 if (BN_is_zero(out_generator) || BN_is_one(out_generator)) {
295 BNerror(BN_R_P_IS_NOT_PRIME);
296 goto err;
297 }
298
299 ret = 1;
300
301 err:
302 BN_CTX_end(ctx);
303
304 return ret;
305}
306
307/*
308 * Initialization step for Tonelli-Shanks.
309 *
310 * In the end, b = a^q (mod p) and x = a^[(q+1)/2] (mod p). Cohen optimizes this
311 * to minimize taking powers of a. This is a bit confusing and distracting, so
312 * factor this into a separate function.
313 */
314
315static int
316bn_mod_sqrt_tonelli_shanks_initialize(BIGNUM *b, BIGNUM *x, const BIGNUM *a,
317 const BIGNUM *p, const BIGNUM *q, BN_CTX *ctx)
318{
319 BIGNUM *k;
320 int ret = 0;
321
322 BN_CTX_start(ctx);
323
324 if ((k = BN_CTX_get(ctx)) == NULL)
325 goto err;
326
327 /* k = (q-1)/2. Since q is odd, we can shift. */
328 if (!BN_rshift1(k, q))
329 goto err;
330
331 /* x = a^[(q-1)/2] (mod p). */
332 if (!BN_mod_exp_ct(x, a, k, p, ctx))
333 goto err;
334
335 /* b = ax^2 = a^q (mod p). */
336 if (!BN_mod_sqr(b, x, p, ctx))
337 goto err;
338 if (!BN_mod_mul(b, a, b, p, ctx))
339 goto err;
340
341 /* x = ax = a^[(q+1)/2] (mod p). */
342 if (!BN_mod_mul(x, a, x, p, ctx))
343 goto err;
344
345 ret = 1;
346
347 err:
348 BN_CTX_end(ctx);
349
350 return ret;
351}
352
353/*
354 * Find smallest exponent m such that b^(2^m) = 1 (mod p). Assuming that a
355 * is a quadratic residue and p is a prime, we know that 1 <= m < r.
356 */
357
358static int
359bn_mod_sqrt_tonelli_shanks_find_exponent(int *out_exponent, const BIGNUM *b,
360 const BIGNUM *p, int r, BN_CTX *ctx)
361{
362 BIGNUM *x;
363 int m;
364 int ret = 0;
365
366 BN_CTX_start(ctx);
367
368 if ((x = BN_CTX_get(ctx)) == NULL)
369 goto err;
370
371 /*
372 * If r <= 1, the Tonelli-Shanks iteration should have terminated as
373 * r == 1 implies b == 1.
374 */
375 if (r <= 1) {
376 BNerror(BN_R_P_IS_NOT_PRIME);
377 goto err;
378 }
379
380 /*
381 * Sanity check to ensure taking squares actually does something:
382 * If b is 1, the Tonelli-Shanks iteration should have terminated.
383 * If b is 0, something's very wrong, in particular p can't be prime.
384 */
385 if (BN_is_zero(b) || BN_is_one(b)) {
386 BNerror(BN_R_P_IS_NOT_PRIME);
387 goto err;
388 }
389
390 if (!bn_copy(x, b))
391 goto err;
392
393 for (m = 1; m < r; m++) {
394 if (!BN_mod_sqr(x, x, p, ctx))
395 goto err;
396 if (BN_is_one(x))
397 break;
398 }
399
400 if (m >= r) {
401 /* This means a is not a quadratic residue. As (a/p) = 1, ... */
402 BNerror(BN_R_P_IS_NOT_PRIME);
403 goto err;
404 }
405
406 *out_exponent = m;
407
408 ret = 1;
409
410 err:
411 BN_CTX_end(ctx);
412
413 return ret;
414}
415
416/*
417 * The update step. With the minimal m such that b^(2^m) = 1 (mod m),
418 * set t = y^[2^(r-m-1)] (mod p) and update x = xt, y = t^2, b = by.
419 * This preserves the loop invariants a b = x^2, y^[2^(r-1)] = -1 and
420 * b^[2^(r-1)] = 1.
421 */
422
423static int
424bn_mod_sqrt_tonelli_shanks_update(BIGNUM *b, BIGNUM *x, BIGNUM *y,
425 const BIGNUM *p, int m, int r, BN_CTX *ctx)
426{
427 BIGNUM *t;
428 int ret = 0;
429
430 BN_CTX_start(ctx);
431
432 if ((t = BN_CTX_get(ctx)) == NULL)
433 goto err;
434
435 /* t = y^[2^(r-m-1)] (mod p). */
436 if (!BN_set_bit(t, r - m - 1))
437 goto err;
438 if (!BN_mod_exp_ct(t, y, t, p, ctx))
439 goto err;
440
441 /* x = xt (mod p). */
442 if (!BN_mod_mul(x, x, t, p, ctx))
443 goto err;
444
445 /* y = t^2 = y^[2^(r-m)] (mod p). */
446 if (!BN_mod_sqr(y, t, p, ctx))
447 goto err;
448
449 /* b = by (mod p). */
450 if (!BN_mod_mul(b, b, y, p, ctx))
451 goto err;
452
453 ret = 1;
454
455 err:
456 BN_CTX_end(ctx);
457
458 return ret;
459}
460
461static int
462bn_mod_sqrt_p_is_1_mod_8(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
463 BN_CTX *ctx)
464{
465 BIGNUM *b, *q, *x, *y;
466 int e, m, r;
467 int ret = 0;
468
469 BN_CTX_start(ctx);
470
471 if ((b = BN_CTX_get(ctx)) == NULL)
472 goto err;
473 if ((q = BN_CTX_get(ctx)) == NULL)
474 goto err;
475 if ((x = BN_CTX_get(ctx)) == NULL)
476 goto err;
477 if ((y = BN_CTX_get(ctx)) == NULL)
478 goto err;
479
480 /*
481 * Factor p - 1 = 2^e q with odd q. Since p = 1 (mod 8), we know e >= 3.
482 */
483
484 e = 1;
485 while (!BN_is_bit_set(p, e))
486 e++;
487 if (!BN_rshift(q, p, e))
488 goto err;
489
490 if (!bn_mod_sqrt_find_sylow_generator(y, p, q, ctx))
491 goto err;
492
493 /*
494 * Set b = a^q (mod p) and x = a^[(q+1)/2] (mod p).
495 */
496 if (!bn_mod_sqrt_tonelli_shanks_initialize(b, x, a, p, q, ctx))
497 goto err;
498
499 /*
500 * The Tonelli-Shanks iteration. Starting with r = e, the following loop
501 * invariants hold at the start of the loop.
502 *
503 * a b = x^2 (mod p)
504 * y^[2^(r-1)] = -1 (mod p)
505 * b^[2^(r-1)] = 1 (mod p)
506 *
507 * In particular, if b = 1 (mod p), x is a square root of a.
508 *
509 * Since p - 1 = 2^e q, we have 2^(e-1) q = (p - 1) / 2, so in the first
510 * iteration this follows from (a/p) = 1, (n/p) = -1, y = n^q, b = a^q.
511 *
512 * In subsequent iterations, t = y^[2^(r-m-1)], where m is the smallest
513 * m such that b^(2^m) = 1. With x = xt (mod p) and b = bt^2 (mod p) the
514 * first invariant is preserved, the second and third follow from
515 * y = t^2 (mod p) and r = m as well as the choice of m.
516 *
517 * Finally, r is strictly decreasing in each iteration. If p is prime,
518 * let S be the 2-Sylow subgroup of GF(p)*. We can prove the algorithm
519 * stops: Let S_r be the subgroup of S consisting of elements of order
520 * dividing 2^r. Then S_r = <y> and b is in S_(r-1). The S_r form a
521 * descending filtration of S and when r = 1, then b = 1.
522 */
523
524 for (r = e; r >= 1; r = m) {
525 /*
526 * Termination condition. If b == 1 then x is a square root.
527 */
528 if (BN_is_one(b))
529 goto done;
530
531 /* Find smallest exponent 1 <= m < r such that b^(2^m) == 1. */
532 if (!bn_mod_sqrt_tonelli_shanks_find_exponent(&m, b, p, r, ctx))
533 goto err;
534
535 /*
536 * With t = y^[2^(r-m-1)], update x = xt, y = t^2, b = by.
537 */
538 if (!bn_mod_sqrt_tonelli_shanks_update(b, x, y, p, m, r, ctx))
539 goto err;
540
541 /*
542 * Sanity check to make sure we don't loop indefinitely.
543 * bn_mod_sqrt_tonelli_shanks_find_exponent() ensures m < r.
544 */
545 if (r <= m)
546 goto err;
547 }
548
549 /*
550 * If p is prime, we should not get here.
551 */
552
553 BNerror(BN_R_NOT_A_SQUARE);
554 goto err;
555
556 done:
557 if (!bn_copy(out_sqrt, x))
558 goto err;
559
560 ret = 1;
561
562 err:
563 BN_CTX_end(ctx);
564
565 return ret;
566}
567
568/*
569 * Choose the smaller of sqrt and |p| - sqrt.
570 */
571
572static int
573bn_mod_sqrt_normalize(BIGNUM *sqrt, const BIGNUM *p, BN_CTX *ctx)
574{
575 BIGNUM *x;
576 int ret = 0;
577
578 BN_CTX_start(ctx);
579
580 if ((x = BN_CTX_get(ctx)) == NULL)
581 goto err;
582
583 if (!BN_lshift1(x, sqrt))
584 goto err;
585
586 if (BN_ucmp(x, p) > 0) {
587 if (!BN_usub(sqrt, p, sqrt))
588 goto err;
589 }
590
591 ret = 1;
592
593 err:
594 BN_CTX_end(ctx);
595
596 return ret;
597}
598
599/*
600 * Verify that a = (sqrt_a)^2 (mod p). Requires that a is reduced (mod p).
601 */
602
603static int
604bn_mod_sqrt_verify(const BIGNUM *a, const BIGNUM *sqrt_a, const BIGNUM *p,
605 BN_CTX *ctx)
606{
607 BIGNUM *x;
608 int ret = 0;
609
610 BN_CTX_start(ctx);
611
612 if ((x = BN_CTX_get(ctx)) == NULL)
613 goto err;
614
615 if (!BN_mod_sqr(x, sqrt_a, p, ctx))
616 goto err;
617
618 if (BN_cmp(x, a) != 0) {
619 BNerror(BN_R_NOT_A_SQUARE);
620 goto err;
621 }
622
623 ret = 1;
624
625 err:
626 BN_CTX_end(ctx);
627
628 return ret;
629}
630
631static int
632bn_mod_sqrt_internal(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
633 BN_CTX *ctx)
634{
635 BIGNUM *a_mod_p, *sqrt;
636 BN_ULONG lsw;
637 int done;
638 int kronecker;
639 int ret = 0;
640
641 BN_CTX_start(ctx);
642
643 if ((a_mod_p = BN_CTX_get(ctx)) == NULL)
644 goto err;
645 if ((sqrt = BN_CTX_get(ctx)) == NULL)
646 goto err;
647
648 if (!BN_nnmod(a_mod_p, a, p, ctx))
649 goto err;
650
651 if (!bn_mod_sqrt_trivial_cases(&done, sqrt, a_mod_p, p, ctx))
652 goto err;
653 if (done)
654 goto verify;
655
656 /*
657 * Make sure that the Kronecker symbol (a/p) == 1. In case p is prime
658 * this is equivalent to a having a square root (mod p). The cost of
659 * BN_kronecker() is O(log^2(n)). This is small compared to the cost
660 * O(log^4(n)) of Tonelli-Shanks.
661 */
662
663 if ((kronecker = BN_kronecker(a_mod_p, p, ctx)) == -2)
664 goto err;
665 if (kronecker <= 0) {
666 /* This error is only accurate if p is known to be a prime. */
667 BNerror(BN_R_NOT_A_SQUARE);
668 goto err;
669 }
670
671 lsw = BN_lsw(p);
672
673 if (lsw % 4 == 3) {
674 if (!bn_mod_sqrt_p_is_3_mod_4(sqrt, a_mod_p, p, ctx))
675 goto err;
676 } else if (lsw % 8 == 5) {
677 if (!bn_mod_sqrt_p_is_5_mod_8(sqrt, a_mod_p, p, ctx))
678 goto err;
679 } else if (lsw % 8 == 1) {
680 if (!bn_mod_sqrt_p_is_1_mod_8(sqrt, a_mod_p, p, ctx))
681 goto err;
682 } else {
683 /* Impossible to hit since the trivial cases ensure p is odd. */
684 BNerror(BN_R_P_IS_NOT_PRIME);
685 goto err;
686 }
687
688 if (!bn_mod_sqrt_normalize(sqrt, p, ctx))
689 goto err;
690
691 verify:
692 if (!bn_mod_sqrt_verify(a_mod_p, sqrt, p, ctx))
693 goto err;
694
695 if (!bn_copy(out_sqrt, sqrt))
696 goto err;
697
698 ret = 1;
699
700 err:
701 BN_CTX_end(ctx);
702
703 return ret;
704}
705
706BIGNUM *
707BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
708{
709 BIGNUM *out_sqrt;
710
711 if ((out_sqrt = in) == NULL)
712 out_sqrt = BN_new();
713 if (out_sqrt == NULL)
714 goto err;
715
716 if (!bn_mod_sqrt_internal(out_sqrt, a, p, ctx))
717 goto err;
718
719 return out_sqrt;
720
721 err:
722 if (out_sqrt != in)
723 BN_free(out_sqrt);
724
725 return NULL;
726}
diff --git a/src/lib/libcrypto/bn/bn_sqrt.c b/src/lib/libcrypto/bn/bn_sqrt.c
deleted file mode 100644
index 3d9f017f59..0000000000
--- a/src/lib/libcrypto/bn/bn_sqrt.c
+++ /dev/null
@@ -1,409 +0,0 @@
1/* $OpenBSD: bn_sqrt.c,v 1.16 2023/03/27 10:25:02 tb Exp $ */
2/* Written by Lenka Fibikova <fibikova@exp-math.uni-essen.de>
3 * and Bodo Moeller for the OpenSSL project. */
4/* ====================================================================
5 * Copyright (c) 1998-2000 The OpenSSL Project. All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 *
11 * 1. Redistributions of source code must retain the above copyright
12 * notice, this list of conditions and the following disclaimer.
13 *
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in
16 * the documentation and/or other materials provided with the
17 * distribution.
18 *
19 * 3. All advertising materials mentioning features or use of this
20 * software must display the following acknowledgment:
21 * "This product includes software developed by the OpenSSL Project
22 * for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
23 *
24 * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
25 * endorse or promote products derived from this software without
26 * prior written permission. For written permission, please contact
27 * openssl-core@openssl.org.
28 *
29 * 5. Products derived from this software may not be called "OpenSSL"
30 * nor may "OpenSSL" appear in their names without prior written
31 * permission of the OpenSSL Project.
32 *
33 * 6. Redistributions of any form whatsoever must retain the following
34 * acknowledgment:
35 * "This product includes software developed by the OpenSSL Project
36 * for use in the OpenSSL Toolkit (http://www.openssl.org/)"
37 *
38 * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
39 * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
40 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
41 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE OpenSSL PROJECT OR
42 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
43 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
44 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
45 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
46 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
47 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
48 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
49 * OF THE POSSIBILITY OF SUCH DAMAGE.
50 * ====================================================================
51 *
52 * This product includes cryptographic software written by Eric Young
53 * (eay@cryptsoft.com). This product includes software written by Tim
54 * Hudson (tjh@cryptsoft.com).
55 *
56 */
57
58#include <openssl/err.h>
59
60#include "bn_local.h"
61
62/*
63 * Returns 'ret' such that ret^2 == a (mod p), if it exists, using the
64 * Tonelli-Shanks algorithm following Henri Cohen, "A Course in Computational
65 * Algebraic Number Theory", algorithm 1.5.1, Springer, Berlin, 1996.
66 *
67 * Note: 'p' must be prime!
68 */
69
70BIGNUM *
71BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
72{
73 BIGNUM *ret = in;
74 int err = 1;
75 int r;
76 BIGNUM *A, *b, *q, *t, *x, *y;
77 int e, i, j;
78
79 if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) {
80 if (BN_abs_is_word(p, 2)) {
81 if (ret == NULL)
82 ret = BN_new();
83 if (ret == NULL)
84 goto end;
85 if (!BN_set_word(ret, BN_is_bit_set(a, 0))) {
86 if (ret != in)
87 BN_free(ret);
88 return NULL;
89 }
90 return ret;
91 }
92
93 BNerror(BN_R_P_IS_NOT_PRIME);
94 return (NULL);
95 }
96
97 if (BN_is_zero(a) || BN_is_one(a)) {
98 if (ret == NULL)
99 ret = BN_new();
100 if (ret == NULL)
101 goto end;
102 if (!BN_set_word(ret, BN_is_one(a))) {
103 if (ret != in)
104 BN_free(ret);
105 return NULL;
106 }
107 return ret;
108 }
109
110 BN_CTX_start(ctx);
111 if ((A = BN_CTX_get(ctx)) == NULL)
112 goto end;
113 if ((b = BN_CTX_get(ctx)) == NULL)
114 goto end;
115 if ((q = BN_CTX_get(ctx)) == NULL)
116 goto end;
117 if ((t = BN_CTX_get(ctx)) == NULL)
118 goto end;
119 if ((x = BN_CTX_get(ctx)) == NULL)
120 goto end;
121 if ((y = BN_CTX_get(ctx)) == NULL)
122 goto end;
123
124 if (ret == NULL)
125 ret = BN_new();
126 if (ret == NULL)
127 goto end;
128
129 /* A = a mod p */
130 if (!BN_nnmod(A, a, p, ctx))
131 goto end;
132
133 /* now write |p| - 1 as 2^e*q where q is odd */
134 e = 1;
135 while (!BN_is_bit_set(p, e))
136 e++;
137 /* we'll set q later (if needed) */
138
139 if (e == 1) {
140 /* The easy case: (|p|-1)/2 is odd, so 2 has an inverse
141 * modulo (|p|-1)/2, and square roots can be computed
142 * directly by modular exponentiation.
143 * We have
144 * 2 * (|p|+1)/4 == 1 (mod (|p|-1)/2),
145 * so we can use exponent (|p|+1)/4, i.e. (|p|-3)/4 + 1.
146 */
147 if (!BN_rshift(q, p, 2))
148 goto end;
149 q->neg = 0;
150 if (!BN_add_word(q, 1))
151 goto end;
152 if (!BN_mod_exp_ct(ret, A, q, p, ctx))
153 goto end;
154 err = 0;
155 goto vrfy;
156 }
157
158 if (e == 2) {
159 /* |p| == 5 (mod 8)
160 *
161 * In this case 2 is always a non-square since
162 * Legendre(2,p) = (-1)^((p^2-1)/8) for any odd prime.
163 * So if a really is a square, then 2*a is a non-square.
164 * Thus for
165 * b := (2*a)^((|p|-5)/8),
166 * i := (2*a)*b^2
167 * we have
168 * i^2 = (2*a)^((1 + (|p|-5)/4)*2)
169 * = (2*a)^((p-1)/2)
170 * = -1;
171 * so if we set
172 * x := a*b*(i-1),
173 * then
174 * x^2 = a^2 * b^2 * (i^2 - 2*i + 1)
175 * = a^2 * b^2 * (-2*i)
176 * = a*(-i)*(2*a*b^2)
177 * = a*(-i)*i
178 * = a.
179 *
180 * (This is due to A.O.L. Atkin,
181 * <URL: http://listserv.nodak.edu/scripts/wa.exe?A2=ind9211&L=nmbrthry&O=T&P=562>,
182 * November 1992.)
183 */
184
185 /* t := 2*a */
186 if (!BN_mod_lshift1_quick(t, A, p))
187 goto end;
188
189 /* b := (2*a)^((|p|-5)/8) */
190 if (!BN_rshift(q, p, 3))
191 goto end;
192 q->neg = 0;
193 if (!BN_mod_exp_ct(b, t, q, p, ctx))
194 goto end;
195
196 /* y := b^2 */
197 if (!BN_mod_sqr(y, b, p, ctx))
198 goto end;
199
200 /* t := (2*a)*b^2 - 1*/
201 if (!BN_mod_mul(t, t, y, p, ctx))
202 goto end;
203 if (!BN_sub_word(t, 1))
204 goto end;
205
206 /* x = a*b*t */
207 if (!BN_mod_mul(x, A, b, p, ctx))
208 goto end;
209 if (!BN_mod_mul(x, x, t, p, ctx))
210 goto end;
211
212 if (!bn_copy(ret, x))
213 goto end;
214 err = 0;
215 goto vrfy;
216 }
217
218 /* e > 2, so we really have to use the Tonelli/Shanks algorithm.
219 * First, find some y that is not a square. */
220 if (!bn_copy(q, p)) /* use 'q' as temp */
221 goto end;
222 q->neg = 0;
223 i = 2;
224 do {
225 /* For efficiency, try small numbers first;
226 * if this fails, try random numbers.
227 */
228 if (i < 22) {
229 if (!BN_set_word(y, i))
230 goto end;
231 } else {
232 if (!BN_pseudo_rand(y, BN_num_bits(p), 0, 0))
233 goto end;
234 if (BN_ucmp(y, p) >= 0) {
235 if (p->neg) {
236 if (!BN_add(y, y, p))
237 goto end;
238 } else {
239 if (!BN_sub(y, y, p))
240 goto end;
241 }
242 }
243 /* now 0 <= y < |p| */
244 if (BN_is_zero(y))
245 if (!BN_set_word(y, i))
246 goto end;
247 }
248
249 r = BN_kronecker(y, q, ctx); /* here 'q' is |p| */
250 if (r < -1)
251 goto end;
252 if (r == 0) {
253 /* m divides p */
254 BNerror(BN_R_P_IS_NOT_PRIME);
255 goto end;
256 }
257 } while (r == 1 && ++i < 82);
258
259 if (r != -1) {
260 /* Many rounds and still no non-square -- this is more likely
261 * a bug than just bad luck.
262 * Even if p is not prime, we should have found some y
263 * such that r == -1.
264 */
265 BNerror(BN_R_TOO_MANY_ITERATIONS);
266 goto end;
267 }
268
269 /* Here's our actual 'q': */
270 if (!BN_rshift(q, q, e))
271 goto end;
272
273 /* Now that we have some non-square, we can find an element
274 * of order 2^e by computing its q'th power. */
275 if (!BN_mod_exp_ct(y, y, q, p, ctx))
276 goto end;
277 if (BN_is_one(y)) {
278 BNerror(BN_R_P_IS_NOT_PRIME);
279 goto end;
280 }
281
282 /* Now we know that (if p is indeed prime) there is an integer
283 * k, 0 <= k < 2^e, such that
284 *
285 * a^q * y^k == 1 (mod p).
286 *
287 * As a^q is a square and y is not, k must be even.
288 * q+1 is even, too, so there is an element
289 *
290 * X := a^((q+1)/2) * y^(k/2),
291 *
292 * and it satisfies
293 *
294 * X^2 = a^q * a * y^k
295 * = a,
296 *
297 * so it is the square root that we are looking for.
298 */
299
300 /* t := (q-1)/2 (note that q is odd) */
301 if (!BN_rshift1(t, q))
302 goto end;
303
304 /* x := a^((q-1)/2) */
305 if (BN_is_zero(t)) { /* special case: p = 2^e + 1 */
306 if (!BN_nnmod(t, A, p, ctx))
307 goto end;
308 if (BN_is_zero(t)) {
309 /* special case: a == 0 (mod p) */
310 BN_zero(ret);
311 err = 0;
312 goto end;
313 } else if (!BN_one(x))
314 goto end;
315 } else {
316 if (!BN_mod_exp_ct(x, A, t, p, ctx))
317 goto end;
318 if (BN_is_zero(x)) {
319 /* special case: a == 0 (mod p) */
320 BN_zero(ret);
321 err = 0;
322 goto end;
323 }
324 }
325
326 /* b := a*x^2 (= a^q) */
327 if (!BN_mod_sqr(b, x, p, ctx))
328 goto end;
329 if (!BN_mod_mul(b, b, A, p, ctx))
330 goto end;
331
332 /* x := a*x (= a^((q+1)/2)) */
333 if (!BN_mod_mul(x, x, A, p, ctx))
334 goto end;
335
336 while (1) {
337 /* Now b is a^q * y^k for some even k (0 <= k < 2^E
338 * where E refers to the original value of e, which we
339 * don't keep in a variable), and x is a^((q+1)/2) * y^(k/2).
340 *
341 * We have a*b = x^2,
342 * y^2^(e-1) = -1,
343 * b^2^(e-1) = 1.
344 */
345
346 if (BN_is_one(b)) {
347 if (!bn_copy(ret, x))
348 goto end;
349 err = 0;
350 goto vrfy;
351 }
352
353 /* Find the smallest i with 0 < i < e such that b^(2^i) = 1. */
354 for (i = 1; i < e; i++) {
355 if (i == 1) {
356 if (!BN_mod_sqr(t, b, p, ctx))
357 goto end;
358 } else {
359 if (!BN_mod_sqr(t, t, p, ctx))
360 goto end;
361 }
362 if (BN_is_one(t))
363 break;
364 }
365 if (i >= e) {
366 BNerror(BN_R_NOT_A_SQUARE);
367 goto end;
368 }
369
370 /* t := y^2^(e - i - 1) */
371 if (!bn_copy(t, y))
372 goto end;
373 for (j = e - i - 1; j > 0; j--) {
374 if (!BN_mod_sqr(t, t, p, ctx))
375 goto end;
376 }
377 if (!BN_mod_mul(y, t, t, p, ctx))
378 goto end;
379 if (!BN_mod_mul(x, x, t, p, ctx))
380 goto end;
381 if (!BN_mod_mul(b, b, y, p, ctx))
382 goto end;
383 e = i;
384 }
385
386vrfy:
387 if (!err) {
388 /* verify the result -- the input might have been not a square
389 * (test added in 0.9.8) */
390
391 if (!BN_mod_sqr(x, ret, p, ctx))
392 err = 1;
393
394 if (!err && 0 != BN_cmp(x, A)) {
395 BNerror(BN_R_NOT_A_SQUARE);
396 err = 1;
397 }
398 }
399
400end:
401 if (err) {
402 if (ret != NULL && ret != in) {
403 BN_free(ret);
404 }
405 ret = NULL;
406 }
407 BN_CTX_end(ctx);
408 return ret;
409}
generated by cgit v1.2.3-55-g6feb (git 2.39.1) at 2025-10-28 18:27:16 +0000