5eda65cfbb048105a697c5eaca64387454642fd5
[openssl.git] / crypto / bn / bn_mul.c
1 /*
2  * Copyright 1995-2018 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 <assert.h>
11 #include "internal/cryptlib.h"
12 #include "bn_lcl.h"
13
14 #if defined(OPENSSL_NO_ASM) || !defined(OPENSSL_BN_ASM_PART_WORDS)
15 /*
16  * Here follows specialised variants of bn_add_words() and bn_sub_words().
17  * They have the property performing operations on arrays of different sizes.
18  * The sizes of those arrays is expressed through cl, which is the common
19  * length ( basically, min(len(a),len(b)) ), and dl, which is the delta
20  * between the two lengths, calculated as len(a)-len(b). All lengths are the
21  * number of BN_ULONGs...  For the operations that require a result array as
22  * parameter, it must have the length cl+abs(dl). These functions should
23  * probably end up in bn_asm.c as soon as there are assembler counterparts
24  * for the systems that use assembler files.
25  */
26
27 BN_ULONG bn_sub_part_words(BN_ULONG *r,
28                            const BN_ULONG *a, const BN_ULONG *b,
29                            int cl, int dl)
30 {
31     BN_ULONG c, t;
32
33     assert(cl >= 0);
34     c = bn_sub_words(r, a, b, cl);
35
36     if (dl == 0)
37         return c;
38
39     r += cl;
40     a += cl;
41     b += cl;
42
43     if (dl < 0) {
44         for (;;) {
45             t = b[0];
46             r[0] = (0 - t - c) & BN_MASK2;
47             if (t != 0)
48                 c = 1;
49             if (++dl >= 0)
50                 break;
51
52             t = b[1];
53             r[1] = (0 - t - c) & BN_MASK2;
54             if (t != 0)
55                 c = 1;
56             if (++dl >= 0)
57                 break;
58
59             t = b[2];
60             r[2] = (0 - t - c) & BN_MASK2;
61             if (t != 0)
62                 c = 1;
63             if (++dl >= 0)
64                 break;
65
66             t = b[3];
67             r[3] = (0 - t - c) & BN_MASK2;
68             if (t != 0)
69                 c = 1;
70             if (++dl >= 0)
71                 break;
72
73             b += 4;
74             r += 4;
75         }
76     } else {
77         int save_dl = dl;
78         while (c) {
79             t = a[0];
80             r[0] = (t - c) & BN_MASK2;
81             if (t != 0)
82                 c = 0;
83             if (--dl <= 0)
84                 break;
85
86             t = a[1];
87             r[1] = (t - c) & BN_MASK2;
88             if (t != 0)
89                 c = 0;
90             if (--dl <= 0)
91                 break;
92
93             t = a[2];
94             r[2] = (t - c) & BN_MASK2;
95             if (t != 0)
96                 c = 0;
97             if (--dl <= 0)
98                 break;
99
100             t = a[3];
101             r[3] = (t - c) & BN_MASK2;
102             if (t != 0)
103                 c = 0;
104             if (--dl <= 0)
105                 break;
106
107             save_dl = dl;
108             a += 4;
109             r += 4;
110         }
111         if (dl > 0) {
112             if (save_dl > dl) {
113                 switch (save_dl - dl) {
114                 case 1:
115                     r[1] = a[1];
116                     if (--dl <= 0)
117                         break;
118                     /* fall thru */
119                 case 2:
120                     r[2] = a[2];
121                     if (--dl <= 0)
122                         break;
123                     /* fall thru */
124                 case 3:
125                     r[3] = a[3];
126                     if (--dl <= 0)
127                         break;
128                 }
129                 a += 4;
130                 r += 4;
131             }
132         }
133         if (dl > 0) {
134             for (;;) {
135                 r[0] = a[0];
136                 if (--dl <= 0)
137                     break;
138                 r[1] = a[1];
139                 if (--dl <= 0)
140                     break;
141                 r[2] = a[2];
142                 if (--dl <= 0)
143                     break;
144                 r[3] = a[3];
145                 if (--dl <= 0)
146                     break;
147
148                 a += 4;
149                 r += 4;
150             }
151         }
152     }
153     return c;
154 }
155 #endif
156
157 #ifdef BN_RECURSION
158 /*
159  * Karatsuba recursive multiplication algorithm (cf. Knuth, The Art of
160  * Computer Programming, Vol. 2)
161  */
162
163 /*-
164  * r is 2*n2 words in size,
165  * a and b are both n2 words in size.
166  * n2 must be a power of 2.
167  * We multiply and return the result.
168  * t must be 2*n2 words in size
169  * We calculate
170  * a[0]*b[0]
171  * a[0]*b[0]+a[1]*b[1]+(a[0]-a[1])*(b[1]-b[0])
172  * a[1]*b[1]
173  */
174 /* dnX may not be positive, but n2/2+dnX has to be */
175 void bn_mul_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n2,
176                       int dna, int dnb, BN_ULONG *t)
177 {
178     int n = n2 / 2, c1, c2;
179     int tna = n + dna, tnb = n + dnb;
180     unsigned int neg, zero;
181     BN_ULONG ln, lo, *p;
182
183 # ifdef BN_MUL_COMBA
184 #  if 0
185     if (n2 == 4) {
186         bn_mul_comba4(r, a, b);
187         return;
188     }
189 #  endif
190     /*
191      * Only call bn_mul_comba 8 if n2 == 8 and the two arrays are complete
192      * [steve]
193      */
194     if (n2 == 8 && dna == 0 && dnb == 0) {
195         bn_mul_comba8(r, a, b);
196         return;
197     }
198 # endif                         /* BN_MUL_COMBA */
199     /* Else do normal multiply */
200     if (n2 < BN_MUL_RECURSIVE_SIZE_NORMAL) {
201         bn_mul_normal(r, a, n2 + dna, b, n2 + dnb);
202         if ((dna + dnb) < 0)
203             memset(&r[2 * n2 + dna + dnb], 0,
204                    sizeof(BN_ULONG) * -(dna + dnb));
205         return;
206     }
207     /* r=(a[0]-a[1])*(b[1]-b[0]) */
208     c1 = bn_cmp_part_words(a, &(a[n]), tna, n - tna);
209     c2 = bn_cmp_part_words(&(b[n]), b, tnb, tnb - n);
210     zero = neg = 0;
211     switch (c1 * 3 + c2) {
212     case -4:
213         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
214         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
215         break;
216     case -3:
217         zero = 1;
218         break;
219     case -2:
220         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
221         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n); /* + */
222         neg = 1;
223         break;
224     case -1:
225     case 0:
226     case 1:
227         zero = 1;
228         break;
229     case 2:
230         bn_sub_part_words(t, a, &(a[n]), tna, n - tna); /* + */
231         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
232         neg = 1;
233         break;
234     case 3:
235         zero = 1;
236         break;
237     case 4:
238         bn_sub_part_words(t, a, &(a[n]), tna, n - tna);
239         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n);
240         break;
241     }
242
243 # ifdef BN_MUL_COMBA
244     if (n == 4 && dna == 0 && dnb == 0) { /* XXX: bn_mul_comba4 could take
245                                            * extra args to do this well */
246         if (!zero)
247             bn_mul_comba4(&(t[n2]), t, &(t[n]));
248         else
249             memset(&t[n2], 0, sizeof(*t) * 8);
250
251         bn_mul_comba4(r, a, b);
252         bn_mul_comba4(&(r[n2]), &(a[n]), &(b[n]));
253     } else if (n == 8 && dna == 0 && dnb == 0) { /* XXX: bn_mul_comba8 could
254                                                   * take extra args to do
255                                                   * this well */
256         if (!zero)
257             bn_mul_comba8(&(t[n2]), t, &(t[n]));
258         else
259             memset(&t[n2], 0, sizeof(*t) * 16);
260
261         bn_mul_comba8(r, a, b);
262         bn_mul_comba8(&(r[n2]), &(a[n]), &(b[n]));
263     } else
264 # endif                         /* BN_MUL_COMBA */
265     {
266         p = &(t[n2 * 2]);
267         if (!zero)
268             bn_mul_recursive(&(t[n2]), t, &(t[n]), n, 0, 0, p);
269         else
270             memset(&t[n2], 0, sizeof(*t) * n2);
271         bn_mul_recursive(r, a, b, n, 0, 0, p);
272         bn_mul_recursive(&(r[n2]), &(a[n]), &(b[n]), n, dna, dnb, p);
273     }
274
275     /*-
276      * t[32] holds (a[0]-a[1])*(b[1]-b[0]), c1 is the sign
277      * r[10] holds (a[0]*b[0])
278      * r[32] holds (b[1]*b[1])
279      */
280
281     c1 = (int)(bn_add_words(t, r, &(r[n2]), n2));
282
283     if (neg) {                  /* if t[32] is negative */
284         c1 -= (int)(bn_sub_words(&(t[n2]), t, &(t[n2]), n2));
285     } else {
286         /* Might have a carry */
287         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), t, n2));
288     }
289
290     /*-
291      * t[32] holds (a[0]-a[1])*(b[1]-b[0])+(a[0]*b[0])+(a[1]*b[1])
292      * r[10] holds (a[0]*b[0])
293      * r[32] holds (b[1]*b[1])
294      * c1 holds the carry bits
295      */
296     c1 += (int)(bn_add_words(&(r[n]), &(r[n]), &(t[n2]), n2));
297     if (c1) {
298         p = &(r[n + n2]);
299         lo = *p;
300         ln = (lo + c1) & BN_MASK2;
301         *p = ln;
302
303         /*
304          * The overflow will stop before we over write words we should not
305          * overwrite
306          */
307         if (ln < (BN_ULONG)c1) {
308             do {
309                 p++;
310                 lo = *p;
311                 ln = (lo + 1) & BN_MASK2;
312                 *p = ln;
313             } while (ln == 0);
314         }
315     }
316 }
317
318 /*
319  * n+tn is the word length t needs to be n*4 is size, as does r
320  */
321 /* tnX may not be negative but less than n */
322 void bn_mul_part_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n,
323                            int tna, int tnb, BN_ULONG *t)
324 {
325     int i, j, n2 = n * 2;
326     int c1, c2, neg;
327     BN_ULONG ln, lo, *p;
328
329     if (n < 8) {
330         bn_mul_normal(r, a, n + tna, b, n + tnb);
331         return;
332     }
333
334     /* r=(a[0]-a[1])*(b[1]-b[0]) */
335     c1 = bn_cmp_part_words(a, &(a[n]), tna, n - tna);
336     c2 = bn_cmp_part_words(&(b[n]), b, tnb, tnb - n);
337     neg = 0;
338     switch (c1 * 3 + c2) {
339     case -4:
340         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
341         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
342         break;
343     case -3:
344     case -2:
345         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
346         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n); /* + */
347         neg = 1;
348         break;
349     case -1:
350     case 0:
351     case 1:
352     case 2:
353         bn_sub_part_words(t, a, &(a[n]), tna, n - tna); /* + */
354         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
355         neg = 1;
356         break;
357     case 3:
358     case 4:
359         bn_sub_part_words(t, a, &(a[n]), tna, n - tna);
360         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n);
361         break;
362     }
363     /*
364      * The zero case isn't yet implemented here. The speedup would probably
365      * be negligible.
366      */
367 # if 0
368     if (n == 4) {
369         bn_mul_comba4(&(t[n2]), t, &(t[n]));
370         bn_mul_comba4(r, a, b);
371         bn_mul_normal(&(r[n2]), &(a[n]), tn, &(b[n]), tn);
372         memset(&r[n2 + tn * 2], 0, sizeof(*r) * (n2 - tn * 2));
373     } else
374 # endif
375     if (n == 8) {
376         bn_mul_comba8(&(t[n2]), t, &(t[n]));
377         bn_mul_comba8(r, a, b);
378         bn_mul_normal(&(r[n2]), &(a[n]), tna, &(b[n]), tnb);
379         memset(&r[n2 + tna + tnb], 0, sizeof(*r) * (n2 - tna - tnb));
380     } else {
381         p = &(t[n2 * 2]);
382         bn_mul_recursive(&(t[n2]), t, &(t[n]), n, 0, 0, p);
383         bn_mul_recursive(r, a, b, n, 0, 0, p);
384         i = n / 2;
385         /*
386          * If there is only a bottom half to the number, just do it
387          */
388         if (tna > tnb)
389             j = tna - i;
390         else
391             j = tnb - i;
392         if (j == 0) {
393             bn_mul_recursive(&(r[n2]), &(a[n]), &(b[n]),
394                              i, tna - i, tnb - i, p);
395             memset(&r[n2 + i * 2], 0, sizeof(*r) * (n2 - i * 2));
396         } else if (j > 0) {     /* eg, n == 16, i == 8 and tn == 11 */
397             bn_mul_part_recursive(&(r[n2]), &(a[n]), &(b[n]),
398                                   i, tna - i, tnb - i, p);
399             memset(&(r[n2 + tna + tnb]), 0,
400                    sizeof(BN_ULONG) * (n2 - tna - tnb));
401         } else {                /* (j < 0) eg, n == 16, i == 8 and tn == 5 */
402
403             memset(&r[n2], 0, sizeof(*r) * n2);
404             if (tna < BN_MUL_RECURSIVE_SIZE_NORMAL
405                 && tnb < BN_MUL_RECURSIVE_SIZE_NORMAL) {
406                 bn_mul_normal(&(r[n2]), &(a[n]), tna, &(b[n]), tnb);
407             } else {
408                 for (;;) {
409                     i /= 2;
410                     /*
411                      * these simplified conditions work exclusively because
412                      * difference between tna and tnb is 1 or 0
413                      */
414                     if (i < tna || i < tnb) {
415                         bn_mul_part_recursive(&(r[n2]),
416                                               &(a[n]), &(b[n]),
417                                               i, tna - i, tnb - i, p);
418                         break;
419                     } else if (i == tna || i == tnb) {
420                         bn_mul_recursive(&(r[n2]),
421                                          &(a[n]), &(b[n]),
422                                          i, tna - i, tnb - i, p);
423                         break;
424                     }
425                 }
426             }
427         }
428     }
429
430     /*-
431      * t[32] holds (a[0]-a[1])*(b[1]-b[0]), c1 is the sign
432      * r[10] holds (a[0]*b[0])
433      * r[32] holds (b[1]*b[1])
434      */
435
436     c1 = (int)(bn_add_words(t, r, &(r[n2]), n2));
437
438     if (neg) {                  /* if t[32] is negative */
439         c1 -= (int)(bn_sub_words(&(t[n2]), t, &(t[n2]), n2));
440     } else {
441         /* Might have a carry */
442         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), t, n2));
443     }
444
445     /*-
446      * t[32] holds (a[0]-a[1])*(b[1]-b[0])+(a[0]*b[0])+(a[1]*b[1])
447      * r[10] holds (a[0]*b[0])
448      * r[32] holds (b[1]*b[1])
449      * c1 holds the carry bits
450      */
451     c1 += (int)(bn_add_words(&(r[n]), &(r[n]), &(t[n2]), n2));
452     if (c1) {
453         p = &(r[n + n2]);
454         lo = *p;
455         ln = (lo + c1) & BN_MASK2;
456         *p = ln;
457
458         /*
459          * The overflow will stop before we over write words we should not
460          * overwrite
461          */
462         if (ln < (BN_ULONG)c1) {
463             do {
464                 p++;
465                 lo = *p;
466                 ln = (lo + 1) & BN_MASK2;
467                 *p = ln;
468             } while (ln == 0);
469         }
470     }
471 }
472
473 /*-
474  * a and b must be the same size, which is n2.
475  * r needs to be n2 words and t needs to be n2*2
476  */
477 void bn_mul_low_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n2,
478                           BN_ULONG *t)
479 {
480     int n = n2 / 2;
481
482     bn_mul_recursive(r, a, b, n, 0, 0, &(t[0]));
483     if (n >= BN_MUL_LOW_RECURSIVE_SIZE_NORMAL) {
484         bn_mul_low_recursive(&(t[0]), &(a[0]), &(b[n]), n, &(t[n2]));
485         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
486         bn_mul_low_recursive(&(t[0]), &(a[n]), &(b[0]), n, &(t[n2]));
487         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
488     } else {
489         bn_mul_low_normal(&(t[0]), &(a[0]), &(b[n]), n);
490         bn_mul_low_normal(&(t[n]), &(a[n]), &(b[0]), n);
491         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
492         bn_add_words(&(r[n]), &(r[n]), &(t[n]), n);
493     }
494 }
495 #endif                          /* BN_RECURSION */
496
497 int BN_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
498 {
499     int ret = bn_mul_fixed_top(r, a, b, ctx);
500
501     bn_correct_top(r);
502     bn_check_top(r);
503
504     return ret;
505 }
506
507 int bn_mul_fixed_top(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
508 {
509     int ret = 0;
510     int top, al, bl;
511     BIGNUM *rr;
512 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
513     int i;
514 #endif
515 #ifdef BN_RECURSION
516     BIGNUM *t = NULL;
517     int j = 0, k;
518 #endif
519
520     bn_check_top(a);
521     bn_check_top(b);
522     bn_check_top(r);
523
524     al = a->top;
525     bl = b->top;
526
527     if ((al == 0) || (bl == 0)) {
528         BN_zero(r);
529         return 1;
530     }
531     top = al + bl;
532
533     BN_CTX_start(ctx);
534     if ((r == a) || (r == b)) {
535         if ((rr = BN_CTX_get(ctx)) == NULL)
536             goto err;
537     } else
538         rr = r;
539
540 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
541     i = al - bl;
542 #endif
543 #ifdef BN_MUL_COMBA
544     if (i == 0) {
545 # if 0
546         if (al == 4) {
547             if (bn_wexpand(rr, 8) == NULL)
548                 goto err;
549             rr->top = 8;
550             bn_mul_comba4(rr->d, a->d, b->d);
551             goto end;
552         }
553 # endif
554         if (al == 8) {
555             if (bn_wexpand(rr, 16) == NULL)
556                 goto err;
557             rr->top = 16;
558             bn_mul_comba8(rr->d, a->d, b->d);
559             goto end;
560         }
561     }
562 #endif                          /* BN_MUL_COMBA */
563 #ifdef BN_RECURSION
564     if ((al >= BN_MULL_SIZE_NORMAL) && (bl >= BN_MULL_SIZE_NORMAL)) {
565         if (i >= -1 && i <= 1) {
566             /*
567              * Find out the power of two lower or equal to the longest of the
568              * two numbers
569              */
570             if (i >= 0) {
571                 j = BN_num_bits_word((BN_ULONG)al);
572             }
573             if (i == -1) {
574                 j = BN_num_bits_word((BN_ULONG)bl);
575             }
576             j = 1 << (j - 1);
577             assert(j <= al || j <= bl);
578             k = j + j;
579             t = BN_CTX_get(ctx);
580             if (t == NULL)
581                 goto err;
582             if (al > j || bl > j) {
583                 if (bn_wexpand(t, k * 4) == NULL)
584                     goto err;
585                 if (bn_wexpand(rr, k * 4) == NULL)
586                     goto err;
587                 bn_mul_part_recursive(rr->d, a->d, b->d,
588                                       j, al - j, bl - j, t->d);
589             } else {            /* al <= j || bl <= j */
590
591                 if (bn_wexpand(t, k * 2) == NULL)
592                     goto err;
593                 if (bn_wexpand(rr, k * 2) == NULL)
594                     goto err;
595                 bn_mul_recursive(rr->d, a->d, b->d, j, al - j, bl - j, t->d);
596             }
597             rr->top = top;
598             goto end;
599         }
600     }
601 #endif                          /* BN_RECURSION */
602     if (bn_wexpand(rr, top) == NULL)
603         goto err;
604     rr->top = top;
605     bn_mul_normal(rr->d, a->d, al, b->d, bl);
606
607 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
608  end:
609 #endif
610     rr->neg = a->neg ^ b->neg;
611     rr->flags |= BN_FLG_FIXED_TOP;
612     if (r != rr && BN_copy(r, rr) == NULL)
613         goto err;
614
615     ret = 1;
616  err:
617     bn_check_top(r);
618     BN_CTX_end(ctx);
619     return ret;
620 }
621
622 void bn_mul_normal(BN_ULONG *r, BN_ULONG *a, int na, BN_ULONG *b, int nb)
623 {
624     BN_ULONG *rr;
625
626     if (na < nb) {
627         int itmp;
628         BN_ULONG *ltmp;
629
630         itmp = na;
631         na = nb;
632         nb = itmp;
633         ltmp = a;
634         a = b;
635         b = ltmp;
636
637     }
638     rr = &(r[na]);
639     if (nb <= 0) {
640         (void)bn_mul_words(r, a, na, 0);
641         return;
642     } else
643         rr[0] = bn_mul_words(r, a, na, b[0]);
644
645     for (;;) {
646         if (--nb <= 0)
647             return;
648         rr[1] = bn_mul_add_words(&(r[1]), a, na, b[1]);
649         if (--nb <= 0)
650             return;
651         rr[2] = bn_mul_add_words(&(r[2]), a, na, b[2]);
652         if (--nb <= 0)
653             return;
654         rr[3] = bn_mul_add_words(&(r[3]), a, na, b[3]);
655         if (--nb <= 0)
656             return;
657         rr[4] = bn_mul_add_words(&(r[4]), a, na, b[4]);
658         rr += 4;
659         r += 4;
660         b += 4;
661     }
662 }
663
664 void bn_mul_low_normal(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n)
665 {
666     bn_mul_words(r, a, n, b[0]);
667
668     for (;;) {
669         if (--n <= 0)
670             return;
671         bn_mul_add_words(&(r[1]), a, n, b[1]);
672         if (--n <= 0)
673             return;
674         bn_mul_add_words(&(r[2]), a, n, b[2]);
675         if (--n <= 0)
676             return;
677         bn_mul_add_words(&(r[3]), a, n, b[3]);
678         if (--n <= 0)
679             return;
680         bn_mul_add_words(&(r[4]), a, n, b[4]);
681         r += 4;
682         b += 4;
683     }
684 }