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         }