RT4593: Add space after comma (doc nits)
[openssl.git] / crypto / bn / bn_mul.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 <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                 case 2:
119                     r[2] = a[2];
120                     if (--dl <= 0)
121                         break;
122                 case 3:
123                     r[3] = a[3];
124                     if (--dl <= 0)
125                         break;
126                 }
127                 a += 4;
128                 r += 4;
129             }
130         }
131         if (dl > 0) {
132             for (;;) {
133                 r[0] = a[0];
134                 if (--dl <= 0)
135                     break;
136                 r[1] = a[1];
137                 if (--dl <= 0)
138                     break;
139                 r[2] = a[2];
140                 if (--dl <= 0)
141                     break;
142                 r[3] = a[3];
143                 if (--dl <= 0)
144                     break;
145
146                 a += 4;
147                 r += 4;
148             }
149         }
150     }
151     return c;
152 }
153 #endif
154
155 BN_ULONG bn_add_part_words(BN_ULONG *r,
156                            const BN_ULONG *a, const BN_ULONG *b,
157                            int cl, int dl)
158 {
159     BN_ULONG c, l, t;
160
161     assert(cl >= 0);
162     c = bn_add_words(r, a, b, cl);
163
164     if (dl == 0)
165         return c;
166
167     r += cl;
168     a += cl;
169     b += cl;
170
171     if (dl < 0) {
172         int save_dl = dl;
173         while (c) {
174             l = (c + b[0]) & BN_MASK2;
175             c = (l < c);
176             r[0] = l;
177             if (++dl >= 0)
178                 break;
179
180             l = (c + b[1]) & BN_MASK2;
181             c = (l < c);
182             r[1] = l;
183             if (++dl >= 0)
184                 break;
185
186             l = (c + b[2]) & BN_MASK2;
187             c = (l < c);
188             r[2] = l;
189             if (++dl >= 0)
190                 break;
191
192             l = (c + b[3]) & BN_MASK2;
193             c = (l < c);
194             r[3] = l;
195             if (++dl >= 0)
196                 break;
197
198             save_dl = dl;
199             b += 4;
200             r += 4;
201         }
202         if (dl < 0) {
203             if (save_dl < dl) {
204                 switch (dl - save_dl) {
205                 case 1:
206                     r[1] = b[1];
207                     if (++dl >= 0)
208                         break;
209                 case 2:
210                     r[2] = b[2];
211                     if (++dl >= 0)
212                         break;
213                 case 3:
214                     r[3] = b[3];
215                     if (++dl >= 0)
216                         break;
217                 }
218                 b += 4;
219                 r += 4;
220             }
221         }
222         if (dl < 0) {
223             for (;;) {
224                 r[0] = b[0];
225                 if (++dl >= 0)
226                     break;
227                 r[1] = b[1];
228                 if (++dl >= 0)
229                     break;
230                 r[2] = b[2];
231                 if (++dl >= 0)
232                     break;
233                 r[3] = b[3];
234                 if (++dl >= 0)
235                     break;
236
237                 b += 4;
238                 r += 4;
239             }
240         }
241     } else {
242         int save_dl = dl;
243         while (c) {
244             t = (a[0] + c) & BN_MASK2;
245             c = (t < c);
246             r[0] = t;
247             if (--dl <= 0)
248                 break;
249
250             t = (a[1] + c) & BN_MASK2;
251             c = (t < c);
252             r[1] = t;
253             if (--dl <= 0)
254                 break;
255
256             t = (a[2] + c) & BN_MASK2;
257             c = (t < c);
258             r[2] = t;
259             if (--dl <= 0)
260                 break;
261
262             t = (a[3] + c) & BN_MASK2;
263             c = (t < c);
264             r[3] = t;
265             if (--dl <= 0)
266                 break;
267
268             save_dl = dl;
269             a += 4;
270             r += 4;
271         }
272         if (dl > 0) {
273             if (save_dl > dl) {
274                 switch (save_dl - dl) {
275                 case 1:
276                     r[1] = a[1];
277                     if (--dl <= 0)
278                         break;
279                 case 2:
280                     r[2] = a[2];
281                     if (--dl <= 0)
282                         break;
283                 case 3:
284                     r[3] = a[3];
285                     if (--dl <= 0)
286                         break;
287                 }
288                 a += 4;
289                 r += 4;
290             }
291         }
292         if (dl > 0) {
293             for (;;) {
294                 r[0] = a[0];
295                 if (--dl <= 0)
296                     break;
297                 r[1] = a[1];
298                 if (--dl <= 0)
299                     break;
300                 r[2] = a[2];
301                 if (--dl <= 0)
302                     break;
303                 r[3] = a[3];
304                 if (--dl <= 0)
305                     break;
306
307                 a += 4;
308                 r += 4;
309             }
310         }
311     }
312     return c;
313 }
314
315 #ifdef BN_RECURSION
316 /*
317  * Karatsuba recursive multiplication algorithm (cf. Knuth, The Art of
318  * Computer Programming, Vol. 2)
319  */
320
321 /*-
322  * r is 2*n2 words in size,
323  * a and b are both n2 words in size.
324  * n2 must be a power of 2.
325  * We multiply and return the result.
326  * t must be 2*n2 words in size
327  * We calculate
328  * a[0]*b[0]
329  * a[0]*b[0]+a[1]*b[1]+(a[0]-a[1])*(b[1]-b[0])
330  * a[1]*b[1]
331  */
332 /* dnX may not be positive, but n2/2+dnX has to be */
333 void bn_mul_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n2,
334                       int dna, int dnb, BN_ULONG *t)
335 {
336     int n = n2 / 2, c1, c2;
337     int tna = n + dna, tnb = n + dnb;
338     unsigned int neg, zero;
339     BN_ULONG ln, lo, *p;
340
341 # ifdef BN_MUL_COMBA
342 #  if 0
343     if (n2 == 4) {
344         bn_mul_comba4(r, a, b);
345         return;
346     }
347 #  endif
348     /*
349      * Only call bn_mul_comba 8 if n2 == 8 and the two arrays are complete
350      * [steve]
351      */
352     if (n2 == 8 && dna == 0 && dnb == 0) {
353         bn_mul_comba8(r, a, b);
354         return;
355     }
356 # endif                         /* BN_MUL_COMBA */
357     /* Else do normal multiply */
358     if (n2 < BN_MUL_RECURSIVE_SIZE_NORMAL) {
359         bn_mul_normal(r, a, n2 + dna, b, n2 + dnb);
360         if ((dna + dnb) < 0)
361             memset(&r[2 * n2 + dna + dnb], 0,
362                    sizeof(BN_ULONG) * -(dna + dnb));
363         return;
364     }
365     /* r=(a[0]-a[1])*(b[1]-b[0]) */
366     c1 = bn_cmp_part_words(a, &(a[n]), tna, n - tna);
367     c2 = bn_cmp_part_words(&(b[n]), b, tnb, tnb - n);
368     zero = neg = 0;
369     switch (c1 * 3 + c2) {
370     case -4:
371         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
372         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
373         break;
374     case -3:
375         zero = 1;
376         break;
377     case -2:
378         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
379         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n); /* + */
380         neg = 1;
381         break;
382     case -1:
383     case 0:
384     case 1:
385         zero = 1;
386         break;
387     case 2:
388         bn_sub_part_words(t, a, &(a[n]), tna, n - tna); /* + */
389         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
390         neg = 1;
391         break;
392     case 3:
393         zero = 1;
394         break;
395     case 4:
396         bn_sub_part_words(t, a, &(a[n]), tna, n - tna);
397         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n);
398         break;
399     }
400
401 # ifdef BN_MUL_COMBA
402     if (n == 4 && dna == 0 && dnb == 0) { /* XXX: bn_mul_comba4 could take
403                                            * extra args to do this well */
404         if (!zero)
405             bn_mul_comba4(&(t[n2]), t, &(t[n]));
406         else
407             memset(&t[n2], 0, sizeof(*t) * 8);
408
409         bn_mul_comba4(r, a, b);
410         bn_mul_comba4(&(r[n2]), &(a[n]), &(b[n]));
411     } else if (n == 8 && dna == 0 && dnb == 0) { /* XXX: bn_mul_comba8 could
412                                                   * take extra args to do
413                                                   * this well */
414         if (!zero)
415             bn_mul_comba8(&(t[n2]), t, &(t[n]));
416         else
417             memset(&t[n2], 0, sizeof(*t) * 16);
418
419         bn_mul_comba8(r, a, b);
420         bn_mul_comba8(&(r[n2]), &(a[n]), &(b[n]));
421     } else
422 # endif                         /* BN_MUL_COMBA */
423     {
424         p = &(t[n2 * 2]);
425         if (!zero)
426             bn_mul_recursive(&(t[n2]), t, &(t[n]), n, 0, 0, p);
427         else
428             memset(&t[n2], 0, sizeof(*t) * n2);
429         bn_mul_recursive(r, a, b, n, 0, 0, p);
430         bn_mul_recursive(&(r[n2]), &(a[n]), &(b[n]), n, dna, dnb, p);
431     }
432
433     /*-
434      * t[32] holds (a[0]-a[1])*(b[1]-b[0]), c1 is the sign
435      * r[10] holds (a[0]*b[0])
436      * r[32] holds (b[1]*b[1])
437      */
438
439     c1 = (int)(bn_add_words(t, r, &(r[n2]), n2));
440
441     if (neg) {                  /* if t[32] is negative */
442         c1 -= (int)(bn_sub_words(&(t[n2]), t, &(t[n2]), n2));
443     } else {
444         /* Might have a carry */
445         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), t, n2));
446     }
447
448     /*-
449      * t[32] holds (a[0]-a[1])*(b[1]-b[0])+(a[0]*b[0])+(a[1]*b[1])
450      * r[10] holds (a[0]*b[0])
451      * r[32] holds (b[1]*b[1])
452      * c1 holds the carry bits
453      */
454     c1 += (int)(bn_add_words(&(r[n]), &(r[n]), &(t[n2]), n2));
455     if (c1) {
456         p = &(r[n + n2]);
457         lo = *p;
458         ln = (lo + c1) & BN_MASK2;
459         *p = ln;
460
461         /*
462          * The overflow will stop before we over write words we should not
463          * overwrite
464          */
465         if (ln < (BN_ULONG)c1) {
466             do {
467                 p++;
468                 lo = *p;
469                 ln = (lo + 1) & BN_MASK2;
470                 *p = ln;
471             } while (ln == 0);
472         }
473     }
474 }
475
476 /*
477  * n+tn is the word length t needs to be n*4 is size, as does r
478  */
479 /* tnX may not be negative but less than n */
480 void bn_mul_part_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n,
481                            int tna, int tnb, BN_ULONG *t)
482 {
483     int i, j, n2 = n * 2;
484     int c1, c2, neg;
485     BN_ULONG ln, lo, *p;
486
487     if (n < 8) {
488         bn_mul_normal(r, a, n + tna, b, n + tnb);
489         return;
490     }
491
492     /* r=(a[0]-a[1])*(b[1]-b[0]) */
493     c1 = bn_cmp_part_words(a, &(a[n]), tna, n - tna);
494     c2 = bn_cmp_part_words(&(b[n]), b, tnb, tnb - n);
495     neg = 0;
496     switch (c1 * 3 + c2) {
497     case -4:
498         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
499         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
500         break;
501     case -3:
502         /* break; */
503     case -2:
504         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
505         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n); /* + */
506         neg = 1;
507         break;
508     case -1:
509     case 0:
510     case 1:
511         /* break; */
512     case 2:
513         bn_sub_part_words(t, a, &(a[n]), tna, n - tna); /* + */
514         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
515         neg = 1;
516         break;
517     case 3:
518         /* break; */
519     case 4:
520         bn_sub_part_words(t, a, &(a[n]), tna, n - tna);
521         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n);
522         break;
523     }
524     /*
525      * The zero case isn't yet implemented here. The speedup would probably
526      * be negligible.
527      */
528 # if 0
529     if (n == 4) {
530         bn_mul_comba4(&(t[n2]), t, &(t[n]));
531         bn_mul_comba4(r, a, b);
532         bn_mul_normal(&(r[n2]), &(a[n]), tn, &(b[n]), tn);
533         memset(&r[n2 + tn * 2], 0, sizeof(*r) * (n2 - tn * 2));
534     } else
535 # endif
536     if (n == 8) {
537         bn_mul_comba8(&(t[n2]), t, &(t[n]));
538         bn_mul_comba8(r, a, b);
539         bn_mul_normal(&(r[n2]), &(a[n]), tna, &(b[n]), tnb);
540         memset(&r[n2 + tna + tnb], 0, sizeof(*r) * (n2 - tna - tnb));
541     } else {
542         p = &(t[n2 * 2]);
543         bn_mul_recursive(&(t[n2]), t, &(t[n]), n, 0, 0, p);
544         bn_mul_recursive(r, a, b, n, 0, 0, p);
545         i = n / 2;
546         /*
547          * If there is only a bottom half to the number, just do it
548          */
549         if (tna > tnb)
550             j = tna - i;
551         else
552             j = tnb - i;
553         if (j == 0) {
554             bn_mul_recursive(&(r[n2]), &(a[n]), &(b[n]),
555                              i, tna - i, tnb - i, p);
556             memset(&r[n2 + i * 2], 0, sizeof(*r) * (n2 - i * 2));
557         } else if (j > 0) {     /* eg, n == 16, i == 8 and tn == 11 */
558             bn_mul_part_recursive(&(r[n2]), &(a[n]), &(b[n]),
559                                   i, tna - i, tnb - i, p);
560             memset(&(r[n2 + tna + tnb]), 0,
561                    sizeof(BN_ULONG) * (n2 - tna - tnb));
562         } else {                /* (j < 0) eg, n == 16, i == 8 and tn == 5 */
563
564             memset(&r[n2], 0, sizeof(*r) * n2);
565             if (tna < BN_MUL_RECURSIVE_SIZE_NORMAL
566                 && tnb < BN_MUL_RECURSIVE_SIZE_NORMAL) {
567                 bn_mul_normal(&(r[n2]), &(a[n]), tna, &(b[n]), tnb);
568             } else {
569                 for (;;) {
570                     i /= 2;
571                     /*
572                      * these simplified conditions work exclusively because
573                      * difference between tna and tnb is 1 or 0
574                      */
575                     if (i < tna || i < tnb) {
576                         bn_mul_part_recursive(&(r[n2]),
577                                               &(a[n]), &(b[n]),
578                                               i, tna - i, tnb - i, p);
579                         break;
580                     } else if (i == tna || i == tnb) {
581                         bn_mul_recursive(&(r[n2]),
582                                          &(a[n]), &(b[n]),
583                                          i, tna - i, tnb - i, p);
584                         break;
585                     }
586                 }
587             }
588         }
589     }
590
591     /*-
592      * t[32] holds (a[0]-a[1])*(b[1]-b[0]), c1 is the sign
593      * r[10] holds (a[0]*b[0])
594      * r[32] holds (b[1]*b[1])
595      */
596
597     c1 = (int)(bn_add_words(t, r, &(r[n2]), n2));
598
599     if (neg) {                  /* if t[32] is negative */
600         c1 -= (int)(bn_sub_words(&(t[n2]), t, &(t[n2]), n2));
601     } else {
602         /* Might have a carry */
603         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), t, n2));
604     }
605
606     /*-
607      * t[32] holds (a[0]-a[1])*(b[1]-b[0])+(a[0]*b[0])+(a[1]*b[1])
608      * r[10] holds (a[0]*b[0])
609      * r[32] holds (b[1]*b[1])
610      * c1 holds the carry bits
611      */
612     c1 += (int)(bn_add_words(&(r[n]), &(r[n]), &(t[n2]), n2));
613     if (c1) {
614         p = &(r[n + n2]);
615         lo = *p;
616         ln = (lo + c1) & BN_MASK2;
617         *p = ln;
618
619         /*
620          * The overflow will stop before we over write words we should not
621          * overwrite
622          */
623         if (ln < (BN_ULONG)c1) {
624             do {
625                 p++;
626                 lo = *p;
627                 ln = (lo + 1) & BN_MASK2;
628                 *p = ln;
629             } while (ln == 0);
630         }
631     }
632 }
633
634 /*-
635  * a and b must be the same size, which is n2.
636  * r needs to be n2 words and t needs to be n2*2
637  */
638 void bn_mul_low_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n2,
639                           BN_ULONG *t)
640 {
641     int n = n2 / 2;
642
643     bn_mul_recursive(r, a, b, n, 0, 0, &(t[0]));
644     if (n >= BN_MUL_LOW_RECURSIVE_SIZE_NORMAL) {
645         bn_mul_low_recursive(&(t[0]), &(a[0]), &(b[n]), n, &(t[n2]));
646         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
647         bn_mul_low_recursive(&(t[0]), &(a[n]), &(b[0]), n, &(t[n2]));
648         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
649     } else {
650         bn_mul_low_normal(&(t[0]), &(a[0]), &(b[n]), n);
651         bn_mul_low_normal(&(t[n]), &(a[n]), &(b[0]), n);
652         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
653         bn_add_words(&(r[n]), &(r[n]), &(t[n]), n);
654     }
655 }
656
657 /*-
658  * a and b must be the same size, which is n2.
659  * r needs to be n2 words and t needs to be n2*2
660  * l is the low words of the output.
661  * t needs to be n2*3
662  */
663 void bn_mul_high(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, BN_ULONG *l, int n2,
664                  BN_ULONG *t)
665 {
666     int i, n;
667     int c1, c2;
668     int neg, oneg, zero;
669     BN_ULONG ll, lc, *lp, *mp;
670
671     n = n2 / 2;
672
673     /* Calculate (al-ah)*(bh-bl) */
674     neg = zero = 0;
675     c1 = bn_cmp_words(&(a[0]), &(a[n]), n);
676     c2 = bn_cmp_words(&(b[n]), &(b[0]), n);
677     switch (c1 * 3 + c2) {
678     case -4:
679         bn_sub_words(&(r[0]), &(a[n]), &(a[0]), n);
680         bn_sub_words(&(r[n]), &(b[0]), &(b[n]), n);
681         break;
682     case -3:
683         zero = 1;
684         break;
685     case -2:
686         bn_sub_words(&(r[0]), &(a[n]), &(a[0]), n);
687         bn_sub_words(&(r[n]), &(b[n]), &(b[0]), n);
688         neg = 1;
689         break;
690     case -1:
691     case 0:
692     case 1:
693         zero = 1;
694         break;
695     case 2:
696         bn_sub_words(&(r[0]), &(a[0]), &(a[n]), n);
697         bn_sub_words(&(r[n]), &(b[0]), &(b[n]), n);
698         neg = 1;
699         break;
700     case 3:
701         zero = 1;
702         break;
703     case 4:
704         bn_sub_words(&(r[0]), &(a[0]), &(a[n]), n);
705         bn_sub_words(&(r[n]), &(b[n]), &(b[0]), n);
706         break;
707     }
708
709     oneg = neg;
710     /* t[10] = (a[0]-a[1])*(b[1]-b[0]) */
711     /* r[10] = (a[1]*b[1]) */
712 # ifdef BN_MUL_COMBA
713     if (n == 8) {
714         bn_mul_comba8(&(t[0]), &(r[0]), &(r[n]));
715         bn_mul_comba8(r, &(a[n]), &(b[n]));
716     } else
717 # endif
718     {
719         bn_mul_recursive(&(t[0]), &(r[0]), &(r[n]), n, 0, 0, &(t[n2]));
720         bn_mul_recursive(r, &(a[n]), &(b[n]), n, 0, 0, &(t[n2]));
721     }
722
723     /*-
724      * s0 == low(al*bl)
725      * s1 == low(ah*bh)+low((al-ah)*(bh-bl))+low(al*bl)+high(al*bl)
726      * We know s0 and s1 so the only unknown is high(al*bl)
727      * high(al*bl) == s1 - low(ah*bh+s0+(al-ah)*(bh-bl))
728      * high(al*bl) == s1 - (r[0]+l[0]+t[0])
729      */
730     if (l != NULL) {
731         lp = &(t[n2 + n]);
732         c1 = (int)(bn_add_words(lp, &(r[0]), &(l[0]), n));
733     } else {
734         c1 = 0;
735         lp = &(r[0]);
736     }
737
738     if (neg)
739         neg = (int)(bn_sub_words(&(t[n2]), lp, &(t[0]), n));
740     else {
741         bn_add_words(&(t[n2]), lp, &(t[0]), n);
742         neg = 0;
743     }
744
745     if (l != NULL) {
746         bn_sub_words(&(t[n2 + n]), &(l[n]), &(t[n2]), n);
747     } else {
748         lp = &(t[n2 + n]);
749         mp = &(t[n2]);
750         for (i = 0; i < n; i++)
751             lp[i] = ((~mp[i]) + 1) & BN_MASK2;
752     }
753
754     /*-
755      * s[0] = low(al*bl)
756      * t[3] = high(al*bl)
757      * t[10] = (a[0]-a[1])*(b[1]-b[0]) neg is the sign
758      * r[10] = (a[1]*b[1])
759      */
760     /*-
761      * R[10] = al*bl
762      * R[21] = al*bl + ah*bh + (a[0]-a[1])*(b[1]-b[0])
763      * R[32] = ah*bh
764      */
765     /*-
766      * R[1]=t[3]+l[0]+r[0](+-)t[0] (have carry/borrow)
767      * R[2]=r[0]+t[3]+r[1](+-)t[1] (have carry/borrow)
768      * R[3]=r[1]+(carry/borrow)
769      */
770     if (l != NULL) {
771         lp = &(t[n2]);
772         c1 = (int)(bn_add_words(lp, &(t[n2 + n]), &(l[0]), n));
773     } else {
774         lp = &(t[n2 + n]);
775         c1 = 0;
776     }
777     c1 += (int)(bn_add_words(&(t[n2]), lp, &(r[0]), n));
778     if (oneg)
779         c1 -= (int)(bn_sub_words(&(t[n2]), &(t[n2]), &(t[0]), n));
780     else
781         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), &(t[0]), n));
782
783     c2 = (int)(bn_add_words(&(r[0]), &(r[0]), &(t[n2 + n]), n));
784     c2 += (int)(bn_add_words(&(r[0]), &(r[0]), &(r[n]), n));
785     if (oneg)
786         c2 -= (int)(bn_sub_words(&(r[0]), &(r[0]), &(t[n]), n));
787     else
788         c2 += (int)(bn_add_words(&(r[0]), &(r[0]), &(t[n]), n));
789
790     if (c1 != 0) {              /* Add starting at r[0], could be +ve or -ve */
791         i = 0;
792         if (c1 > 0) {
793             lc = c1;
794             do {
795                 ll = (r[i] + lc) & BN_MASK2;
796                 r[i++] = ll;
797                 lc = (lc > ll);
798             } while (lc);
799         } else {
800             lc = -c1;
801             do {
802                 ll = r[i];
803                 r[i++] = (ll - lc) & BN_MASK2;
804                 lc = (lc > ll);
805             } while (lc);
806         }
807     }
808     if (c2 != 0) {              /* Add starting at r[1] */
809         i = n;
810         if (c2 > 0) {
811             lc = c2;
812             do {
813                 ll = (r[i] + lc) & BN_MASK2;
814                 r[i++] = ll;
815                 lc = (lc > ll);
816             } while (lc);
817         } else {
818             lc = -c2;
819             do {
820                 ll = r[i];
821                 r[i++] = (ll - lc) & BN_MASK2;
822                 lc = (lc > ll);
823             } while (lc);
824         }
825     }
826 }
827 #endif                          /* BN_RECURSION */
828
829 int BN_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
830 {
831     int ret = 0;
832     int top, al, bl;
833     BIGNUM *rr;
834 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
835     int i;
836 #endif
837 #ifdef BN_RECURSION
838     BIGNUM *t = NULL;
839     int j = 0, k;
840 #endif
841
842     bn_check_top(a);
843     bn_check_top(b);
844     bn_check_top(r);
845
846     al = a->top;
847     bl = b->top;
848
849     if ((al == 0) || (bl == 0)) {
850         BN_zero(r);
851         return (1);
852     }
853     top = al + bl;
854
855     BN_CTX_start(ctx);
856     if ((r == a) || (r == b)) {
857         if ((rr = BN_CTX_get(ctx)) == NULL)
858             goto err;
859     } else
860         rr = r;
861     rr->neg = a->neg ^ b->neg;
862
863 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
864     i = al - bl;
865 #endif
866 #ifdef BN_MUL_COMBA
867     if (i == 0) {
868 # if 0
869         if (al == 4) {
870             if (bn_wexpand(rr, 8) == NULL)
871                 goto err;
872             rr->top = 8;
873             bn_mul_comba4(rr->d, a->d, b->d);
874             goto end;
875         }
876 # endif
877         if (al == 8) {
878             if (bn_wexpand(rr, 16) == NULL)
879                 goto err;
880             rr->top = 16;
881             bn_mul_comba8(rr->d, a->d, b->d);
882             goto end;
883         }
884     }
885 #endif                          /* BN_MUL_COMBA */
886 #ifdef BN_RECURSION
887     if ((al >= BN_MULL_SIZE_NORMAL) && (bl >= BN_MULL_SIZE_NORMAL)) {
888         if (i >= -1 && i <= 1) {
889             /*
890              * Find out the power of two lower or equal to the longest of the
891              * two numbers
892              */
893             if (i >= 0) {
894                 j = BN_num_bits_word((BN_ULONG)al);
895             }
896             if (i == -1) {
897                 j = BN_num_bits_word((BN_ULONG)bl);
898             }
899             j = 1 << (j - 1);
900             assert(j <= al || j <= bl);
901             k = j + j;
902             t = BN_CTX_get(ctx);
903             if (t == NULL)
904                 goto err;
905             if (al > j || bl > j) {
906                 if (bn_wexpand(t, k * 4) == NULL)
907                     goto err;
908                 if (bn_wexpand(rr, k * 4) == NULL)
909                     goto err;
910                 bn_mul_part_recursive(rr->d, a->d, b->d,
911                                       j, al - j, bl - j, t->d);
912             } else {            /* al <= j || bl <= j */
913
914                 if (bn_wexpand(t, k * 2) == NULL)
915                     goto err;
916                 if (bn_wexpand(rr, k * 2) == NULL)
917                     goto err;
918                 bn_mul_recursive(rr->d, a->d, b->d, j, al - j, bl - j, t->d);
919             }
920             rr->top = top;
921             goto end;
922         }
923 # if 0
924         if (i == 1 && !BN_get_flags(b, BN_FLG_STATIC_DATA)) {
925             BIGNUM *tmp_bn = (BIGNUM *)b;
926             if (bn_wexpand(tmp_bn, al) == NULL)
927                 goto err;
928             tmp_bn->d[bl] = 0;
929             bl++;
930             i--;
931         } else if (i == -1 && !BN_get_flags(a, BN_FLG_STATIC_DATA)) {
932             BIGNUM *tmp_bn = (BIGNUM *)a;
933             if (bn_wexpand(tmp_bn, bl) == NULL)
934                 goto err;
935             tmp_bn->d[al] = 0;
936             al++;
937             i++;
938         }
939         if (i == 0) {
940             /* symmetric and > 4 */
941             /* 16 or larger */
942             j = BN_num_bits_word((BN_ULONG)al);
943             j = 1 << (j - 1);
944             k = j + j;
945             t = BN_CTX_get(ctx);
946             if (al == j) {      /* exact multiple */
947                 if (bn_wexpand(t, k * 2) == NULL)
948                     goto err;
949                 if (bn_wexpand(rr, k * 2) == NULL)
950                     goto err;
951                 bn_mul_recursive(rr->d, a->d, b->d, al, t->d);
952             } else {
953                 if (bn_wexpand(t, k * 4) == NULL)
954                     goto err;
955                 if (bn_wexpand(rr, k * 4) == NULL)
956                     goto err;
957                 bn_mul_part_recursive(rr->d, a->d, b->d, al - j, j, t->d);
958             }
959             rr->top = top;
960             goto end;
961         }
962 # endif
963     }
964 #endif                          /* BN_RECURSION */
965     if (bn_wexpand(rr, top) == NULL)
966         goto err;
967     rr->top = top;
968     bn_mul_normal(rr->d, a->d, al, b->d, bl);
969
970 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
971  end:
972 #endif
973     bn_correct_top(rr);
974     if (r != rr)
975         BN_copy(r, rr);
976     ret = 1;
977  err:
978     bn_check_top(r);
979     BN_CTX_end(ctx);
980     return (ret);
981 }
982
983 void bn_mul_normal(BN_ULONG *r, BN_ULONG *a, int na, BN_ULONG *b, int nb)
984 {
985     BN_ULONG *rr;
986
987     if (na < nb) {
988         int itmp;
989         BN_ULONG *ltmp;
990
991         itmp = na;
992         na = nb;
993         nb = itmp;
994         ltmp = a;
995         a = b;
996         b = ltmp;
997
998     }
999     rr = &(r[na]);
1000     if (nb <= 0) {
1001         (void)bn_mul_words(r, a, na, 0);
1002         return;
1003     } else
1004         rr[0] = bn_mul_words(r, a, na, b[0]);
1005
1006     for (;;) {
1007         if (--nb <= 0)
1008             return;
1009         rr[1] = bn_mul_add_words(&(r[1]), a, na, b[1]);
1010         if (--nb <= 0)
1011             return;
1012         rr[2] = bn_mul_add_words(&(r[2]), a, na, b[2]);
1013         if (--nb <= 0)
1014             return;
1015         rr[3] = bn_mul_add_words(&(r[3]), a, na, b[3]);
1016         if (--nb <= 0)
1017             return;
1018         rr[4] = bn_mul_add_words(&(r[4]), a, na, b[4]);
1019         rr += 4;
1020         r += 4;
1021         b += 4;
1022     }
1023 }
1024
1025 void bn_mul_low_normal(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n)
1026 {
1027     bn_mul_words(r, a, n, b[0]);
1028
1029     for (;;) {
1030         if (--n <= 0)
1031             return;
1032         bn_mul_add_words(&(r[1]), a, n, b[1]);
1033         if (--n <= 0)
1034             return;
1035         bn_mul_add_words(&(r[2]), a, n, b[2]);
1036         if (--n <= 0)
1037             return;
1038         bn_mul_add_words(&(r[3]), a, n, b[3]);
1039         if (--n <= 0)
1040             return;
1041         bn_mul_add_words(&(r[4]), a, n, b[4]);
1042         r += 4;
1043         b += 4;
1044     }
1045 }