Import of old SSLeay release: SSLeay 0.8.1b
[openssl.git] / crypto / bn / stuff / div.c
1 /* crypto/bn/div.c */
2
3 #include <stdio.h>
4 #include "cryptlib.h"
5 #include "bn.h"
6
7 BN_ULONG bn_div_2word();
8
9 int BN_div2(dv, rm, num, div,ctx)
10 BIGNUM *dv;
11 BIGNUM *rm;
12 BIGNUM *num;
13 BIGNUM *div;
14 BN_CTX *ctx;
15         {
16         int norm_shift,i,j,nm,nd,loop;
17         BIGNUM *tmp,wnum,*snum,*sdiv,*res;
18         BN_ULONG *resp,*wnump;
19         BN_ULONG d0,d1;
20         int num_n,div_n;
21
22 #ifdef DEBUG
23 BN_print(stdout,num); printf(" number\n");
24 BN_print(stdout,div); printf(" divisor\n");
25 #endif
26         if (BN_is_zero(num))
27                 {
28                 BNerr(BN_F_BN_DIV,BN_R_DIV_BY_ZERO);
29                 return(0);
30                 }
31
32         if (BN_cmp(num,div) < 0)
33                 {
34                 if (rm != NULL)
35                         { if (BN_copy(rm,num) == NULL) return(0); }
36                 if (dv != NULL) BN_zero(dv);
37                 return(1);
38                 }
39
40         tmp=ctx->bn[ctx->tos]; 
41         snum=ctx->bn[ctx->tos+1];
42         sdiv=ctx->bn[ctx->tos+2];
43         if (dv == NULL)
44                 res=ctx->bn[ctx->tos+3];
45         else    res=dv;
46
47         /* First we normalise the numbers */
48         norm_shift=BN_BITS2-((BN_num_bits(div))%BN_BITS2);
49         BN_lshift(sdiv,div,norm_shift);
50         norm_shift+=BN_BITS2;
51         BN_lshift(snum,num,norm_shift);
52         div_n=sdiv->top;
53         num_n=snum->top;
54         loop=num_n-div_n;
55 #ifdef DEBUG
56 BN_print(stdout,snum); printf(" shifted num, forget last word\n");
57 BN_print(stdout,sdiv); printf(" shifted div\n");
58 #endif
59
60         /* Lets setup a 'win'dow into snum
61          * This is the part that corresponds to the current
62          * 'area' being divided */
63         wnum.d=  &(snum->d[loop]);
64         wnum.top= div_n;
65         wnum.max= snum->max; /* a bit of a lie */
66         wnum.neg= 0;
67
68         /* Get the top 2 words of sdiv */
69         i=sdiv->top;
70         d0=sdiv->d[div_n-1];
71         d1=sdiv->d[div_n-2];
72
73         /* pointer to the 'top' of snum */
74         wnump= &(snum->d[num_n-1]);
75
76         /* Setup to 'res' */
77         res->neg=0;
78         res->top=loop;
79         resp= &(res->d[loop-1]);
80         bn_expand(res,(loop+1)*BN_BITS2);
81
82         /* space for temp */
83         bn_expand(tmp,(div_n+1)*BN_BITS2);
84
85 #ifdef DEBUG
86 printf("wnum="); BN_print(stdout,&wnum); printf(" initial sub check\n");
87 printf("div ="); BN_print(stdout,sdiv); printf(" loop=%d\n",loop);
88 #endif
89         if (BN_cmp(&wnum,sdiv) >= 0)
90                 {
91                 BN_sub(&wnum,&wnum,sdiv);
92                 *resp=1;
93                 res->d[res->top-1]=1;
94                 }
95         else
96                 res->top--;
97         resp--;
98 #ifdef DEBUG
99 BN_print(stdout,res); printf(" initial result\n");
100 BN_print(stdout,&wnum); printf(" wnum\n");
101 #endif
102
103         for (i=0; i<loop-1; i++)
104                 {
105                 BN_ULONG q,n0;
106                 BN_ULLONG t1,t2,t3;
107                 BN_ULONG l0;
108
109                 wnum.d--;
110                 wnum.top++;
111
112 #ifdef DEBUG
113 BN_print(stderr,&wnum); printf(" to divide\n");
114 #endif
115
116                 q=0;
117                 n0=wnump[0];
118                 t1=((BN_ULLONG)n0<<BN_BITS2)|wnump[-1];
119                 if (n0 == d0)
120                         q=BN_MASK2;
121                 else
122                         {
123                         t2=(t1/d0);
124                         q=(t2&BN_MASK2);
125 #ifdef DEBUG
126 printf("t1=%08X / d0=%08X = %X (%X)\n",t1,d0,q,t2);
127 #endif
128                         }
129                 for (;;)
130                         {
131                         t2=(BN_ULLONG)d1*q;
132                         t3=t1-(BN_ULLONG)q*d0;
133 #ifdef DEBUG
134 printf("d1*q= %X    n01-q*d0 = %X\n",t2,t3);
135 #endif
136                         if ((t3>>BN_BITS2) ||
137                                 (t2 <= ((t3<<BN_BITS2)+wnump[-2])))
138                                 break;
139 #ifdef DEBUG
140 printf("reduce q\n");
141 #endif
142                         q--;
143                         }
144                 l0=bn_mul_word(tmp->d,sdiv->d,div_n,q);
145                 if (l0)
146                         tmp->d[div_n]=l0;
147                 else
148                         tmp->d[div_n]=0;
149                 for (j=div_n+1; j>0; j--)
150                         if (tmp->d[j-1]) break;
151                 tmp->top=j;
152
153 #ifdef DEBUG
154 printf("q=%08X\n",q);
155 BN_print(stdout,&wnum); printf(" number\n");
156 BN_print(stdout,tmp); printf(" subtract\n");
157
158 BN_print(stdout,snum); printf(" shifted number before\n");
159 BN_print(stdout,&wnum); printf(" wnum before\n");
160 #endif
161                 j=wnum.top;
162                 BN_sub(&wnum,&wnum,tmp);
163                 snum->top=snum->top+wnum.top-j;
164
165 #ifdef DEBUG
166 BN_print(stdout,&wnum); printf(" wnum after\n");
167 BN_print(stdout,snum); printf(" shifted number after\n");
168 #endif
169
170                 if (wnum.neg)
171                         {
172                         q--;
173                         j=wnum.top;
174                         BN_add(&wnum,&wnum,sdiv);
175                         snum->top+=wnum.top-j;
176                         fprintf(stderr,"addback\n");
177 #ifdef DEBUG
178 BN_print(stdout,snum); printf("after addback************************:\n");
179 #endif
180                         }
181                 *(resp--)=q;
182 #ifdef DEBUG
183 BN_print(stdout,res); printf(" result\n");
184 #endif
185                 wnump--;
186                 }
187         if (rm != NULL)
188                 BN_rshift(rm,snum,norm_shift);
189         return(1);
190         }
191
192 main()
193         {
194         BIGNUM *a,*b,*c,*d;
195         BIGNUM *cc,*dd;
196         BN_CTX *ctx;
197         int i,x;
198
199         a=BN_new();
200         b=BN_new();
201         c=BN_new();
202         d=BN_new();
203         cc=BN_new();
204         dd=BN_new();
205         ctx=BN_CTX_new();
206
207 for (i=0; i<10240; i++)
208         {
209         BN_rand(a,80,0,0);
210         BN_rand(b,60,0,0);
211         
212         BN_div2(d,c,a,b,ctx);
213         BN_div(dd,cc,a,b,ctx);
214         if ((BN_cmp(d,dd) != 0) || (BN_cmp(c,cc) != 0))
215                 {
216                 BN_print(stderr,a); fprintf(stderr," / ");
217                 BN_print(stderr,b); fprintf(stderr," d=");
218                 BN_print(stderr,d); fprintf(stderr," r= ");
219                 BN_print(stderr,c); fprintf(stderr,"\nd=");
220                 BN_print(stderr,dd); fprintf(stderr," r= ");
221                 BN_print(stderr,cc); fprintf(stderr,"\n");
222                 }
223
224         }
225
226 #ifdef undef
227 /*
228         BN_rand(a,600,0,0);
229         BN_rand(b,400,0,0);
230         for (i=0; i<2000000; i++)
231                 {
232                 BN_div2(d,c,a,b,ctx);
233                 }
234 */
235 /*      for (i=0;;) */
236 /*      for (i=0; i<0xffffffff; i++)
237                 {
238                 BN_ULONG rr,r,a,b,c;
239                 BN_ULLONG l;
240
241                 a=rand()&BN_MASK2;
242                 b=rand()&BN_MASK2;
243                 for (;;)
244                         {
245                         c=rand()&BN_MASK2;
246                         if (c) break;
247                         }
248 /*              for (x=1; x<256*256; x++) */
249                         {
250                         c=x;
251                         a=i>>8;
252                         b=i&0xff;
253                         a&= ~(0xFFFFFF<<(BN_num_bits_word(c)));
254
255                         r=bn_div_2word(a,b,c);
256
257                         rr=(BN_ULONG)((((BN_ULLONG)a<<BN_BITS2)|b)/c);
258
259                         if ((i & 0xfffff) == 0) fprintf(stderr,"%d\n",i,r,rr); 
260 /*if (x == 255)
261         fprintf(stderr,"%6d/%3d = %4d %4d\n",(a<<8)|b,c,r,rr); */
262                         if (rr != r)
263                                 {
264                                 fprintf(stderr,"%8d %02X%02X / %02X = %02X %02X\n",
265                                         i,a,b,c,rr,r);
266                                 abort();
267                                 }
268                         }
269                 }
270 #endif
271         }
272
273 /* Divide h-l by d and return the result. */
274 BN_ULONG bn_div_2word(l,h,d)
275 BN_ULONG l,h,d;
276         {
277         BN_ULONG dh,dl,q,ret=0,th,tl,t,top;
278         int i,count=2;
279
280         if (d == 0) return(-1);
281
282         i=BN_num_bits_word(d);
283         if ((i != BN_BITS2) && (h > 1<<i))
284                 {
285                 fprintf(stderr,"Division would overflow\n");
286                 abort();
287                 }
288         i=BN_BITS2-i;
289         if (h >= d) h-=d;
290
291         if (i)
292                 {
293                 d<<=i;
294                 h=(h<<i)|(l>>(BN_BITS2-i));
295                 l<<=i;
296                 }
297         dh=(d&BN_MASK2h)>>BN_BITS4;
298         dl=(d&BN_MASK2l);
299         for (;;)
300                 {
301                 if ((h>>BN_BITS4) == dh)
302                         q=BN_MASK2l;
303                 else
304                         q=h/dh;
305
306                 for (;;)
307                         {
308                         t=(h-q*dh);
309                         if ((t&BN_MASK2h) ||
310                                 ((dl*q) <= (
311                                         (t<<BN_BITS4)+
312                                         ((l&BN_MASK2h)>>BN_BITS4))))
313                                 break;
314                         q--;
315                         }
316                 th=q*dh;
317                 tl=q*dl;
318                 t=(tl>>BN_BITS4);
319                 tl=(tl<<BN_BITS4)&BN_MASK2h;
320                 th+=t;
321
322                 if (l < tl) th++;
323                 l-=tl;
324                 if (h < th)
325                         {
326                         fprintf(stderr,"add back\n");
327                         h+=d;
328                         q--;
329                         }
330                 h-=th;
331
332                 if (--count == 0) break;
333
334                 ret=q<<BN_BITS4;
335                 h=((h<<BN_BITS4)|(l>>BN_BITS4))&BN_MASK2;
336                 l=(l&BN_MASK2l)<<BN_BITS4;
337                 }
338         ret|=q;
339         return(ret);
340         }