xref: /openbsd/lib/libcrypto/bn/bn_gcd.c (revision 3bef86f7)
1 /* $OpenBSD: bn_gcd.c,v 1.28 2023/06/02 17:15:30 tb Exp $ */
2 /* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
3  * All rights reserved.
4  *
5  * This package is an SSL implementation written
6  * by Eric Young (eay@cryptsoft.com).
7  * The implementation was written so as to conform with Netscapes SSL.
8  *
9  * This library is free for commercial and non-commercial use as long as
10  * the following conditions are aheared to.  The following conditions
11  * apply to all code found in this distribution, be it the RC4, RSA,
12  * lhash, DES, etc., code; not just the SSL code.  The SSL documentation
13  * included with this distribution is covered by the same copyright terms
14  * except that the holder is Tim Hudson (tjh@cryptsoft.com).
15  *
16  * Copyright remains Eric Young's, and as such any Copyright notices in
17  * the code are not to be removed.
18  * If this package is used in a product, Eric Young should be given attribution
19  * as the author of the parts of the library used.
20  * This can be in the form of a textual message at program startup or
21  * in documentation (online or textual) provided with the package.
22  *
23  * Redistribution and use in source and binary forms, with or without
24  * modification, are permitted provided that the following conditions
25  * are met:
26  * 1. Redistributions of source code must retain the copyright
27  *    notice, this list of conditions and the following disclaimer.
28  * 2. Redistributions in binary form must reproduce the above copyright
29  *    notice, this list of conditions and the following disclaimer in the
30  *    documentation and/or other materials provided with the distribution.
31  * 3. All advertising materials mentioning features or use of this software
32  *    must display the following acknowledgement:
33  *    "This product includes cryptographic software written by
34  *     Eric Young (eay@cryptsoft.com)"
35  *    The word 'cryptographic' can be left out if the rouines from the library
36  *    being used are not cryptographic related :-).
37  * 4. If you include any Windows specific code (or a derivative thereof) from
38  *    the apps directory (application code) you must include an acknowledgement:
39  *    "This product includes software written by Tim Hudson (tjh@cryptsoft.com)"
40  *
41  * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``AS IS'' AND
42  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
43  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
44  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
45  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
46  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
47  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
48  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
49  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
50  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
51  * SUCH DAMAGE.
52  *
53  * The licence and distribution terms for any publically available version or
54  * derivative of this code cannot be changed.  i.e. this code cannot simply be
55  * copied and put under another distribution licence
56  * [including the GNU Public Licence.]
57  */
58 /* ====================================================================
59  * Copyright (c) 1998-2001 The OpenSSL Project.  All rights reserved.
60  *
61  * Redistribution and use in source and binary forms, with or without
62  * modification, are permitted provided that the following conditions
63  * are met:
64  *
65  * 1. Redistributions of source code must retain the above copyright
66  *    notice, this list of conditions and the following disclaimer.
67  *
68  * 2. Redistributions in binary form must reproduce the above copyright
69  *    notice, this list of conditions and the following disclaimer in
70  *    the documentation and/or other materials provided with the
71  *    distribution.
72  *
73  * 3. All advertising materials mentioning features or use of this
74  *    software must display the following acknowledgment:
75  *    "This product includes software developed by the OpenSSL Project
76  *    for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
77  *
78  * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
79  *    endorse or promote products derived from this software without
80  *    prior written permission. For written permission, please contact
81  *    openssl-core@openssl.org.
82  *
83  * 5. Products derived from this software may not be called "OpenSSL"
84  *    nor may "OpenSSL" appear in their names without prior written
85  *    permission of the OpenSSL Project.
86  *
87  * 6. Redistributions of any form whatsoever must retain the following
88  *    acknowledgment:
89  *    "This product includes software developed by the OpenSSL Project
90  *    for use in the OpenSSL Toolkit (http://www.openssl.org/)"
91  *
92  * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
93  * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
94  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
95  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE OpenSSL PROJECT OR
96  * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
97  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
98  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
99  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
100  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
101  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
102  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
103  * OF THE POSSIBILITY OF SUCH DAMAGE.
104  * ====================================================================
105  *
106  * This product includes cryptographic software written by Eric Young
107  * (eay@cryptsoft.com).  This product includes software written by Tim
108  * Hudson (tjh@cryptsoft.com).
109  *
110  */
111 
112 #include <openssl/err.h>
113 
114 #include "bn_local.h"
115 
116 static BIGNUM *
117 euclid(BIGNUM *a, BIGNUM *b)
118 {
119 	BIGNUM *t;
120 	int shifts = 0;
121 
122 	/* Loop invariant: 0 <= b <= a. */
123 	while (!BN_is_zero(b)) {
124 		if (BN_is_odd(a) && BN_is_odd(b)) {
125 			if (!BN_sub(a, a, b))
126 				goto err;
127 			if (!BN_rshift1(a, a))
128 				goto err;
129 		} else if (BN_is_odd(a) && !BN_is_odd(b)) {
130 			if (!BN_rshift1(b, b))
131 				goto err;
132 		} else if (!BN_is_odd(a) && BN_is_odd(b)) {
133 			if (!BN_rshift1(a, a))
134 				goto err;
135 		} else {
136 			if (!BN_rshift1(a, a))
137 				goto err;
138 			if (!BN_rshift1(b, b))
139 				goto err;
140 			shifts++;
141 			continue;
142 		}
143 
144 		if (BN_cmp(a, b) < 0) {
145 			t = a;
146 			a = b;
147 			b = t;
148 		}
149 	}
150 
151 	if (shifts) {
152 		if (!BN_lshift(a, a, shifts))
153 			goto err;
154 	}
155 
156 	return a;
157 
158  err:
159 	return NULL;
160 }
161 
162 int
163 BN_gcd(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
164 {
165 	BIGNUM *a, *b, *t;
166 	int ret = 0;
167 
168 	BN_CTX_start(ctx);
169 	if ((a = BN_CTX_get(ctx)) == NULL)
170 		goto err;
171 	if ((b = BN_CTX_get(ctx)) == NULL)
172 		goto err;
173 
174 	if (!bn_copy(a, in_a))
175 		goto err;
176 	if (!bn_copy(b, in_b))
177 		goto err;
178 	a->neg = 0;
179 	b->neg = 0;
180 
181 	if (BN_cmp(a, b) < 0) {
182 		t = a;
183 		a = b;
184 		b = t;
185 	}
186 	t = euclid(a, b);
187 	if (t == NULL)
188 		goto err;
189 
190 	if (!bn_copy(r, t))
191 		goto err;
192 	ret = 1;
193 
194  err:
195 	BN_CTX_end(ctx);
196 	return (ret);
197 }
198 
199 int
200 BN_gcd_nonct(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
201 {
202 	return BN_gcd(r, in_a, in_b, ctx);
203 }
204 
205 /*
206  * BN_gcd_no_branch is a special version of BN_mod_inverse_no_branch.
207  * that returns the GCD.
208  */
209 static BIGNUM *
210 BN_gcd_no_branch(BIGNUM *in, const BIGNUM *a, const BIGNUM *n,
211     BN_CTX *ctx)
212 {
213 	BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
214 	BIGNUM local_A, local_B;
215 	BIGNUM *pA, *pB;
216 	BIGNUM *ret = NULL;
217 	int sign;
218 
219 	if (in == NULL)
220 		goto err;
221 	R = in;
222 
223 	BN_init(&local_A);
224 	BN_init(&local_B);
225 
226 	BN_CTX_start(ctx);
227 	if ((A = BN_CTX_get(ctx)) == NULL)
228 		goto err;
229 	if ((B = BN_CTX_get(ctx)) == NULL)
230 		goto err;
231 	if ((X = BN_CTX_get(ctx)) == NULL)
232 		goto err;
233 	if ((D = BN_CTX_get(ctx)) == NULL)
234 		goto err;
235 	if ((M = BN_CTX_get(ctx)) == NULL)
236 		goto err;
237 	if ((Y = BN_CTX_get(ctx)) == NULL)
238 		goto err;
239 	if ((T = BN_CTX_get(ctx)) == NULL)
240 		goto err;
241 
242 	if (!BN_one(X))
243 		goto err;
244 	BN_zero(Y);
245 	if (!bn_copy(B, a))
246 		goto err;
247 	if (!bn_copy(A, n))
248 		goto err;
249 	A->neg = 0;
250 
251 	if (B->neg || (BN_ucmp(B, A) >= 0)) {
252 		/*
253 		 * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
254 		 * BN_div_no_branch will be called eventually.
255 		 */
256 		pB = &local_B;
257 		/* BN_init() done at the top of the function. */
258 		BN_with_flags(pB, B, BN_FLG_CONSTTIME);
259 		if (!BN_nnmod(B, pB, A, ctx))
260 			goto err;
261 	}
262 	sign = -1;
263 	/* From  B = a mod |n|,  A = |n|  it follows that
264 	 *
265 	 *      0 <= B < A,
266 	 *     -sign*X*a  ==  B   (mod |n|),
267 	 *      sign*Y*a  ==  A   (mod |n|).
268 	 */
269 
270 	while (!BN_is_zero(B)) {
271 		BIGNUM *tmp;
272 
273 		/*
274 		 *      0 < B < A,
275 		 * (*) -sign*X*a  ==  B   (mod |n|),
276 		 *      sign*Y*a  ==  A   (mod |n|)
277 		 */
278 
279 		/*
280 		 * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
281 		 * BN_div_no_branch will be called eventually.
282 		 */
283 		pA = &local_A;
284 		/* BN_init() done at the top of the function. */
285 		BN_with_flags(pA, A, BN_FLG_CONSTTIME);
286 
287 		/* (D, M) := (A/B, A%B) ... */
288 		if (!BN_div_ct(D, M, pA, B, ctx))
289 			goto err;
290 
291 		/* Now
292 		 *      A = D*B + M;
293 		 * thus we have
294 		 * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
295 		 */
296 		tmp = A; /* keep the BIGNUM object, the value does not matter */
297 
298 		/* (A, B) := (B, A mod B) ... */
299 		A = B;
300 		B = M;
301 		/* ... so we have  0 <= B < A  again */
302 
303 		/* Since the former  M  is now  B  and the former  B  is now  A,
304 		 * (**) translates into
305 		 *       sign*Y*a  ==  D*A + B    (mod |n|),
306 		 * i.e.
307 		 *       sign*Y*a - D*A  ==  B    (mod |n|).
308 		 * Similarly, (*) translates into
309 		 *      -sign*X*a  ==  A          (mod |n|).
310 		 *
311 		 * Thus,
312 		 *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
313 		 * i.e.
314 		 *        sign*(Y + D*X)*a  ==  B  (mod |n|).
315 		 *
316 		 * So if we set  (X, Y, sign) := (Y + D*X, X, -sign),  we arrive back at
317 		 *      -sign*X*a  ==  B   (mod |n|),
318 		 *       sign*Y*a  ==  A   (mod |n|).
319 		 * Note that  X  and  Y  stay non-negative all the time.
320 		 */
321 
322 		if (!BN_mul(tmp, D, X, ctx))
323 			goto err;
324 		if (!BN_add(tmp, tmp, Y))
325 			goto err;
326 
327 		M = Y; /* keep the BIGNUM object, the value does not matter */
328 		Y = X;
329 		X = tmp;
330 		sign = -sign;
331 	}
332 
333 	/*
334 	 * The while loop (Euclid's algorithm) ends when
335 	 *      A == gcd(a,n);
336 	 */
337 
338 	if (!bn_copy(R, A))
339 		goto err;
340 	ret = R;
341  err:
342 	if ((ret == NULL) && (in == NULL))
343 		BN_free(R);
344 	BN_CTX_end(ctx);
345 	return (ret);
346 }
347 
348 int
349 BN_gcd_ct(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
350 {
351 	if (BN_gcd_no_branch(r, in_a, in_b, ctx) == NULL)
352 		return 0;
353 	return 1;
354 }
355 
356 /* BN_mod_inverse_no_branch is a special version of BN_mod_inverse.
357  * It does not contain branches that may leak sensitive information.
358  */
359 static BIGNUM *
360 BN_mod_inverse_no_branch(BIGNUM *in, const BIGNUM *a, const BIGNUM *n,
361     BN_CTX *ctx)
362 {
363 	BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
364 	BIGNUM local_A, local_B;
365 	BIGNUM *pA, *pB;
366 	BIGNUM *ret = NULL;
367 	int sign;
368 
369 	BN_init(&local_A);
370 	BN_init(&local_B);
371 
372 	BN_CTX_start(ctx);
373 	if ((A = BN_CTX_get(ctx)) == NULL)
374 		goto err;
375 	if ((B = BN_CTX_get(ctx)) == NULL)
376 		goto err;
377 	if ((X = BN_CTX_get(ctx)) == NULL)
378 		goto err;
379 	if ((D = BN_CTX_get(ctx)) == NULL)
380 		goto err;
381 	if ((M = BN_CTX_get(ctx)) == NULL)
382 		goto err;
383 	if ((Y = BN_CTX_get(ctx)) == NULL)
384 		goto err;
385 	if ((T = BN_CTX_get(ctx)) == NULL)
386 		goto err;
387 
388 	if (in == NULL)
389 		R = BN_new();
390 	else
391 		R = in;
392 	if (R == NULL)
393 		goto err;
394 
395 	if (!BN_one(X))
396 		goto err;
397 	BN_zero(Y);
398 	if (!bn_copy(B, a))
399 		goto err;
400 	if (!bn_copy(A, n))
401 		goto err;
402 	A->neg = 0;
403 
404 	if (B->neg || (BN_ucmp(B, A) >= 0)) {
405 		/*
406 		 * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
407 		 * BN_div_no_branch will be called eventually.
408 		 */
409 		pB = &local_B;
410 		/* BN_init() done at the top of the function. */
411 		BN_with_flags(pB, B, BN_FLG_CONSTTIME);
412 		if (!BN_nnmod(B, pB, A, ctx))
413 			goto err;
414 	}
415 	sign = -1;
416 	/* From  B = a mod |n|,  A = |n|  it follows that
417 	 *
418 	 *      0 <= B < A,
419 	 *     -sign*X*a  ==  B   (mod |n|),
420 	 *      sign*Y*a  ==  A   (mod |n|).
421 	 */
422 
423 	while (!BN_is_zero(B)) {
424 		BIGNUM *tmp;
425 
426 		/*
427 		 *      0 < B < A,
428 		 * (*) -sign*X*a  ==  B   (mod |n|),
429 		 *      sign*Y*a  ==  A   (mod |n|)
430 		 */
431 
432 		/*
433 		 * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
434 		 * BN_div_no_branch will be called eventually.
435 		 */
436 		pA = &local_A;
437 		/* BN_init() done at the top of the function. */
438 		BN_with_flags(pA, A, BN_FLG_CONSTTIME);
439 
440 		/* (D, M) := (A/B, A%B) ... */
441 		if (!BN_div_ct(D, M, pA, B, ctx))
442 			goto err;
443 
444 		/* Now
445 		 *      A = D*B + M;
446 		 * thus we have
447 		 * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
448 		 */
449 		tmp = A; /* keep the BIGNUM object, the value does not matter */
450 
451 		/* (A, B) := (B, A mod B) ... */
452 		A = B;
453 		B = M;
454 		/* ... so we have  0 <= B < A  again */
455 
456 		/* Since the former  M  is now  B  and the former  B  is now  A,
457 		 * (**) translates into
458 		 *       sign*Y*a  ==  D*A + B    (mod |n|),
459 		 * i.e.
460 		 *       sign*Y*a - D*A  ==  B    (mod |n|).
461 		 * Similarly, (*) translates into
462 		 *      -sign*X*a  ==  A          (mod |n|).
463 		 *
464 		 * Thus,
465 		 *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
466 		 * i.e.
467 		 *        sign*(Y + D*X)*a  ==  B  (mod |n|).
468 		 *
469 		 * So if we set  (X, Y, sign) := (Y + D*X, X, -sign),  we arrive back at
470 		 *      -sign*X*a  ==  B   (mod |n|),
471 		 *       sign*Y*a  ==  A   (mod |n|).
472 		 * Note that  X  and  Y  stay non-negative all the time.
473 		 */
474 
475 		if (!BN_mul(tmp, D, X, ctx))
476 			goto err;
477 		if (!BN_add(tmp, tmp, Y))
478 			goto err;
479 
480 		M = Y; /* keep the BIGNUM object, the value does not matter */
481 		Y = X;
482 		X = tmp;
483 		sign = -sign;
484 	}
485 
486 	/*
487 	 * The while loop (Euclid's algorithm) ends when
488 	 *      A == gcd(a,n);
489 	 * we have
490 	 *       sign*Y*a  ==  A  (mod |n|),
491 	 * where  Y  is non-negative.
492 	 */
493 
494 	if (sign < 0) {
495 		if (!BN_sub(Y, n, Y))
496 			goto err;
497 	}
498 	/* Now  Y*a  ==  A  (mod |n|).  */
499 
500 	if (!BN_is_one(A)) {
501 		BNerror(BN_R_NO_INVERSE);
502 		goto err;
503 	}
504 
505 	if (!BN_nnmod(Y, Y, n, ctx))
506 		goto err;
507 	if (!bn_copy(R, Y))
508 		goto err;
509 
510 	ret = R;
511 
512  err:
513 	if ((ret == NULL) && (in == NULL))
514 		BN_free(R);
515 	BN_CTX_end(ctx);
516 	return (ret);
517 }
518 
519 /* solves ax == 1 (mod n) */
520 static BIGNUM *
521 BN_mod_inverse_internal(BIGNUM *in, const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx,
522     int ct)
523 {
524 	BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
525 	BIGNUM *ret = NULL;
526 	int sign;
527 
528 	if (ct)
529 		return BN_mod_inverse_no_branch(in, a, n, ctx);
530 
531 	BN_CTX_start(ctx);
532 	if ((A = BN_CTX_get(ctx)) == NULL)
533 		goto err;
534 	if ((B = BN_CTX_get(ctx)) == NULL)
535 		goto err;
536 	if ((X = BN_CTX_get(ctx)) == NULL)
537 		goto err;
538 	if ((D = BN_CTX_get(ctx)) == NULL)
539 		goto err;
540 	if ((M = BN_CTX_get(ctx)) == NULL)
541 		goto err;
542 	if ((Y = BN_CTX_get(ctx)) == NULL)
543 		goto err;
544 	if ((T = BN_CTX_get(ctx)) == NULL)
545 		goto err;
546 
547 	if (in == NULL)
548 		R = BN_new();
549 	else
550 		R = in;
551 	if (R == NULL)
552 		goto err;
553 
554 	if (!BN_one(X))
555 		goto err;
556 	BN_zero(Y);
557 	if (!bn_copy(B, a))
558 		goto err;
559 	if (!bn_copy(A, n))
560 		goto err;
561 	A->neg = 0;
562 	if (B->neg || (BN_ucmp(B, A) >= 0)) {
563 		if (!BN_nnmod(B, B, A, ctx))
564 			goto err;
565 	}
566 	sign = -1;
567 	/* From  B = a mod |n|,  A = |n|  it follows that
568 	 *
569 	 *      0 <= B < A,
570 	 *     -sign*X*a  ==  B   (mod |n|),
571 	 *      sign*Y*a  ==  A   (mod |n|).
572 	 */
573 
574 	if (BN_is_odd(n) && (BN_num_bits(n) <= (BN_BITS <= 32 ? 450 : 2048))) {
575 		/* Binary inversion algorithm; requires odd modulus.
576 		 * This is faster than the general algorithm if the modulus
577 		 * is sufficiently small (about 400 .. 500 bits on 32-bit
578 		 * systems, but much more on 64-bit systems) */
579 		int shift;
580 
581 		while (!BN_is_zero(B)) {
582 			/*
583 			 *      0 < B < |n|,
584 			 *      0 < A <= |n|,
585 			 * (1) -sign*X*a  ==  B   (mod |n|),
586 			 * (2)  sign*Y*a  ==  A   (mod |n|)
587 			 */
588 
589 			/* Now divide  B  by the maximum possible power of two in the integers,
590 			 * and divide  X  by the same value mod |n|.
591 			 * When we're done, (1) still holds. */
592 			shift = 0;
593 			while (!BN_is_bit_set(B, shift)) /* note that 0 < B */
594 			{
595 				shift++;
596 
597 				if (BN_is_odd(X)) {
598 					if (!BN_uadd(X, X, n))
599 						goto err;
600 				}
601 				/* now X is even, so we can easily divide it by two */
602 				if (!BN_rshift1(X, X))
603 					goto err;
604 			}
605 			if (shift > 0) {
606 				if (!BN_rshift(B, B, shift))
607 					goto err;
608 			}
609 
610 			/* Same for  A  and  Y.  Afterwards, (2) still holds. */
611 			shift = 0;
612 			while (!BN_is_bit_set(A, shift)) /* note that 0 < A */
613 			{
614 				shift++;
615 
616 				if (BN_is_odd(Y)) {
617 					if (!BN_uadd(Y, Y, n))
618 						goto err;
619 				}
620 				/* now Y is even */
621 				if (!BN_rshift1(Y, Y))
622 					goto err;
623 			}
624 			if (shift > 0) {
625 				if (!BN_rshift(A, A, shift))
626 					goto err;
627 			}
628 
629 			/* We still have (1) and (2).
630 			 * Both  A  and  B  are odd.
631 			 * The following computations ensure that
632 			 *
633 			 *     0 <= B < |n|,
634 			 *      0 < A < |n|,
635 			 * (1) -sign*X*a  ==  B   (mod |n|),
636 			 * (2)  sign*Y*a  ==  A   (mod |n|),
637 			 *
638 			 * and that either  A  or  B  is even in the next iteration.
639 			 */
640 			if (BN_ucmp(B, A) >= 0) {
641 				/* -sign*(X + Y)*a == B - A  (mod |n|) */
642 				if (!BN_uadd(X, X, Y))
643 					goto err;
644 				/* NB: we could use BN_mod_add_quick(X, X, Y, n), but that
645 				 * actually makes the algorithm slower */
646 				if (!BN_usub(B, B, A))
647 					goto err;
648 			} else {
649 				/*  sign*(X + Y)*a == A - B  (mod |n|) */
650 				if (!BN_uadd(Y, Y, X))
651 					goto err;
652 				/* as above, BN_mod_add_quick(Y, Y, X, n) would slow things down */
653 				if (!BN_usub(A, A, B))
654 					goto err;
655 			}
656 		}
657 	} else {
658 		/* general inversion algorithm */
659 
660 		while (!BN_is_zero(B)) {
661 			BIGNUM *tmp;
662 
663 			/*
664 			 *      0 < B < A,
665 			 * (*) -sign*X*a  ==  B   (mod |n|),
666 			 *      sign*Y*a  ==  A   (mod |n|)
667 			 */
668 
669 			/* (D, M) := (A/B, A%B) ... */
670 			if (BN_num_bits(A) == BN_num_bits(B)) {
671 				if (!BN_one(D))
672 					goto err;
673 				if (!BN_sub(M, A, B))
674 					goto err;
675 			} else if (BN_num_bits(A) == BN_num_bits(B) + 1) {
676 				/* A/B is 1, 2, or 3 */
677 				if (!BN_lshift1(T, B))
678 					goto err;
679 				if (BN_ucmp(A, T) < 0) {
680 					/* A < 2*B, so D=1 */
681 					if (!BN_one(D))
682 						goto err;
683 					if (!BN_sub(M, A, B))
684 						goto err;
685 				} else {
686 					/* A >= 2*B, so D=2 or D=3 */
687 					if (!BN_sub(M, A, T))
688 						goto err;
689 					if (!BN_add(D,T,B)) goto err; /* use D (:= 3*B) as temp */
690 						if (BN_ucmp(A, D) < 0) {
691 						/* A < 3*B, so D=2 */
692 						if (!BN_set_word(D, 2))
693 							goto err;
694 						/* M (= A - 2*B) already has the correct value */
695 					} else {
696 						/* only D=3 remains */
697 						if (!BN_set_word(D, 3))
698 							goto err;
699 						/* currently  M = A - 2*B,  but we need  M = A - 3*B */
700 						if (!BN_sub(M, M, B))
701 							goto err;
702 					}
703 				}
704 			} else {
705 				if (!BN_div_nonct(D, M, A, B, ctx))
706 					goto err;
707 			}
708 
709 			/* Now
710 			 *      A = D*B + M;
711 			 * thus we have
712 			 * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
713 			 */
714 			tmp = A; /* keep the BIGNUM object, the value does not matter */
715 
716 			/* (A, B) := (B, A mod B) ... */
717 			A = B;
718 			B = M;
719 			/* ... so we have  0 <= B < A  again */
720 
721 			/* Since the former  M  is now  B  and the former  B  is now  A,
722 			 * (**) translates into
723 			 *       sign*Y*a  ==  D*A + B    (mod |n|),
724 			 * i.e.
725 			 *       sign*Y*a - D*A  ==  B    (mod |n|).
726 			 * Similarly, (*) translates into
727 			 *      -sign*X*a  ==  A          (mod |n|).
728 			 *
729 			 * Thus,
730 			 *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
731 			 * i.e.
732 			 *        sign*(Y + D*X)*a  ==  B  (mod |n|).
733 			 *
734 			 * So if we set  (X, Y, sign) := (Y + D*X, X, -sign),  we arrive back at
735 			 *      -sign*X*a  ==  B   (mod |n|),
736 			 *       sign*Y*a  ==  A   (mod |n|).
737 			 * Note that  X  and  Y  stay non-negative all the time.
738 			 */
739 
740 			/* most of the time D is very small, so we can optimize tmp := D*X+Y */
741 			if (BN_is_one(D)) {
742 				if (!BN_add(tmp, X, Y))
743 					goto err;
744 			} else {
745 				if (BN_is_word(D, 2)) {
746 					if (!BN_lshift1(tmp, X))
747 						goto err;
748 				} else if (BN_is_word(D, 4)) {
749 					if (!BN_lshift(tmp, X, 2))
750 						goto err;
751 				} else if (D->top == 1) {
752 					if (!bn_copy(tmp, X))
753 						goto err;
754 					if (!BN_mul_word(tmp, D->d[0]))
755 						goto err;
756 				} else {
757 					if (!BN_mul(tmp, D,X, ctx))
758 						goto err;
759 				}
760 				if (!BN_add(tmp, tmp, Y))
761 					goto err;
762 			}
763 
764 			M = Y; /* keep the BIGNUM object, the value does not matter */
765 			Y = X;
766 			X = tmp;
767 			sign = -sign;
768 		}
769 	}
770 
771 	/*
772 	 * The while loop (Euclid's algorithm) ends when
773 	 *      A == gcd(a,n);
774 	 * we have
775 	 *       sign*Y*a  ==  A  (mod |n|),
776 	 * where  Y  is non-negative.
777 	 */
778 
779 	if (sign < 0) {
780 		if (!BN_sub(Y, n, Y))
781 			goto err;
782 	}
783 	/* Now  Y*a  ==  A  (mod |n|).  */
784 
785 	if (!BN_is_one(A)) {
786 		BNerror(BN_R_NO_INVERSE);
787 		goto err;
788 	}
789 
790 	if (!BN_nnmod(Y, Y, n, ctx))
791 		goto err;
792 	if (!bn_copy(R, Y))
793 		goto err;
794 
795 	ret = R;
796 
797  err:
798 	if ((ret == NULL) && (in == NULL))
799 		BN_free(R);
800 	BN_CTX_end(ctx);
801 	return (ret);
802 }
803 
804 BIGNUM *
805 BN_mod_inverse(BIGNUM *in, const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
806 {
807 	int ct = ((BN_get_flags(a, BN_FLG_CONSTTIME) != 0) ||
808 	    (BN_get_flags(n, BN_FLG_CONSTTIME) != 0));
809 	return BN_mod_inverse_internal(in, a, n, ctx, ct);
810 }
811 
812 BIGNUM *
813 BN_mod_inverse_nonct(BIGNUM *in, const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
814 {
815 	return BN_mod_inverse_internal(in, a, n, ctx, 0);
816 }
817 
818 BIGNUM *
819 BN_mod_inverse_ct(BIGNUM *in, const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
820 {
821 	return BN_mod_inverse_internal(in, a, n, ctx, 1);
822 }
823