99402fc4c93998c041b3ee10b87467ab54c047d3
[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 #ifdef BN_RECURSION
156 /*
157  * Karatsuba recursive multiplication algorithm (cf. Knuth, The Art of
158  * Computer Programming, Vol. 2)
159  */
160
161 /*-
162  * r is 2*n2 words in size,
163  * a and b are both n2 words in size.
164  * n2 must be a power of 2.
165  * We multiply and return the result.
166  * t must be 2*n2 words in size
167  * We calculate
168  * a[0]*b[0]
169  * a[0]*b[0]+a[1]*b[1]+(a[0]-a[1])*(b[1]-b[0])
170  * a[1]*b[1]
171  */
172 /* dnX may not be positive, but n2/2+dnX has to be */
173 void bn_mul_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n2,
174                       int dna, int dnb, BN_ULONG *t)
175 {
176     int n = n2 / 2, c1, c2;
177     int tna = n + dna, tnb = n + dnb;
178     unsigned int neg, zero;
179     BN_ULONG ln, lo, *p;
180
181 # ifdef BN_MUL_COMBA
182 #  if 0
183     if (n2 == 4) {
184         bn_mul_comba4(r, a, b);
185         return;
186     }
187 #  endif
188     /*
189      * Only call bn_mul_comba 8 if n2 == 8 and the two arrays are complete
190      * [steve]
191      */
192     if (n2 == 8 && dna == 0 && dnb == 0) {
193         bn_mul_comba8(r, a, b);
194         return;
195     }
196 # endif                         /* BN_MUL_COMBA */
197     /* Else do normal multiply */
198     if (n2 < BN_MUL_RECURSIVE_SIZE_NORMAL) {
199         bn_mul_normal(r, a, n2 + dna, b, n2 + dnb);
200         if ((dna + dnb) < 0)
201             memset(&r[2 * n2 + dna + dnb], 0,
202                    sizeof(BN_ULONG) * -(dna + dnb));
203         return;
204     }
205     /* r=(a[0]-a[1])*(b[1]-b[0]) */
206     c1 = bn_cmp_part_words(a, &(a[n]), tna, n - tna);
207     c2 = bn_cmp_part_words(&(b[n]), b, tnb, tnb - n);
208     zero = neg = 0;
209     switch (c1 * 3 + c2) {
210     case -4:
211         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
212         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
213         break;
214     case -3:
215         zero = 1;
216         break;
217     case -2:
218         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
219         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n); /* + */
220         neg = 1;
221         break;
222     case -1:
223     case 0:
224     case 1:
225         zero = 1;
226         break;
227     case 2:
228         bn_sub_part_words(t, a, &(a[n]), tna, n - tna); /* + */
229         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
230         neg = 1;
231         break;
232     case 3:
233         zero = 1;
234         break;
235     case 4:
236         bn_sub_part_words(t, a, &(a[n]), tna, n - tna);
237         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n);
238         break;
239     }
240
241 # ifdef BN_MUL_COMBA
242     if (n == 4 && dna == 0 && dnb == 0) { /* XXX: bn_mul_comba4 could take
243                                            * extra args to do this well */
244         if (!zero)
245             bn_mul_comba4(&(t[n2]), t, &(t[n]));
246         else
247             memset(&t[n2], 0, sizeof(*t) * 8);
248
249         bn_mul_comba4(r, a, b);
250         bn_mul_comba4(&(r[n2]), &(a[n]), &(b[n]));
251     } else if (n == 8 && dna == 0 && dnb == 0) { /* XXX: bn_mul_comba8 could
252                                                   * take extra args to do
253                                                   * this well */
254         if (!zero)
255             bn_mul_comba8(&(t[n2]), t, &(t[n]));
256         else
257             memset(&t[n2], 0, sizeof(*t) * 16);
258
259         bn_mul_comba8(r, a, b);
260         bn_mul_comba8(&(r[n2]), &(a[n]), &(b[n]));
261     } else
262 # endif                         /* BN_MUL_COMBA */
263     {
264         p = &(t[n2 * 2]);
265         if (!zero)
266             bn_mul_recursive(&(t[n2]), t, &(t[n]), n, 0, 0, p);
267         else
268             memset(&t[n2], 0, sizeof(*t) * n2);
269         bn_mul_recursive(r, a, b, n, 0, 0, p);
270         bn_mul_recursive(&(r[n2]), &(a[n]), &(b[n]), n, dna, dnb, p);
271     }
272
273     /*-
274      * t[32] holds (a[0]-a[1])*(b[1]-b[0]), c1 is the sign
275      * r[10] holds (a[0]*b[0])
276      * r[32] holds (b[1]*b[1])
277      */
278
279     c1 = (int)(bn_add_words(t, r, &(r[n2]), n2));
280
281     if (neg) {                  /* if t[32] is negative */
282         c1 -= (int)(bn_sub_words(&(t[n2]), t, &(t[n2]), n2));
283     } else {
284         /* Might have a carry */
285         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), t, n2));
286     }
287
288     /*-
289      * t[32] holds (a[0]-a[1])*(b[1]-b[0])+(a[0]*b[0])+(a[1]*b[1])
290      * r[10] holds (a[0]*b[0])
291      * r[32] holds (b[1]*b[1])
292      * c1 holds the carry bits
293      */
294     c1 += (int)(bn_add_words(&(r[n]), &(r[n]), &(t[n2]), n2));
295     if (c1) {
296         p = &(r[n + n2]);
297         lo = *p;
298         ln = (lo + c1) & BN_MASK2;
299         *p = ln;
300
301         /*
302          * The overflow will stop before we over write words we should not
303          * overwrite
304          */
305         if (ln < (BN_ULONG)c1) {
306             do {
307                 p++;
308                 lo = *p;
309                 ln = (lo + 1) & BN_MASK2;
310                 *p = ln;
311             } while (ln == 0);
312         }
313     }
314 }
315
316 /*
317  * n+tn is the word length t needs to be n*4 is size, as does r
318  */
319 /* tnX may not be negative but less than n */
320 void bn_mul_part_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n,
321                            int tna, int tnb, BN_ULONG *t)
322 {
323     int i, j, n2 = n * 2;
324     int c1, c2, neg;
325     BN_ULONG ln, lo, *p;
326
327     if (n < 8) {
328         bn_mul_normal(r, a, n + tna, b, n + tnb);
329         return;
330     }
331
332     /* r=(a[0]-a[1])*(b[1]-b[0]) */
333     c1 = bn_cmp_part_words(a, &(a[n]), tna, n - tna);
334     c2 = bn_cmp_part_words(&(b[n]), b, tnb, tnb - n);
335     neg = 0;
336     switch (c1 * 3 + c2) {
337     case -4:
338         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
339         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
340         break;
341     case -3:
342     case -2:
343         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
344         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n); /* + */
345         neg = 1;
346         break;
347     case -1:
348     case 0:
349     case 1:
350     case 2:
351         bn_sub_part_words(t, a, &(a[n]), tna, n - tna); /* + */
352         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
353         neg = 1;
354         break;
355     case 3:
356     case 4:
357         bn_sub_part_words(t, a, &(a[n]), tna, n - tna);
358         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n);
359         break;
360     }
361     /*
362      * The zero case isn't yet implemented here. The speedup would probably
363      * be negligible.
364      */
365 # if 0
366     if (n == 4) {
367         bn_mul_comba4(&(t[n2]), t, &(t[n]));
368         bn_mul_comba4(r, a, b);
369         bn_mul_normal(&(r[n2]), &(a[n]), tn, &(b[n]), tn);
370         memset(&r[n2 + tn * 2], 0, sizeof(*r) * (n2 - tn * 2));
371     } else
372 # endif
373     if (n == 8) {
374         bn_mul_comba8(&(t[n2]), t, &(t[n]));
375         bn_mul_comba8(r, a, b);
376         bn_mul_normal(&(r[n2]), &(a[n]), tna, &(b[n]), tnb);
377         memset(&r[n2 + tna + tnb], 0, sizeof(*r) * (n2 - tna - tnb));
378     } else {
379         p = &(t[n2 * 2]);
380         bn_mul_recursive(&(t[n2]), t, &(t[n]), n, 0, 0, p);
381         bn_mul_recursive(r, a, b, n, 0, 0, p);
382         i = n / 2;
383         /*
384          * If there is only a bottom half to the number, just do it
385          */
386         if (tna > tnb)
387             j = tna - i;
388         else
389             j = tnb - i;
390         if (j == 0) {
391             bn_mul_recursive(&(r[n2]), &(a[n]), &(b[n]),
392                              i, tna - i, tnb - i, p);
393             memset(&r[n2 + i * 2], 0, sizeof(*r) * (n2 - i * 2));
394         } else if (j > 0) {     /* eg, n == 16, i == 8 and tn == 11 */
395             bn_mul_part_recursive(&(r[n2]), &(a[n]), &(b[n]),
396                                   i, tna - i, tnb - i, p);
397             memset(&(r[n2 + tna + tnb]), 0,
398                    sizeof(BN_ULONG) * (n2 - tna - tnb));
399         } else {                /* (j < 0) eg, n == 16, i == 8 and tn == 5 */
400
401             memset(&r[n2], 0, sizeof(*r) * n2);
402             if (tna < BN_MUL_RECURSIVE_SIZE_NORMAL
403                 && tnb < BN_MUL_RECURSIVE_SIZE_NORMAL) {
404                 bn_mul_normal(&(r[n2]), &(a[n]), tna, &(b[n]), tnb);
405             } else {
406                 for (;;) {
407                     i /= 2;
408                     /*
409                      * these simplified conditions work exclusively because
410                      * difference between tna and tnb is 1 or 0
411                      */
412                     if (i < tna || i < tnb) {
413                         bn_mul_part_recursive(&(r[n2]),
414                                               &(a[n]), &(b[n]),
415                                               i, tna - i, tnb - i, p);
416                         break;
417                     } else if (i == tna || i == tnb) {
418                         bn_mul_recursive(&(r[n2]),
419                                          &(a[n]), &(b[n]),
420                                          i, tna - i, tnb - i, p);
421                         break;
422                     }
423                 }
424             }
425         }
426     }
427
428     /*-
429      * t[32] holds (a[0]-a[1])*(b[1]-b[0]), c1 is the sign
430      * r[10] holds (a[0]*b[0])
431      * r[32] holds (b[1]*b[1])
432      */
433
434     c1 = (int)(bn_add_words(t, r, &(r[n2]), n2));
435
436     if (neg) {                  /* if t[32] is negative */
437         c1 -= (int)(bn_sub_words(&(t[n2]), t, &(t[n2]), n2));
438     } else {
439         /* Might have a carry */
440         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), t, n2));
441     }
442
443     /*-
444      * t[32] holds (a[0]-a[1])*(b[1]-b[0])+(a[0]*b[0])+(a[1]*b[1])
445      * r[10] holds (a[0]*b[0])
446      * r[32] holds (b[1]*b[1])
447      * c1 holds the carry bits
448      */
449     c1 += (int)(bn_add_words(&(r[n]), &(r[n]), &(t[n2]), n2));
450     if (c1) {
451         p = &(r[n + n2]);
452         lo = *p;
453         ln = (lo + c1) & BN_MASK2;
454         *p = ln;
455
456         /*
457          * The overflow will stop before we over write words we should not
458          * overwrite
459          */
460         if (ln < (BN_ULONG)c1) {
461             do {
462                 p++;
463                 lo = *p;
464                 ln = (lo + 1) & BN_MASK2;
465                 *p = ln;
466             } while (ln == 0);
467         }
468     }
469 }
470
471 /*-
472  * a and b must be the same size, which is n2.
473  * r needs to be n2 words and t needs to be n2*2
474  */
475 void bn_mul_low_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n2,
476                           BN_ULONG *t)
477 {
478     int n = n2 / 2;
479
480     bn_mul_recursive(r, a, b, n, 0, 0, &(t[0]));
481     if (n >= BN_MUL_LOW_RECURSIVE_SIZE_NORMAL) {
482         bn_mul_low_recursive(&(t[0]), &(a[0]), &(b[n]), n, &(t[n2]));
483         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
484         bn_mul_low_recursive(&(t[0]), &(a[n]), &(b[0]), n, &(t[n2]));
485         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
486     } else {
487         bn_mul_low_normal(&(t[0]), &(a[0]), &(b[n]), n);
488         bn_mul_low_normal(&(t[n]), &(a[n]), &(b[0]), n);
489         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
490         bn_add_words(&(r[n]), &(r[n]), &(t[n]), n);
491     }
492 }
493 #endif                          /* BN_RECURSION */
494
495 int BN_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
496 {
497     int ret = 0;
498     int top, al, bl;
499     BIGNUM *rr;
500 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
501     int i;
502 #endif
503 #ifdef BN_RECURSION
504     BIGNUM *t = NULL;
505     int j = 0, k;
506 #endif
507
508     bn_check_top(a);
509     bn_check_top(b);
510     bn_check_top(r);
511
512     al = a->top;
513     bl = b->top;
514
515     if ((al == 0) || (bl == 0)) {
516         BN_zero(r);
517         return (1);
518     }
519     top = al + bl;
520
521     BN_CTX_start(ctx);
522     if ((r == a) || (r == b)) {
523         if ((rr = BN_CTX_get(ctx)) == NULL)
524             goto err;
525     } else
526         rr = r;
527
528 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
529     i = al - bl;
530 #endif
531 #ifdef BN_MUL_COMBA
532     if (i == 0) {
533 # if 0
534         if (al == 4) {
535             if (bn_wexpand(rr, 8) == NULL)
536                 goto err;
537             rr->top = 8;
538             bn_mul_comba4(rr->d, a->d, b->d);
539             goto end;
540         }
541 # endif
542         if (al == 8) {
543             if (bn_wexpand(rr, 16) == NULL)
544                 goto err;
545             rr->top = 16;
546             bn_mul_comba8(rr->d, a->d, b->d);
547             goto end;
548         }
549     }
550 #endif                          /* BN_MUL_COMBA */
551 #ifdef BN_RECURSION
552     if ((al >= BN_MULL_SIZE_NORMAL) && (bl >= BN_MULL_SIZE_NORMAL)) {
553         if (i >= -1 && i <= 1) {
554             /*
555              * Find out the power of two lower or equal to the longest of the
556              * two numbers
557              */
558             if (i >= 0) {
559                 j = BN_num_bits_word((BN_ULONG)al);
560             }
561             if (i == -1) {
562                 j = BN_num_bits_word((BN_ULONG)bl);
563             }
564             j = 1 << (j - 1);
565             assert(j <= al || j <= bl);
566             k = j + j;
567             t = BN_CTX_get(ctx);
568             if (t == NULL)
569                 goto err;
570             if (al > j || bl > j) {
571                 if (bn_wexpand(t, k * 4) == NULL)
572                     goto err;
573                 if (bn_wexpand(rr, k * 4) == NULL)
574                     goto err;
575                 bn_mul_part_recursive(rr->d, a->d, b->d,
576                                       j, al - j, bl - j, t->d);
577             } else {            /* al <= j || bl <= j */
578
579                 if (bn_wexpand(t, k * 2) == NULL)
580                     goto err;
581                 if (bn_wexpand(rr, k * 2) == NULL)
582                     goto err;
583                 bn_mul_recursive(rr->d, a->d, b->d, j, al - j, bl - j, t->d);
584             }
585             rr->top = top;
586             goto end;
587         }
588 # if 0
589         if (i == 1 && !BN_get_flags(b, BN_FLG_STATIC_DATA)) {
590             BIGNUM *tmp_bn = (BIGNUM *)b;
591             if (bn_wexpand(tmp_bn, al) == NULL)
592                 goto err;
593             tmp_bn->d[bl] = 0;
594             bl++;
595             i--;
596         } else if (i == -1 && !BN_get_flags(a, BN_FLG_STATIC_DATA)) {
597             BIGNUM *tmp_bn = (BIGNUM *)a;
598             if (bn_wexpand(tmp_bn, bl) == NULL)
599                 goto err;
600             tmp_bn->d[al] = 0;
601             al++;
602             i++;
603         }
604         if (i == 0) {
605             /* symmetric and > 4 */
606             /* 16 or larger */
607             j = BN_num_bits_word((BN_ULONG)al);
608             j = 1 << (j - 1);
609             k = j + j;
610             t = BN_CTX_get(ctx);
611             if (al == j) {      /* exact multiple */
612                 if (bn_wexpand(t, k * 2) == NULL)
613                     goto err;
614                 if (bn_wexpand(rr, k * 2) == NULL)
615                     goto err;
616                 bn_mul_recursive(rr->d, a->d, b->d, al, t->d);
617             } else {
618                 if (bn_wexpand(t, k * 4) == NULL)
619                     goto err;
620                 if (bn_wexpand(rr, k * 4) == NULL)
621                     goto err;
622                 bn_mul_part_recursive(rr->d, a->d, b->d, al - j, j, t->d);
623             }
624             rr->top = top;
625             goto end;
626         }
627 # endif
628     }
629 #endif                          /* BN_RECURSION */
630     if (bn_wexpand(rr, top) == NULL)
631         goto err;
632     rr->top = top;
633     bn_mul_normal(rr->d, a->d, al, b->d, bl);
634
635 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
636  end:
637 #endif
638     rr->neg = a->neg ^ b->neg;
639     bn_correct_top(rr);
640     if (r != rr && BN_copy(r, rr) == NULL)
641         goto err;
642
643     ret = 1;
644  err:
645     bn_check_top(r);
646     BN_CTX_end(ctx);
647     return (ret);
648 }
649
650 void bn_mul_normal(BN_ULONG *r, BN_ULONG *a, int na, BN_ULONG *b, int nb)
651 {
652     BN_ULONG *rr;
653
654     if (na < nb) {
655         int itmp;
656         BN_ULONG *ltmp;
657
658         itmp = na;
659         na = nb;
660         nb = itmp;
661         ltmp = a;
662         a = b;
663         b = ltmp;
664
665     }
666     rr = &(r[na]);
667     if (nb <= 0) {
668         (void)bn_mul_words(r, a, na, 0);
669         return;
670     } else
671         rr[0] = bn_mul_words(r, a, na, b[0]);
672
673     for (;;) {
674         if (--nb <= 0)
675             return;
676         rr[1] = bn_mul_add_words(&(r[1]), a, na, b[1]);
677         if (--nb <= 0)
678             return;
679         rr[2] = bn_mul_add_words(&(r[2]), a, na, b[2]);
680         if (--nb <= 0)
681             return;
682         rr[3] = bn_mul_add_words(&(r[3]), a, na, b[3]);
683         if (--nb <= 0)
684             return;
685         rr[4] = bn_mul_add_words(&(r[4]), a, na, b[4]);
686         rr += 4;
687         r += 4;
688         b += 4;
689     }
690 }
691
692 void bn_mul_low_normal(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n)
693 {
694     bn_mul_words(r, a, n, b[0]);
695
696     for (;;) {
697         if (--n <= 0)
698             return;
699         bn_mul_add_words(&(r[1]), a, n, b[1]);
700         if (--n <= 0)
701             return;
702         bn_mul_add_words(&(r[2]), a, n, b[2]);
703         if (--n <= 0)
704             return;
705         bn_mul_add_words(&(r[3]), a, n, b[3]);
706         if (--n <= 0)
707             return;
708         bn_mul_add_words(&(r[4]), a, n, b[4]);
709         r += 4;
710         b += 4;
711     }
712 }