xref: /dragonfly/crypto/libressl/crypto/bn/bn_gcd.c (revision f9993810)
1 /* $OpenBSD: bn_gcd.c,v 1.16 2021/12/26 15:16:50 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_lcl.h"
115 
116 static BIGNUM *euclid(BIGNUM *a, BIGNUM *b);
117 static BIGNUM *BN_gcd_no_branch(BIGNUM *in, const BIGNUM *a, const BIGNUM *n,
118     BN_CTX *ctx);
119 
120 int
121 BN_gcd(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
122 {
123 	BIGNUM *a, *b, *t;
124 	int ret = 0;
125 
126 	bn_check_top(in_a);
127 	bn_check_top(in_b);
128 
129 	BN_CTX_start(ctx);
130 	if ((a = BN_CTX_get(ctx)) == NULL)
131 		goto err;
132 	if ((b = BN_CTX_get(ctx)) == NULL)
133 		goto err;
134 
135 	if (BN_copy(a, in_a) == NULL)
136 		goto err;
137 	if (BN_copy(b, in_b) == NULL)
138 		goto err;
139 	a->neg = 0;
140 	b->neg = 0;
141 
142 	if (BN_cmp(a, b) < 0) {
143 		t = a;
144 		a = b;
145 		b = t;
146 	}
147 	t = euclid(a, b);
148 	if (t == NULL)
149 		goto err;
150 
151 	if (BN_copy(r, t) == NULL)
152 		goto err;
153 	ret = 1;
154 
155 err:
156 	BN_CTX_end(ctx);
157 	bn_check_top(r);
158 	return (ret);
159 }
160 
161 int
162 BN_gcd_ct(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
163 {
164 	if (BN_gcd_no_branch(r, in_a, in_b, ctx) == NULL)
165 		return 0;
166 	return 1;
167 }
168 
169 int
170 BN_gcd_nonct(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
171 {
172 	return BN_gcd(r, in_a, in_b, ctx);
173 }
174 
175 
176 static BIGNUM *
177 euclid(BIGNUM *a, BIGNUM *b)
178 {
179 	BIGNUM *t;
180 	int shifts = 0;
181 
182 	bn_check_top(a);
183 	bn_check_top(b);
184 
185 	/* 0 <= b <= a */
186 	while (!BN_is_zero(b)) {
187 		/* 0 < b <= a */
188 
189 		if (BN_is_odd(a)) {
190 			if (BN_is_odd(b)) {
191 				if (!BN_sub(a, a, b))
192 					goto err;
193 				if (!BN_rshift1(a, a))
194 					goto err;
195 				if (BN_cmp(a, b) < 0) {
196 					t = a;
197 					a = b;
198 					b = t;
199 				}
200 			}
201 			else		/* a odd - b even */
202 			{
203 				if (!BN_rshift1(b, b))
204 					goto err;
205 				if (BN_cmp(a, b) < 0) {
206 					t = a;
207 					a = b;
208 					b = t;
209 				}
210 			}
211 		}
212 		else			/* a is even */
213 		{
214 			if (BN_is_odd(b)) {
215 				if (!BN_rshift1(a, a))
216 					goto err;
217 				if (BN_cmp(a, b) < 0) {
218 					t = a;
219 					a = b;
220 					b = t;
221 				}
222 			}
223 			else		/* a even - b even */
224 			{
225 				if (!BN_rshift1(a, a))
226 					goto err;
227 				if (!BN_rshift1(b, b))
228 					goto err;
229 				shifts++;
230 			}
231 		}
232 		/* 0 <= b <= a */
233 	}
234 
235 	if (shifts) {
236 		if (!BN_lshift(a, a, shifts))
237 			goto err;
238 	}
239 	bn_check_top(a);
240 	return (a);
241 
242 err:
243 	return (NULL);
244 }
245 
246 
247 /* solves ax == 1 (mod n) */
248 static BIGNUM *BN_mod_inverse_no_branch(BIGNUM *in, const BIGNUM *a,
249     const BIGNUM *n, BN_CTX *ctx);
250 
251 static BIGNUM *
252 BN_mod_inverse_internal(BIGNUM *in, const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx,
253     int ct)
254 {
255 	BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
256 	BIGNUM *ret = NULL;
257 	int sign;
258 
259 	if (ct)
260 		return BN_mod_inverse_no_branch(in, a, n, ctx);
261 
262 	bn_check_top(a);
263 	bn_check_top(n);
264 
265 	BN_CTX_start(ctx);
266 	if ((A = BN_CTX_get(ctx)) == NULL)
267 		goto err;
268 	if ((B = BN_CTX_get(ctx)) == NULL)
269 		goto err;
270 	if ((X = BN_CTX_get(ctx)) == NULL)
271 		goto err;
272 	if ((D = BN_CTX_get(ctx)) == NULL)
273 		goto err;
274 	if ((M = BN_CTX_get(ctx)) == NULL)
275 		goto err;
276 	if ((Y = BN_CTX_get(ctx)) == NULL)
277 		goto err;
278 	if ((T = BN_CTX_get(ctx)) == NULL)
279 		goto err;
280 
281 	if (in == NULL)
282 		R = BN_new();
283 	else
284 		R = in;
285 	if (R == NULL)
286 		goto err;
287 
288 	BN_one(X);
289 	BN_zero(Y);
290 	if (BN_copy(B, a) == NULL)
291 		goto err;
292 	if (BN_copy(A, n) == NULL)
293 		goto err;
294 	A->neg = 0;
295 	if (B->neg || (BN_ucmp(B, A) >= 0)) {
296 		if (!BN_nnmod(B, B, A, ctx))
297 			goto err;
298 	}
299 	sign = -1;
300 	/* From  B = a mod |n|,  A = |n|  it follows that
301 	 *
302 	 *      0 <= B < A,
303 	 *     -sign*X*a  ==  B   (mod |n|),
304 	 *      sign*Y*a  ==  A   (mod |n|).
305 	 */
306 
307 	if (BN_is_odd(n) && (BN_num_bits(n) <= (BN_BITS <= 32 ? 450 : 2048))) {
308 		/* Binary inversion algorithm; requires odd modulus.
309 		 * This is faster than the general algorithm if the modulus
310 		 * is sufficiently small (about 400 .. 500 bits on 32-bit
311 		 * sytems, but much more on 64-bit systems) */
312 		int shift;
313 
314 		while (!BN_is_zero(B)) {
315 			/*
316 			 *      0 < B < |n|,
317 			 *      0 < A <= |n|,
318 			 * (1) -sign*X*a  ==  B   (mod |n|),
319 			 * (2)  sign*Y*a  ==  A   (mod |n|)
320 			 */
321 
322 			/* Now divide  B  by the maximum possible power of two in the integers,
323 			 * and divide  X  by the same value mod |n|.
324 			 * When we're done, (1) still holds. */
325 			shift = 0;
326 			while (!BN_is_bit_set(B, shift)) /* note that 0 < B */
327 			{
328 				shift++;
329 
330 				if (BN_is_odd(X)) {
331 					if (!BN_uadd(X, X, n))
332 						goto err;
333 				}
334 				/* now X is even, so we can easily divide it by two */
335 				if (!BN_rshift1(X, X))
336 					goto err;
337 			}
338 			if (shift > 0) {
339 				if (!BN_rshift(B, B, shift))
340 					goto err;
341 			}
342 
343 
344 			/* Same for  A  and  Y.  Afterwards, (2) still holds. */
345 			shift = 0;
346 			while (!BN_is_bit_set(A, shift)) /* note that 0 < A */
347 			{
348 				shift++;
349 
350 				if (BN_is_odd(Y)) {
351 					if (!BN_uadd(Y, Y, n))
352 						goto err;
353 				}
354 				/* now Y is even */
355 				if (!BN_rshift1(Y, Y))
356 					goto err;
357 			}
358 			if (shift > 0) {
359 				if (!BN_rshift(A, A, shift))
360 					goto err;
361 			}
362 
363 
364 			/* We still have (1) and (2).
365 			 * Both  A  and  B  are odd.
366 			 * The following computations ensure that
367 			 *
368 			 *     0 <= B < |n|,
369 			 *      0 < A < |n|,
370 			 * (1) -sign*X*a  ==  B   (mod |n|),
371 			 * (2)  sign*Y*a  ==  A   (mod |n|),
372 			 *
373 			 * and that either  A  or  B  is even in the next iteration.
374 			 */
375 			if (BN_ucmp(B, A) >= 0) {
376 				/* -sign*(X + Y)*a == B - A  (mod |n|) */
377 				if (!BN_uadd(X, X, Y))
378 					goto err;
379 				/* NB: we could use BN_mod_add_quick(X, X, Y, n), but that
380 				 * actually makes the algorithm slower */
381 				if (!BN_usub(B, B, A))
382 					goto err;
383 			} else {
384 				/*  sign*(X + Y)*a == A - B  (mod |n|) */
385 				if (!BN_uadd(Y, Y, X))
386 					goto err;
387 				/* as above, BN_mod_add_quick(Y, Y, X, n) would slow things down */
388 				if (!BN_usub(A, A, B))
389 					goto err;
390 			}
391 		}
392 	} else {
393 		/* general inversion algorithm */
394 
395 		while (!BN_is_zero(B)) {
396 			BIGNUM *tmp;
397 
398 			/*
399 			 *      0 < B < A,
400 			 * (*) -sign*X*a  ==  B   (mod |n|),
401 			 *      sign*Y*a  ==  A   (mod |n|)
402 			 */
403 
404 			/* (D, M) := (A/B, A%B) ... */
405 			if (BN_num_bits(A) == BN_num_bits(B)) {
406 				if (!BN_one(D))
407 					goto err;
408 				if (!BN_sub(M, A, B))
409 					goto err;
410 			} else if (BN_num_bits(A) == BN_num_bits(B) + 1) {
411 				/* A/B is 1, 2, or 3 */
412 				if (!BN_lshift1(T, B))
413 					goto err;
414 				if (BN_ucmp(A, T) < 0) {
415 					/* A < 2*B, so D=1 */
416 					if (!BN_one(D))
417 						goto err;
418 					if (!BN_sub(M, A, B))
419 						goto err;
420 				} else {
421 					/* A >= 2*B, so D=2 or D=3 */
422 					if (!BN_sub(M, A, T))
423 						goto err;
424 					if (!BN_add(D,T,B)) goto err; /* use D (:= 3*B) as temp */
425 						if (BN_ucmp(A, D) < 0) {
426 						/* A < 3*B, so D=2 */
427 						if (!BN_set_word(D, 2))
428 							goto err;
429 						/* M (= A - 2*B) already has the correct value */
430 					} else {
431 						/* only D=3 remains */
432 						if (!BN_set_word(D, 3))
433 							goto err;
434 						/* currently  M = A - 2*B,  but we need  M = A - 3*B */
435 						if (!BN_sub(M, M, B))
436 							goto err;
437 					}
438 				}
439 			} else {
440 				if (!BN_div_nonct(D, M, A, B, ctx))
441 					goto err;
442 			}
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 			/* most of the time D is very small, so we can optimize tmp := D*X+Y */
476 			if (BN_is_one(D)) {
477 				if (!BN_add(tmp, X, Y))
478 					goto err;
479 			} else {
480 				if (BN_is_word(D, 2)) {
481 					if (!BN_lshift1(tmp, X))
482 						goto err;
483 				} else if (BN_is_word(D, 4)) {
484 					if (!BN_lshift(tmp, X, 2))
485 						goto err;
486 				} else if (D->top == 1) {
487 					if (!BN_copy(tmp, X))
488 						goto err;
489 					if (!BN_mul_word(tmp, D->d[0]))
490 						goto err;
491 				} else {
492 					if (!BN_mul(tmp, D,X, ctx))
493 						goto err;
494 				}
495 				if (!BN_add(tmp, tmp, Y))
496 					goto err;
497 			}
498 
499 			M = Y; /* keep the BIGNUM object, the value does not matter */
500 			Y = X;
501 			X = tmp;
502 			sign = -sign;
503 		}
504 	}
505 
506 	/*
507 	 * The while loop (Euclid's algorithm) ends when
508 	 *      A == gcd(a,n);
509 	 * we have
510 	 *       sign*Y*a  ==  A  (mod |n|),
511 	 * where  Y  is non-negative.
512 	 */
513 
514 	if (sign < 0) {
515 		if (!BN_sub(Y, n, Y))
516 			goto err;
517 	}
518 	/* Now  Y*a  ==  A  (mod |n|).  */
519 
520 	if (BN_is_one(A)) {
521 		/* Y*a == 1  (mod |n|) */
522 		if (!Y->neg && BN_ucmp(Y, n) < 0) {
523 			if (!BN_copy(R, Y))
524 				goto err;
525 		} else {
526 			if (!BN_nnmod(R, Y,n, ctx))
527 				goto err;
528 		}
529 	} else {
530 		BNerror(BN_R_NO_INVERSE);
531 		goto err;
532 	}
533 	ret = R;
534 
535 err:
536 	if ((ret == NULL) && (in == NULL))
537 		BN_free(R);
538 	BN_CTX_end(ctx);
539 	bn_check_top(ret);
540 	return (ret);
541 }
542 
543 BIGNUM *
544 BN_mod_inverse(BIGNUM *in, const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
545 {
546 	int ct = ((BN_get_flags(a, BN_FLG_CONSTTIME) != 0) ||
547 	    (BN_get_flags(n, BN_FLG_CONSTTIME) != 0));
548 	return BN_mod_inverse_internal(in, a, n, ctx, ct);
549 }
550 
551 BIGNUM *
552 BN_mod_inverse_nonct(BIGNUM *in, const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
553 {
554 	return BN_mod_inverse_internal(in, a, n, ctx, 0);
555 }
556 
557 BIGNUM *
558 BN_mod_inverse_ct(BIGNUM *in, const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
559 {
560 	return BN_mod_inverse_internal(in, a, n, ctx, 1);
561 }
562 
563 /* BN_mod_inverse_no_branch is a special version of BN_mod_inverse.
564  * It does not contain branches that may leak sensitive information.
565  */
566 static BIGNUM *
567 BN_mod_inverse_no_branch(BIGNUM *in, const BIGNUM *a, const BIGNUM *n,
568     BN_CTX *ctx)
569 {
570 	BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
571 	BIGNUM local_A, local_B;
572 	BIGNUM *pA, *pB;
573 	BIGNUM *ret = NULL;
574 	int sign;
575 
576 	bn_check_top(a);
577 	bn_check_top(n);
578 
579 	BN_init(&local_A);
580 	BN_init(&local_B);
581 
582 	BN_CTX_start(ctx);
583 	if ((A = BN_CTX_get(ctx)) == NULL)
584 		goto err;
585 	if ((B = BN_CTX_get(ctx)) == NULL)
586 		goto err;
587 	if ((X = BN_CTX_get(ctx)) == NULL)
588 		goto err;
589 	if ((D = BN_CTX_get(ctx)) == NULL)
590 		goto err;
591 	if ((M = BN_CTX_get(ctx)) == NULL)
592 		goto err;
593 	if ((Y = BN_CTX_get(ctx)) == NULL)
594 		goto err;
595 	if ((T = BN_CTX_get(ctx)) == NULL)
596 		goto err;
597 
598 	if (in == NULL)
599 		R = BN_new();
600 	else
601 		R = in;
602 	if (R == NULL)
603 		goto err;
604 
605 	BN_one(X);
606 	BN_zero(Y);
607 	if (BN_copy(B, a) == NULL)
608 		goto err;
609 	if (BN_copy(A, n) == NULL)
610 		goto err;
611 	A->neg = 0;
612 
613 	if (B->neg || (BN_ucmp(B, A) >= 0)) {
614 		/*
615 		 * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
616 		 * BN_div_no_branch will be called eventually.
617 		 */
618 		pB = &local_B;
619 		/* BN_init() done at the top of the function. */
620 		BN_with_flags(pB, B, BN_FLG_CONSTTIME);
621 		if (!BN_nnmod(B, pB, A, ctx))
622 			goto err;
623 	}
624 	sign = -1;
625 	/* From  B = a mod |n|,  A = |n|  it follows that
626 	 *
627 	 *      0 <= B < A,
628 	 *     -sign*X*a  ==  B   (mod |n|),
629 	 *      sign*Y*a  ==  A   (mod |n|).
630 	 */
631 
632 	while (!BN_is_zero(B)) {
633 		BIGNUM *tmp;
634 
635 		/*
636 		 *      0 < B < A,
637 		 * (*) -sign*X*a  ==  B   (mod |n|),
638 		 *      sign*Y*a  ==  A   (mod |n|)
639 		 */
640 
641 		/*
642 		 * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
643 		 * BN_div_no_branch will be called eventually.
644 		 */
645 		pA = &local_A;
646 		/* BN_init() done at the top of the function. */
647 		BN_with_flags(pA, A, BN_FLG_CONSTTIME);
648 
649 		/* (D, M) := (A/B, A%B) ... */
650 		if (!BN_div_ct(D, M, pA, B, ctx))
651 			goto err;
652 
653 		/* Now
654 		 *      A = D*B + M;
655 		 * thus we have
656 		 * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
657 		 */
658 		tmp = A; /* keep the BIGNUM object, the value does not matter */
659 
660 		/* (A, B) := (B, A mod B) ... */
661 		A = B;
662 		B = M;
663 		/* ... so we have  0 <= B < A  again */
664 
665 		/* Since the former  M  is now  B  and the former  B  is now  A,
666 		 * (**) translates into
667 		 *       sign*Y*a  ==  D*A + B    (mod |n|),
668 		 * i.e.
669 		 *       sign*Y*a - D*A  ==  B    (mod |n|).
670 		 * Similarly, (*) translates into
671 		 *      -sign*X*a  ==  A          (mod |n|).
672 		 *
673 		 * Thus,
674 		 *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
675 		 * i.e.
676 		 *        sign*(Y + D*X)*a  ==  B  (mod |n|).
677 		 *
678 		 * So if we set  (X, Y, sign) := (Y + D*X, X, -sign),  we arrive back at
679 		 *      -sign*X*a  ==  B   (mod |n|),
680 		 *       sign*Y*a  ==  A   (mod |n|).
681 		 * Note that  X  and  Y  stay non-negative all the time.
682 		 */
683 
684 		if (!BN_mul(tmp, D, X, ctx))
685 			goto err;
686 		if (!BN_add(tmp, tmp, Y))
687 			goto err;
688 
689 		M = Y; /* keep the BIGNUM object, the value does not matter */
690 		Y = X;
691 		X = tmp;
692 		sign = -sign;
693 	}
694 
695 	/*
696 	 * The while loop (Euclid's algorithm) ends when
697 	 *      A == gcd(a,n);
698 	 * we have
699 	 *       sign*Y*a  ==  A  (mod |n|),
700 	 * where  Y  is non-negative.
701 	 */
702 
703 	if (sign < 0) {
704 		if (!BN_sub(Y, n, Y))
705 			goto err;
706 	}
707 	/* Now  Y*a  ==  A  (mod |n|).  */
708 
709 	if (BN_is_one(A)) {
710 		/* Y*a == 1  (mod |n|) */
711 		if (!Y->neg && BN_ucmp(Y, n) < 0) {
712 			if (!BN_copy(R, Y))
713 				goto err;
714 		} else {
715 			if (!BN_nnmod(R, Y, n, ctx))
716 				goto err;
717 		}
718 	} else {
719 		BNerror(BN_R_NO_INVERSE);
720 		goto err;
721 	}
722 	ret = R;
723 
724 err:
725 	if ((ret == NULL) && (in == NULL))
726 		BN_free(R);
727 	BN_CTX_end(ctx);
728 	bn_check_top(ret);
729 	return (ret);
730 }
731 
732 /*
733  * BN_gcd_no_branch is a special version of BN_mod_inverse_no_branch.
734  * that returns the GCD.
735  */
736 static BIGNUM *
737 BN_gcd_no_branch(BIGNUM *in, const BIGNUM *a, const BIGNUM *n,
738     BN_CTX *ctx)
739 {
740 	BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
741 	BIGNUM local_A, local_B;
742 	BIGNUM *pA, *pB;
743 	BIGNUM *ret = NULL;
744 	int sign;
745 
746 	if (in == NULL)
747 		goto err;
748 	R = in;
749 
750 	BN_init(&local_A);
751 	BN_init(&local_B);
752 
753 	bn_check_top(a);
754 	bn_check_top(n);
755 
756 	BN_CTX_start(ctx);
757 	if ((A = BN_CTX_get(ctx)) == NULL)
758 		goto err;
759 	if ((B = BN_CTX_get(ctx)) == NULL)
760 		goto err;
761 	if ((X = BN_CTX_get(ctx)) == NULL)
762 		goto err;
763 	if ((D = BN_CTX_get(ctx)) == NULL)
764 		goto err;
765 	if ((M = BN_CTX_get(ctx)) == NULL)
766 		goto err;
767 	if ((Y = BN_CTX_get(ctx)) == NULL)
768 		goto err;
769 	if ((T = BN_CTX_get(ctx)) == NULL)
770 		goto err;
771 
772 	BN_one(X);
773 	BN_zero(Y);
774 	if (BN_copy(B, a) == NULL)
775 		goto err;
776 	if (BN_copy(A, n) == NULL)
777 		goto err;
778 	A->neg = 0;
779 
780 	if (B->neg || (BN_ucmp(B, A) >= 0)) {
781 		/*
782 		 * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
783 		 * BN_div_no_branch will be called eventually.
784 		 */
785 		pB = &local_B;
786 		/* BN_init() done at the top of the function. */
787 		BN_with_flags(pB, B, BN_FLG_CONSTTIME);
788 		if (!BN_nnmod(B, pB, A, ctx))
789 			goto err;
790 	}
791 	sign = -1;
792 	/* From  B = a mod |n|,  A = |n|  it follows that
793 	 *
794 	 *      0 <= B < A,
795 	 *     -sign*X*a  ==  B   (mod |n|),
796 	 *      sign*Y*a  ==  A   (mod |n|).
797 	 */
798 
799 	while (!BN_is_zero(B)) {
800 		BIGNUM *tmp;
801 
802 		/*
803 		 *      0 < B < A,
804 		 * (*) -sign*X*a  ==  B   (mod |n|),
805 		 *      sign*Y*a  ==  A   (mod |n|)
806 		 */
807 
808 		/*
809 		 * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
810 		 * BN_div_no_branch will be called eventually.
811 		 */
812 		pA = &local_A;
813 		/* BN_init() done at the top of the function. */
814 		BN_with_flags(pA, A, BN_FLG_CONSTTIME);
815 
816 		/* (D, M) := (A/B, A%B) ... */
817 		if (!BN_div_ct(D, M, pA, B, ctx))
818 			goto err;
819 
820 		/* Now
821 		 *      A = D*B + M;
822 		 * thus we have
823 		 * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
824 		 */
825 		tmp = A; /* keep the BIGNUM object, the value does not matter */
826 
827 		/* (A, B) := (B, A mod B) ... */
828 		A = B;
829 		B = M;
830 		/* ... so we have  0 <= B < A  again */
831 
832 		/* Since the former  M  is now  B  and the former  B  is now  A,
833 		 * (**) translates into
834 		 *       sign*Y*a  ==  D*A + B    (mod |n|),
835 		 * i.e.
836 		 *       sign*Y*a - D*A  ==  B    (mod |n|).
837 		 * Similarly, (*) translates into
838 		 *      -sign*X*a  ==  A          (mod |n|).
839 		 *
840 		 * Thus,
841 		 *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
842 		 * i.e.
843 		 *        sign*(Y + D*X)*a  ==  B  (mod |n|).
844 		 *
845 		 * So if we set  (X, Y, sign) := (Y + D*X, X, -sign),  we arrive back at
846 		 *      -sign*X*a  ==  B   (mod |n|),
847 		 *       sign*Y*a  ==  A   (mod |n|).
848 		 * Note that  X  and  Y  stay non-negative all the time.
849 		 */
850 
851 		if (!BN_mul(tmp, D, X, ctx))
852 			goto err;
853 		if (!BN_add(tmp, tmp, Y))
854 			goto err;
855 
856 		M = Y; /* keep the BIGNUM object, the value does not matter */
857 		Y = X;
858 		X = tmp;
859 		sign = -sign;
860 	}
861 
862 	/*
863 	 * The while loop (Euclid's algorithm) ends when
864 	 *      A == gcd(a,n);
865 	 */
866 
867 	if (!BN_copy(R, A))
868 		goto err;
869 	ret = R;
870 err:
871 	if ((ret == NULL) && (in == NULL))
872 		BN_free(R);
873 	BN_CTX_end(ctx);
874 	bn_check_top(ret);
875 	return (ret);
876 }
877