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         }