RT4593: Add space after comma (doc nits)
[openssl.git] / crypto / bn / bn_gcd.c
1 /*
2  * Copyright 1995-2016 The OpenSSL Project Authors. All Rights Reserved.
3  *
4  * Licensed under the OpenSSL license (the "License").  You may not use
5  * this file except in compliance with the License.  You can obtain a copy
6  * in the file LICENSE in the source distribution or at
7  * https://www.openssl.org/source/license.html
8  */
9
10 #include "internal/cryptlib.h"
11 #include "bn_lcl.h"
12
13 static BIGNUM *euclid(BIGNUM *a, BIGNUM *b);
14
15 int BN_gcd(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
16 {
17     BIGNUM *a, *b, *t;
18     int ret = 0;
19
20     bn_check_top(in_a);
21     bn_check_top(in_b);
22
23     BN_CTX_start(ctx);
24     a = BN_CTX_get(ctx);
25     b = BN_CTX_get(ctx);
26     if (a == NULL || b == NULL)
27         goto err;
28
29     if (BN_copy(a, in_a) == NULL)
30         goto err;
31     if (BN_copy(b, in_b) == NULL)
32         goto err;
33     a->neg = 0;
34     b->neg = 0;
35
36     if (BN_cmp(a, b) < 0) {
37         t = a;
38         a = b;
39         b = t;
40     }
41     t = euclid(a, b);
42     if (t == NULL)
43         goto err;
44
45     if (BN_copy(r, t) == NULL)
46         goto err;
47     ret = 1;
48  err:
49     BN_CTX_end(ctx);
50     bn_check_top(r);
51     return (ret);
52 }
53
54 static BIGNUM *euclid(BIGNUM *a, BIGNUM *b)
55 {
56     BIGNUM *t;
57     int shifts = 0;
58
59     bn_check_top(a);
60     bn_check_top(b);
61
62     /* 0 <= b <= a */
63     while (!BN_is_zero(b)) {
64         /* 0 < b <= a */
65
66         if (BN_is_odd(a)) {
67             if (BN_is_odd(b)) {
68                 if (!BN_sub(a, a, b))
69                     goto err;
70                 if (!BN_rshift1(a, a))
71                     goto err;
72                 if (BN_cmp(a, b) < 0) {
73                     t = a;
74                     a = b;
75                     b = t;
76                 }
77             } else {            /* a odd - b even */
78
79                 if (!BN_rshift1(b, b))
80                     goto err;
81                 if (BN_cmp(a, b) < 0) {
82                     t = a;
83                     a = b;
84                     b = t;
85                 }
86             }
87         } else {                /* a is even */
88
89             if (BN_is_odd(b)) {
90                 if (!BN_rshift1(a, a))
91                     goto err;
92                 if (BN_cmp(a, b) < 0) {
93                     t = a;
94                     a = b;
95                     b = t;
96                 }
97             } else {            /* a even - b even */
98
99                 if (!BN_rshift1(a, a))
100                     goto err;
101                 if (!BN_rshift1(b, b))
102                     goto err;
103                 shifts++;
104             }
105         }
106         /* 0 <= b <= a */
107     }
108
109     if (shifts) {
110         if (!BN_lshift(a, a, shifts))
111             goto err;
112     }
113     bn_check_top(a);
114     return (a);
115  err:
116     return (NULL);
117 }
118
119 /* solves ax == 1 (mod n) */
120 static BIGNUM *BN_mod_inverse_no_branch(BIGNUM *in,
121                                         const BIGNUM *a, const BIGNUM *n,
122                                         BN_CTX *ctx);
123
124 BIGNUM *BN_mod_inverse(BIGNUM *in,
125                        const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
126 {
127     BIGNUM *rv;
128     int noinv;
129     rv = int_bn_mod_inverse(in, a, n, ctx, &noinv);
130     if (noinv)
131         BNerr(BN_F_BN_MOD_INVERSE, BN_R_NO_INVERSE);
132     return rv;
133 }
134
135 BIGNUM *int_bn_mod_inverse(BIGNUM *in,
136                            const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx,
137                            int *pnoinv)
138 {
139     BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
140     BIGNUM *ret = NULL;
141     int sign;
142
143     if (pnoinv)
144         *pnoinv = 0;
145
146     if ((BN_get_flags(a, BN_FLG_CONSTTIME) != 0)
147         || (BN_get_flags(n, BN_FLG_CONSTTIME) != 0)) {
148         return BN_mod_inverse_no_branch(in, a, n, ctx);
149     }
150
151     bn_check_top(a);
152     bn_check_top(n);
153
154     BN_CTX_start(ctx);
155     A = BN_CTX_get(ctx);
156     B = BN_CTX_get(ctx);
157     X = BN_CTX_get(ctx);
158     D = BN_CTX_get(ctx);
159     M = BN_CTX_get(ctx);
160     Y = BN_CTX_get(ctx);
161     T = BN_CTX_get(ctx);
162     if (T == NULL)
163         goto err;
164
165     if (in == NULL)
166         R = BN_new();
167     else
168         R = in;
169     if (R == NULL)
170         goto err;
171
172     BN_one(X);
173     BN_zero(Y);
174     if (BN_copy(B, a) == NULL)
175         goto err;
176     if (BN_copy(A, n) == NULL)
177         goto err;
178     A->neg = 0;
179     if (B->neg || (BN_ucmp(B, A) >= 0)) {
180         if (!BN_nnmod(B, B, A, ctx))
181             goto err;
182     }
183     sign = -1;
184     /*-
185      * From  B = a mod |n|,  A = |n|  it follows that
186      *
187      *      0 <= B < A,
188      *     -sign*X*a  ==  B   (mod |n|),
189      *      sign*Y*a  ==  A   (mod |n|).
190      */
191
192     if (BN_is_odd(n) && (BN_num_bits(n) <= 2048)) {
193         /*
194          * Binary inversion algorithm; requires odd modulus. This is faster
195          * than the general algorithm if the modulus is sufficiently small
196          * (about 400 .. 500 bits on 32-bit systems, but much more on 64-bit
197          * systems)
198          */
199         int shift;
200
201         while (!BN_is_zero(B)) {
202             /*-
203              *      0 < B < |n|,
204              *      0 < A <= |n|,
205              * (1) -sign*X*a  ==  B   (mod |n|),
206              * (2)  sign*Y*a  ==  A   (mod |n|)
207              */
208
209             /*
210              * Now divide B by the maximum possible power of two in the
211              * integers, and divide X by the same value mod |n|. When we're
212              * done, (1) still holds.
213              */
214             shift = 0;
215             while (!BN_is_bit_set(B, shift)) { /* note that 0 < B */
216                 shift++;
217
218                 if (BN_is_odd(X)) {
219                     if (!BN_uadd(X, X, n))
220                         goto err;
221                 }
222                 /*
223                  * now X is even, so we can easily divide it by two
224                  */
225                 if (!BN_rshift1(X, X))
226                     goto err;
227             }
228             if (shift > 0) {
229                 if (!BN_rshift(B, B, shift))
230                     goto err;
231             }
232
233             /*
234              * Same for A and Y.  Afterwards, (2) still holds.
235              */
236             shift = 0;
237             while (!BN_is_bit_set(A, shift)) { /* note that 0 < A */
238                 shift++;
239
240                 if (BN_is_odd(Y)) {
241                     if (!BN_uadd(Y, Y, n))
242                         goto err;
243                 }
244                 /* now Y is even */
245                 if (!BN_rshift1(Y, Y))
246                     goto err;
247             }
248             if (shift > 0) {
249                 if (!BN_rshift(A, A, shift))
250                     goto err;
251             }
252
253             /*-
254              * We still have (1) and (2).
255              * Both  A  and  B  are odd.
256              * The following computations ensure that
257              *
258              *     0 <= B < |n|,
259              *      0 < A < |n|,
260              * (1) -sign*X*a  ==  B   (mod |n|),
261              * (2)  sign*Y*a  ==  A   (mod |n|),
262              *
263              * and that either  A  or  B  is even in the next iteration.
264              */
265             if (BN_ucmp(B, A) >= 0) {
266                 /* -sign*(X + Y)*a == B - A  (mod |n|) */
267                 if (!BN_uadd(X, X, Y))
268                     goto err;
269                 /*
270                  * NB: we could use BN_mod_add_quick(X, X, Y, n), but that
271                  * actually makes the algorithm slower
272                  */
273                 if (!BN_usub(B, B, A))
274                     goto err;
275             } else {
276                 /*  sign*(X + Y)*a == A - B  (mod |n|) */
277                 if (!BN_uadd(Y, Y, X))
278                     goto err;
279                 /*
280                  * as above, BN_mod_add_quick(Y, Y, X, n) would slow things
281                  * down
282                  */
283                 if (!BN_usub(A, A, B))
284                     goto err;
285             }
286         }
287     } else {
288         /* general inversion algorithm */
289
290         while (!BN_is_zero(B)) {
291             BIGNUM *tmp;
292
293             /*-
294              *      0 < B < A,
295              * (*) -sign*X*a  ==  B   (mod |n|),
296              *      sign*Y*a  ==  A   (mod |n|)
297              */
298
299             /* (D, M) := (A/B, A%B) ... */
300             if (BN_num_bits(A) == BN_num_bits(B)) {
301                 if (!BN_one(D))
302                     goto err;
303                 if (!BN_sub(M, A, B))
304                     goto err;
305             } else if (BN_num_bits(A) == BN_num_bits(B) + 1) {
306                 /* A/B is 1, 2, or 3 */
307                 if (!BN_lshift1(T, B))
308                     goto err;
309                 if (BN_ucmp(A, T) < 0) {
310                     /* A < 2*B, so D=1 */
311                     if (!BN_one(D))
312                         goto err;
313                     if (!BN_sub(M, A, B))
314                         goto err;
315                 } else {
316                     /* A >= 2*B, so D=2 or D=3 */
317                     if (!BN_sub(M, A, T))
318                         goto err;
319                     if (!BN_add(D, T, B))
320                         goto err; /* use D (:= 3*B) as temp */
321                     if (BN_ucmp(A, D) < 0) {
322                         /* A < 3*B, so D=2 */
323                         if (!BN_set_word(D, 2))
324                             goto err;
325                         /*
326                          * M (= A - 2*B) already has the correct value
327                          */
328                     } else {
329                         /* only D=3 remains */
330                         if (!BN_set_word(D, 3))
331                             goto err;
332                         /*
333                          * currently M = A - 2*B, but we need M = A - 3*B
334                          */
335                         if (!BN_sub(M, M, B))
336                             goto err;
337                     }
338                 }
339             } else {
340                 if (!BN_div(D, M, A, B, ctx))
341                     goto err;
342             }
343
344             /*-
345              * Now
346              *      A = D*B + M;
347              * thus we have
348              * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
349              */
350
351             tmp = A;            /* keep the BIGNUM object, the value does not
352                                  * matter */
353
354             /* (A, B) := (B, A mod B) ... */
355             A = B;
356             B = M;
357             /* ... so we have  0 <= B < A  again */
358
359             /*-
360              * Since the former  M  is now  B  and the former  B  is now  A,
361              * (**) translates into
362              *       sign*Y*a  ==  D*A + B    (mod |n|),
363              * i.e.
364              *       sign*Y*a - D*A  ==  B    (mod |n|).
365              * Similarly, (*) translates into
366              *      -sign*X*a  ==  A          (mod |n|).
367              *
368              * Thus,
369              *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
370              * i.e.
371              *        sign*(Y + D*X)*a  ==  B  (mod |n|).
372              *
373              * So if we set  (X, Y, sign) := (Y + D*X, X, -sign), we arrive back at
374              *      -sign*X*a  ==  B   (mod |n|),
375              *       sign*Y*a  ==  A   (mod |n|).
376              * Note that  X  and  Y  stay non-negative all the time.
377              */
378
379             /*
380              * most of the time D is very small, so we can optimize tmp :=
381              * D*X+Y
382              */
383             if (BN_is_one(D)) {
384                 if (!BN_add(tmp, X, Y))
385                     goto err;
386             } else {
387                 if (BN_is_word(D, 2)) {
388                     if (!BN_lshift1(tmp, X))
389                         goto err;
390                 } else if (BN_is_word(D, 4)) {
391                     if (!BN_lshift(tmp, X, 2))
392                         goto err;
393                 } else if (D->top == 1) {
394                     if (!BN_copy(tmp, X))
395                         goto err;
396                     if (!BN_mul_word(tmp, D->d[0]))
397                         goto err;
398                 } else {
399                     if (!BN_mul(tmp, D, X, ctx))
400                         goto err;
401                 }
402                 if (!BN_add(tmp, tmp, Y))
403                     goto err;
404             }
405
406             M = Y;              /* keep the BIGNUM object, the value does not
407                                  * matter */
408             Y = X;
409             X = tmp;
410             sign = -sign;
411         }
412     }
413
414     /*-
415      * The while loop (Euclid's algorithm) ends when
416      *      A == gcd(a,n);
417      * we have
418      *       sign*Y*a  ==  A  (mod |n|),
419      * where  Y  is non-negative.
420      */
421
422     if (sign < 0) {
423         if (!BN_sub(Y, n, Y))
424             goto err;
425     }
426     /* Now  Y*a  ==  A  (mod |n|).  */
427
428     if (BN_is_one(A)) {
429         /* Y*a == 1  (mod |n|) */
430         if (!Y->neg && BN_ucmp(Y, n) < 0) {
431             if (!BN_copy(R, Y))
432                 goto err;
433         } else {
434             if (!BN_nnmod(R, Y, n, ctx))
435                 goto err;
436         }
437     } else {
438         if (pnoinv)
439             *pnoinv = 1;
440         goto err;
441     }
442     ret = R;
443  err:
444     if ((ret == NULL) && (in == NULL))
445         BN_free(R);
446     BN_CTX_end(ctx);
447     bn_check_top(ret);
448     return (ret);
449 }
450
451 /*
452  * BN_mod_inverse_no_branch is a special version of BN_mod_inverse. It does
453  * not contain branches that may leak sensitive information.
454  */
455 static BIGNUM *BN_mod_inverse_no_branch(BIGNUM *in,
456                                         const BIGNUM *a, const BIGNUM *n,
457                                         BN_CTX *ctx)
458 {
459     BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
460     BIGNUM *ret = NULL;
461     int sign;
462
463     bn_check_top(a);
464     bn_check_top(n);
465
466     BN_CTX_start(ctx);
467     A = BN_CTX_get(ctx);
468     B = BN_CTX_get(ctx);
469     X = BN_CTX_get(ctx);
470     D = BN_CTX_get(ctx);
471     M = BN_CTX_get(ctx);
472     Y = BN_CTX_get(ctx);
473     T = BN_CTX_get(ctx);
474     if (T == NULL)
475         goto err;
476
477     if (in == NULL)
478         R = BN_new();
479     else
480         R = in;
481     if (R == NULL)
482         goto err;
483
484     BN_one(X);
485     BN_zero(Y);
486     if (BN_copy(B, a) == NULL)
487         goto err;
488     if (BN_copy(A, n) == NULL)
489         goto err;
490     A->neg = 0;
491
492     if (B->neg || (BN_ucmp(B, A) >= 0)) {
493         /*
494          * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
495          * BN_div_no_branch will be called eventually.
496          */
497          {
498             BIGNUM local_B;
499             bn_init(&local_B);
500             BN_with_flags(&local_B, B, BN_FLG_CONSTTIME);
501             if (!BN_nnmod(B, &local_B, A, ctx))
502                 goto err;
503             /* Ensure local_B goes out of scope before any further use of B */
504         }
505     }
506     sign = -1;
507     /*-
508      * From  B = a mod |n|,  A = |n|  it follows that
509      *
510      *      0 <= B < A,
511      *     -sign*X*a  ==  B   (mod |n|),
512      *      sign*Y*a  ==  A   (mod |n|).
513      */
514
515     while (!BN_is_zero(B)) {
516         BIGNUM *tmp;
517
518         /*-
519          *      0 < B < A,
520          * (*) -sign*X*a  ==  B   (mod |n|),
521          *      sign*Y*a  ==  A   (mod |n|)
522          */
523
524         /*
525          * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
526          * BN_div_no_branch will be called eventually.
527          */
528         {
529             BIGNUM local_A;
530             bn_init(&local_A);
531             BN_with_flags(&local_A, A, BN_FLG_CONSTTIME);
532
533             /* (D, M) := (A/B, A%B) ... */
534             if (!BN_div(D, M, &local_A, B, ctx))
535                 goto err;
536             /* Ensure local_A goes out of scope before any further use of A */
537         }
538
539         /*-
540          * Now
541          *      A = D*B + M;
542          * thus we have
543          * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
544          */
545
546         tmp = A;                /* keep the BIGNUM object, the value does not
547                                  * matter */
548
549         /* (A, B) := (B, A mod B) ... */
550         A = B;
551         B = M;
552         /* ... so we have  0 <= B < A  again */
553
554         /*-
555          * Since the former  M  is now  B  and the former  B  is now  A,
556          * (**) translates into
557          *       sign*Y*a  ==  D*A + B    (mod |n|),
558          * i.e.
559          *       sign*Y*a - D*A  ==  B    (mod |n|).
560          * Similarly, (*) translates into
561          *      -sign*X*a  ==  A          (mod |n|).
562          *
563          * Thus,
564          *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
565          * i.e.
566          *        sign*(Y + D*X)*a  ==  B  (mod |n|).
567          *
568          * So if we set  (X, Y, sign) := (Y + D*X, X, -sign), we arrive back at
569          *      -sign*X*a  ==  B   (mod |n|),
570          *       sign*Y*a  ==  A   (mod |n|).
571          * Note that  X  and  Y  stay non-negative all the time.
572          */
573
574         if (!BN_mul(tmp, D, X, ctx))
575             goto err;
576         if (!BN_add(tmp, tmp, Y))
577             goto err;
578
579         M = Y;                  /* keep the BIGNUM object, the value does not
580                                  * matter */
581         Y = X;
582         X = tmp;
583         sign = -sign;
584     }
585
586     /*-
587      * The while loop (Euclid's algorithm) ends when
588      *      A == gcd(a,n);
589      * we have
590      *       sign*Y*a  ==  A  (mod |n|),
591      * where  Y  is non-negative.
592      */
593
594     if (sign < 0) {
595         if (!BN_sub(Y, n, Y))
596             goto err;
597     }
598     /* Now  Y*a  ==  A  (mod |n|).  */
599
600     if (BN_is_one(A)) {
601         /* Y*a == 1  (mod |n|) */
602         if (!Y->neg && BN_ucmp(Y, n) < 0) {
603             if (!BN_copy(R, Y))
604                 goto err;
605         } else {
606             if (!BN_nnmod(R, Y, n, ctx))
607                 goto err;
608         }
609     } else {
610         BNerr(BN_F_BN_MOD_INVERSE_NO_BRANCH, BN_R_NO_INVERSE);
611         goto err;
612     }
613     ret = R;
614  err:
615     if ((ret == NULL) && (in == NULL))
616         BN_free(R);
617     BN_CTX_end(ctx);
618     bn_check_top(ret);
619     return (ret);
620 }