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