summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authortb <>2023-04-11 10:08:44 +0000
committertb <>2023-04-11 10:08:44 +0000
commit2389f1d8a42806b852cf81082f3eb70ecbfdd8ae (patch)
tree1a8331f19d48335fd135b93854c36f2ea3244918 /src
parented39abf2dcf0e5422ecdb755974b914bf9225ed6 (diff)
downloadopenbsd-2389f1d8a42806b852cf81082f3eb70ecbfdd8ae.tar.gz
openbsd-2389f1d8a42806b852cf81082f3eb70ecbfdd8ae.tar.bz2
openbsd-2389f1d8a42806b852cf81082f3eb70ecbfdd8ae.zip
Add a new implementation of BN_mod_sqrt()
This is a reimplementation from scratch of the Tonelli-Shanks algorithm based on Henri Cohen "A Course in Computational Algebraic Number Theory", Springer GTM 138, section 1.5.1. It is API compatible with the previous implementation, so no documentation change is required. Contrary to the old implementation, this does not have any infinite loops and has various additional sanity checks to prevent misbehavior in case the input modulus is not a prime. It contains extensive comments and the individual parts of the algorithm are split into digestible chunks instead of having one huge function. One difference of note is that it BN_mod_sqrt() now always returns the smaller of the two possible answers. In other words, while its core is non-deterministic, its answer is not. ok jsing
Diffstat (limited to 'src')
-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-03-22 20:15:14 +0000