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