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