Import of old SSLeay release: SSLeay 0.8.1b
[openssl.git] / crypto / bn / stuff / wei_mulw.c
1 /* crypto/bn/wei_mulw.c */
2
3 #include <stdio.h>
4 #include "cryptlib.h"
5 #include "bn.h"
6 #include "bn_lcl.h"
7
8 BN_ULONG bn_add_word(BN_ULONG *a,BN_ULONG c,int num);
9 BN_ULONG bn_add_words(BN_ULONG *ret,BN_ULONG *a,BN_ULONG *b,int num);
10 BN_ULONG bn_sub_words(BN_ULONG *ret,BN_ULONG *a,BN_ULONG *b,int num);
11
12 void BN_mul_4words(BN_ULONG *ret,BN_ULONG a0,BN_ULONG a1,
13         BN_ULONG b0,BN_ULONG b1);
14
15 void pr(a,n,s)
16 BN_ULONG *a;
17 int n;
18         {
19         while (n--)
20                 fprintf(stdout,"%02X",a[n]);
21         fprintf(stdout,"%s",s);
22         }
23
24
25 BN_ULONG bn_add_word(a,w,num)
26 BN_ULONG *a;
27 BN_ULONG w;
28 int num;
29         {
30         BN_ULONG t;
31
32 #ifdef DEBUG
33 { BN_ULONG *aa=a; int i; for (i=num; i>0; i--) fprintf(stdout,"%02X",aa[i-1]);
34 fprintf(stdout," + %X - ",w); i=num;
35 #endif
36         
37 loop:
38         t= *a;
39         t=(t+w)&BN_MASK2;
40         *(a++)=t;
41         w=(t < w);
42         if (w && --num) goto loop;
43
44 #ifdef DEBUG
45 for (; i>0; i--) fprintf(stdout,"%02X",aa[i-1]);
46 fprintf(stdout,"\n");
47 }
48 #endif
49
50         return(w);
51         }
52
53 BN_ULONG bn_add_words(r,a,b,num)
54 BN_ULONG *r;
55 BN_ULONG *a;
56 BN_ULONG *b;
57 int num;
58         {
59 #if defined(BN_LLONG)
60         BN_ULLONG t;
61         BN_ULONG c=0;
62         int i;
63
64         if (num&1) abort();
65
66         for (i=0; i<num; i+=2)
67                 {
68                 t=(BN_ULLONG)a[i]+b[i]+c;
69                 r[i+0]=L(t);
70                 t=(BN_ULLONG) H(t)+a[i+1]+b[i+1];
71                 r[i+1]=L(t);
72                 c=H(t);
73                 }
74         return(c);
75 #else
76         BN_ULONG c=0,t1,t2;
77
78         for ( ; num; num--)
79                 {
80                 t1= *(a++);
81                 t2= *(b++);
82
83                 if (c)
84                         {
85                         c=(t2 >= ((~t1)&BN_MASK2));
86                         (*r++)=(t1+t2+1)&BN_MASK2;
87                         }
88                 else
89                         {
90                         t2=(t1+t2)&BN_MASK2;
91                         c=(t2 < t1);
92                         (*r++)=t2;
93                         }
94                 }
95         return(c);
96 #endif
97         }
98
99 BN_ULONG bn_sub_words(r,a,b,num)
100 BN_ULONG *r;
101 BN_ULONG *a;
102 BN_ULONG *b;
103 int num;
104         {
105 #if defined(BN_LLONG)
106         BN_ULLONG t;
107         BN_ULONG c=0;
108         int i;
109
110         if (num&1) abort();
111
112         for (i=0; i<num; i+=2)
113                 {
114                 t=(BN_ULLONG)a[i]-b[i]-c;
115                 r[i+0]=L(t);
116                 t=(BN_ULLONG)a[i+1]-b[i+1]-(0-H(t))&BN_MASK2;
117                 r[i+1]=L(t);
118                 c=H(t);
119                 }
120         return(c);
121 #else
122         BN_ULONG c=0,t1,t2;
123
124         for ( ; num; num--)
125                 {
126                 t1= *(a++);
127                 t2= *(b++);
128
129                 if (c)
130                         {
131                         c=(t1 <= t2);
132                         t1=(t1-t2-1);
133                         }
134                 else
135                         {
136                         c=(t1 < t2);
137                         t1=(t1-t2);
138                         }
139                 (*r++)=t1&BN_MASK2;
140                 }
141         return(c);
142 #endif
143         }
144
145
146 /* ret[3,2,1,0] = a1,a0 * b1,b0 */
147 void BN_mul_4words(ret,a0,a1,b0,b1)
148 BN_ULONG *ret;
149 BN_ULONG a0,a1,b0,b1;
150         {
151         BN_ULONG s,u;
152         BN_ULLONG fix,a0b0,a1b1,tmp;
153
154         if (a1 >= a0)
155                 {
156                 s=(a1-a0);
157                 u=(b0-b1);
158                 fix=(BN_ULLONG)s*u;
159                 if (b0 >= b1) s=0;
160                 }
161         else
162                 {
163                 BN_ULONG u;
164
165                 if (b0 > b1)
166                         {
167                         s=(b0-b1);
168                         u=(a1-a0);
169                         fix=(BN_ULLONG)s*u;
170                         }
171                 else
172                         {
173                         u=(a0-a1);
174                         s=(b1-b0);
175                         fix=(BN_ULLONG)s*u;
176                         s=0;
177                         }
178                 }
179         
180         a0b0=(BN_ULLONG)a0*b0;
181         ret[0]=L(a0b0);
182
183         a1b1=(BN_ULLONG)a1*b1;
184         tmp=(BN_ULLONG) H(a0b0) + L(a0b0) + L(fix) + L(a1b1);
185         ret[1]=L(tmp);
186
187         tmp=(BN_ULLONG) a1b1 + H(tmp) + H(a0b0) + H(fix) + H(a1b1) - s;
188         ret[2]=L(tmp);
189         ret[3]=H(tmp);
190         }
191
192 /* ret[3,2,1,0] += a1,a0 * b1,b0 */
193 BN_ULONG BN_mul_add_4words(ret,a0,a1,b0,b1)
194 BN_ULONG *ret;
195 BN_ULONG a0,a1,b0,b1;
196         {
197         BN_ULONG s,u;
198         BN_ULLONG fix,a0b0,a1b1,tmp;
199
200 #ifdef DEBUG
201 fprintf(stdout,"%02X%02X%02X%02X",ret[3],ret[2],ret[1],ret[0]);
202 fprintf(stdout," + ( %02X%02X * %02X%02X ) - ",a1,a0,b1,b0);
203 #endif
204         if (a1 >= a0)
205                 {
206                 s=(a1-a0);
207                 u=(b0-b1);
208                 fix=(BN_ULLONG)s*u;
209                 if (b0 >= b1) s=0;
210                 }
211         else
212                 {
213                 if (b0 > b1)
214                         {
215                         s=(b0-b1);
216                         u=(a1-a0);
217                         fix=(BN_ULLONG)s*u;
218                         }
219                 else
220                         {
221                         u=(a0-a1);
222                         s=(b1-b0);
223                         fix=(BN_ULLONG)s*u;
224                         s=0;
225                         }
226                 }
227         
228         a0b0=(BN_ULLONG)a0*b0;
229         tmp=a0b0+ret[0];
230         ret[0]=L(tmp);
231
232         a1b1=(BN_ULLONG)a1*b1;
233         tmp=(BN_ULLONG) H(tmp) + L(a0b0) + L(fix) + L(a1b1) + ret[1];
234         ret[1]=L(tmp);
235
236         tmp=(BN_ULLONG) H(tmp) + L(a1b1) + H(a0b0) +
237                 H(fix) + H(a1b1) -s + ret[2];
238         ret[2]=L(tmp);
239
240         tmp=(BN_ULLONG) H(tmp) + H(a1b1) + ret[3];
241         ret[3]=L(tmp);
242 #ifdef DEBUG
243 fprintf(stdout,"%02X%02X%02X%02X%02X\n",H(tmp),ret[3],ret[2],ret[1],ret[0]);
244 #endif
245         return(H(tmp));
246         }
247
248 /* ret[3,2,1,0] += a1,a0 * a1,a0 */
249 void BN_sqr_4words(ret,a0,a1)
250 BN_ULONG *ret;
251 BN_ULONG a0,a1;
252         {
253         BN_ULONG s,u;
254         BN_ULLONG tmp,tmp2;
255
256         tmp=(BN_ULLONG)a0*a0;
257         ret[0]=L(tmp);
258
259         tmp2=(BN_ULLONG)a0*a1;
260         tmp=(BN_ULLONG)H(tmp)+L(tmp2)*2;
261         ret[1]=L(tmp);
262
263         tmp=(BN_ULLONG)a1*a1+H(tmp)+H(tmp2)*2;
264         ret[2]=L(tmp);
265         ret[3]=L(tmp);
266         }
267
268 #define N0      (0)
269 #define N1      (half)
270 #define N2      (num)
271 #define N3      (num+half)
272
273 #define word_cmp(r,a,b,num) \
274         { \
275         int n=num; \
276 \
277         (r)=0; \
278         while (n--) \
279                 { \
280                 if ((a)[(n)] > (b)[(n)]) \
281                         { (r)=1; break; } \
282                 else if ((a)[(n)] < (b)[(n)]) \
283                         { (r)= -1; break; } \
284                 } \
285         }
286
287
288 /* (a->top == b->top) && (a->top >= 2) && !(a->top & 1) */
289 void bn_recursize_mul(r,t,a,b,num)
290 BN_ULONG *r,*t,*a,*b;
291 int num;
292         {
293         if ((num < 2) || (num&1))
294                 abort();
295
296 /* fprintf(stderr,"num=%d half=%d\n",num,num/2);*/
297         if (num == 2)
298                 BN_mul_4words(r,a[0],a[1],b[0],b[1]);
299         else if (num == 4)
300                 {
301                 BN_ULONG c,tmp;
302
303                 BN_mul_4words(&(r[0]),a[0],a[1],b[0],b[1]);
304                 BN_mul_4words(&(r[4]),a[2],a[3],b[2],b[3]);
305
306                 c =BN_mul_add_4words(&(r[2]),a[0],a[1],b[2],b[3]);
307                 c+=BN_mul_add_4words(&(r[2]),a[2],a[3],b[0],b[1]);
308
309                 bn_add_word(&(r[6]),c,2);
310                 }
311         else
312                 {
313                 int half=num/2;
314                 int carry,cmp_a,cmp_b;
315
316                 word_cmp(cmp_a,&(a[0]),&(a[half]),half);
317                 word_cmp(cmp_b,&(b[0]),&(b[half]),half);
318
319                 switch (cmp_a*2+cmp_a+cmp_b)
320                         {
321                 case -4:
322                         bn_sub_words(&(t[N0]),&(a[N1]),&(a[N0]),half);
323                         bn_sub_words(&(t[N1]),&(b[N0]),&(b[N1]),half);
324                         bn_recursize_mul(&(r[N1]),&(t[N2]),
325                                 &(t[N0]),&(t[N1]),half);
326                         bn_sub_words(&(r[N2]),&(r[N2]),&(t[N0]),half);
327                         carry= -1;
328                         break;
329                 case -2:
330                         bn_sub_words(&(t[N0]),&(a[N1]),&(a[N0]),half);
331                         bn_sub_words(&(t[N1]),&(b[N0]),&(b[N1]),half);
332                         bn_recursize_mul(&(r[N1]),&(t[N2]),
333                                 &(t[N0]),&(t[N1]),half);
334                         carry=0;
335                         break;
336                 case 2:
337                         bn_sub_words(&(t[N0]),&(a[N0]),&(a[N1]),half);
338                         bn_sub_words(&(t[N1]),&(b[N1]),&(b[N0]),half);
339                         bn_recursize_mul(&(r[N1]),&(t[N2]),
340                                 &(t[N0]),&(t[N1]),half);
341                         carry=0;
342                         break;
343                 case 4:
344                         bn_sub_words(&(t[N0]),&(a[N1]),&(a[N0]),half);
345                         bn_sub_words(&(t[N1]),&(b[N0]),&(b[N1]),half);
346                         bn_recursize_mul(&(r[N1]),&(t[N2]),
347                                 &(t[N0]),&(t[N1]),half);
348                         bn_sub_words(&(r[N2]),&(r[N2]),&(t[N1]),half);
349                         carry= -1;
350                         break;
351                 default:
352                         memset(&(r[N1]),0,sizeof(BN_ULONG)*num);
353                         break;
354                         }
355                 
356                 bn_recursize_mul(&(t[N0]),&(t[N2]),&(a[N0]),&(b[N0]),half);
357 #ifdef DEBUG
358         pr(a,half," * ");
359         pr(b,half," - ");
360         pr(t,num," - 0\n");
361 #endif
362                 memcpy(&(r[N0]),&(t[N0]),half*sizeof(BN_ULONG));
363                 if (bn_add_words(&(r[N1]),&(r[N1]),&(t[N1]),half))
364                         { bn_add_word(&(t[N1]),1,half); }
365
366                 carry+=bn_add_words(&(r[N1]),&(r[N1]),&(t[N0]),num);
367
368                 bn_recursize_mul(&(t[N0]),&(t[N2]),&(a[N1]),&(b[N1]),half);
369
370                 carry+=bn_add_words(&(r[N1]),&(r[N1]),&(t[N0]),num);
371                 carry+=bn_add_words(&(r[N2]),&(r[N2]),&(t[N0]),half);
372                 memcpy(&(r[N3]),&(t[N1]),half*sizeof(BN_ULONG));
373
374                 bn_add_word(&(r[N3]),carry,half);
375                 }
376         }
377
378 main()
379         {
380         BIGNUM *a,*b,*r,*t;
381         int i,j;
382
383         a=BN_new();
384         b=BN_new();
385         r=BN_new();
386         t=BN_new();
387
388 #define BITS 1024
389         bn_expand(r,BITS*2);
390         bn_expand(t,BITS*2);
391         fprintf(stdout,"obase=16\n");
392         fprintf(stdout,"ibase=16\n");
393         for (i=0; i<10; i++)
394                 {
395                 BN_rand(a,BITS,0,0);
396                 BN_rand(b,BITS,0,0);
397                 r->top=(BITS*2)/BN_BITS2;
398                 memset(r->d,0,sizeof(r->top)*sizeof(BN_ULONG));
399                 memset(t->d,0,sizeof(r->top)*sizeof(BN_ULONG));
400                 for (j=0; j<1000; j++)
401                         {
402
403 /*                      BN_mul(r,a,b); /**/
404                         bn_recursize_mul(r->d,t->d,a->d,b->d,a->top); /**/
405                         }
406                 BN_print(stdout,a); fprintf(stdout," * ");
407                 BN_print(stdout,b); fprintf(stdout," - ");
408                 BN_print(stdout,r); fprintf(stdout,"\n");
409                 }
410         }