util/gdss/lib/crypto/bignum/c/KerN.c

   1 /* Copyright     Digital Equipment Corporation & INRIA     1988, 1989 */
   2 /* Last modified_on Fri Mar  2 16:49:23 GMT+1:00 1990 by herve */
   3 /*      modified_on Mon Feb 19 20:22:20 GMT+1:00 1990 by shand */
   4
   5
   6 /* KerN.c: the kernel written in C */
   7
   8 /*
   9  *      Description of types and constants.
  10  *
  11  * Several conventions are used in the commentary:
  12  *    A "BigNum" is the name for an infinite-precision number.
  13  *    Capital letters (e.g., "N") are used to refer to the value of BigNums.
  14  *    The word "digit" refers to a single BigNum digit.
  15  *    The notation "Size(N)" refers to the number of digits in N,
  16  *       which is typically passed to the subroutine as "nl".
  17  *    The notation "Length(N)" refers to the number of digits in N,
  18  *       not including any leading zeros.
  19  *    The word "Base" is used for the number 2 ** BN_DIGIT_SIZE, where
  20  *       BN_DIGIT_SIZE is the number of bits in a single BigNum digit.
  21  *    The expression "BBase(N)" is used for Base ** NumDigits(N).
  22  *    The term "leading zeros" refers to any zeros before the most
  23  *       significant digit of a number.
  24  *
  25  *
  26  * In the code, we have:
  27  *
  28  *    "nn" is a pointer to a big number,
  29  *    "nl" is the number of digits from nn,
  30  *    "d" is a digit.
  31  *
  32  */
  33
  34
  35 /*\f*/
  36
  37
  38 #define BNNMACROS_OFF
  39 #include "BigNum.h"
  40
  41
  42                 /*** copyright ***/
  43
  44 static char copyright[]="@(#)KerN.c: copyright Digital Equipment Corporation & INRIA 1988, 1989\n";
  45
  46
  47         /******* non arithmetic access to digits ********/
  48
  49
  50 void BnnSetToZero (nn, nl)
  51
  52 register BigNum         nn;
  53 register int            nl;
  54
  55 /*
  56  * Sets all the specified digits of the BigNum to 0
  57  */
  58
  59 {
  60 #ifdef NOMEM
  61     while (--nl >= 0)
  62         *(nn++) = 0;
  63 #else
  64     memset (nn, 0, nl*BN_DIGIT_SIZE/BN_BYTE_SIZE);
  65 #endif
  66 }
  67
  68                 /***************************************/
  69
  70
  71 void BnnAssign (mm, nn, nl)
  72
  73 register BigNum mm, nn;
  74 register int    nl;
  75
  76 /*
  77  * Copies N => M
  78  */
  79
  80 {
  81     if (mm < nn || mm > nn+nl)
  82 #ifdef NOMEM
  83         while (--nl >= 0)
  84             *mm++ = *nn++;
  85 #else
  86         /* be care: bcopy (SRC, DEST, L): SRC-->DEST !!! */
  87 #ifdef SOLARIS
  88         memmove (mm, nn, nl*BN_DIGIT_SIZE/BN_BYTE_SIZE);
  89 #else
  90         bcopy (nn, mm, nl*BN_DIGIT_SIZE/BN_BYTE_SIZE);
  91 #endif
  92 #endif
  93     else
  94     if (mm > nn)
  95     {
  96         nn += nl;
  97         mm += nl;
  98         while (--nl >= 0)
  99             *--mm = *--nn;
 100     }
 101 }
 102
 103                 /***************************************/
 104 /*\f*/
 105
 106
 107 void BnnSetDigit (nn, d)
 108
 109 BigNum  nn;
 110 int     d;
 111
 112 /*
 113  * Sets a single digit of N to the passed value
 114  */
 115
 116 {
 117     *nn = d;
 118 }
 119
 120                 /***************************************/
 121
 122
 123 BigNumDigit BnnGetDigit (nn)
 124
 125 BigNum  nn;
 126
 127 /*
 128  * Returns the single digit pointed by N
 129  */
 130
 131 {
 132     return (*nn);
 133 }
 134
 135                 /***************************************/
 136 /*\f*/
 137
 138
 139 BigNumLength BnnNumDigits (nn, nl)
 140
 141 register BigNum nn;
 142 register int    nl;
 143
 144 /*
 145  * Returns the number of digits of N, not counting leading zeros
 146  */
 147
 148 {
 149     nn += nl;
 150
 151     while (nl != 0 && *--nn == 0)
 152         nl--;
 153
 154     return (nl == 0 ? 1 : nl);
 155 }
 156
 157                 /***************************************/
 158
 159
 160 BigNumDigit BnnNumLeadingZeroBitsInDigit (d)
 161
 162 BigNumDigit d;
 163
 164 /*
 165  * Returns the number of leading zero bits in a digit
 166  */
 167
 168 {
 169     register BigNumDigit mask = 1 << (BN_DIGIT_SIZE - 1);
 170     register int        p = 0;
 171
 172
 173     if (d == 0)
 174         return (BN_DIGIT_SIZE);
 175
 176     while ((d & mask) == 0)
 177     {
 178         p++;
 179         mask >>= 1;
 180     }
 181
 182     return (p);
 183 }
 184
 185                 /***************************************/
 186 /*\f*/
 187
 188         /************** Predicates on one digit ***************/
 189
 190
 191 Boolean BnnDoesDigitFitInWord (d)
 192
 193 BigNumDigit d;
 194
 195 /*
 196  * Returns TRUE iff the digit can be represented in just BN_WORD_SIZE bits
 197  */
 198 {
 199     /* The C compiler must evaluate the predicate at compile time */
 200     if (BN_DIGIT_SIZE > BN_WORD_SIZE)
 201         return (d >= 1 << BN_WORD_SIZE ? FALSE : TRUE);
 202     else
 203         return (TRUE);
 204 }
 205
 206                 /***************************************/
 207
 208
 209 Boolean BnnIsDigitZero (d)
 210
 211 BigNumDigit d;
 212
 213 /* Returns TRUE iff digit = 0 */
 214
 215 {
 216     return (d == 0);
 217 }
 218
 219                 /***************************************/
 220
 221
 222 Boolean BnnIsDigitNormalized (d)
 223
 224 BigNumDigit d;
 225
 226 /*
 227  * Returns TRUE iff Base/2 <= digit < Base
 228  * i.e., if digit's leading bit is 1
 229  */
 230
 231 {
 232     return (d & (1 << (BN_DIGIT_SIZE - 1)) ? TRUE : FALSE);
 233 }
 234
 235                 /***************************************/
 236
 237
 238 Boolean BnnIsDigitOdd (d)
 239
 240 BigNumDigit d;
 241
 242 /*
 243  * Returns TRUE iff digit is odd
 244  */
 245
 246 {
 247     return (d & 1 ? TRUE : FALSE);
 248 }
 249
 250                 /***************************************/
 251
 252
 253 BigNumCmp BnnCompareDigits (d1, d2)
 254
 255 BigNumDigit d1, d2;
 256
 257 /*
 258  * Returns      BN_GREATER      if digit1 > digit2
 259  *              BN_EQUAL        if digit1 = digit2
 260  *              BN_LESS if digit1 < digit2
 261  */
 262
 263 {
 264     return (d1 > d2 ? BN_GT : (d1 == d2 ? BN_EQ : BN_LT));
 265 }
 266
 267         /***************** Logical operations ********************/
 268
 269
 270 void BnnComplement (nn, nl)
 271
 272 register BigNum nn;
 273 register int    nl;
 274
 275 /*
 276  * Performs the computation BBase(N) - N - 1 => N
 277  */
 278
 279 {
 280     while (--nl >= 0)
 281         *(nn++) ^= -1;
 282 }
 283
 284                 /***************************************/
 285 /*\f*/
 286
 287
 288 void BnnAndDigits (n, d)
 289
 290 BigNum          n;
 291 BigNumDigit     d;
 292
 293 /*
 294  * Returns the logical computation n[0] AND d in n[0]
 295  */
 296
 297 {
 298     *n &= d;
 299 }
 300
 301                 /***************************************/
 302
 303
 304 void BnnOrDigits (n, d)
 305
 306 BigNum          n;
 307 BigNumDigit     d;
 308
 309 /*
 310  * Returns the logical computation n[0] OR d2 in n[0].
 311  */
 312
 313 {
 314     *n |= d;
 315 }
 316
 317                 /***************************************/
 318
 319
 320 void BnnXorDigits (n, d)
 321
 322 BigNum          n;
 323 BigNumDigit     d;
 324
 325 /*
 326  * Returns the logical computation n[0] XOR d in n[0].
 327  */
 328
 329 {
 330     *n ^= d;
 331 }
 332
 333                 /***************************************/
 334 /*\f*/
 335
 336         /****************** Shift operations *******************/
 337
 338
 339 BigNumDigit BnnShiftLeft (mm, ml, nbits)
 340
 341 register BigNum mm;
 342 register int    ml;
 343          int    nbits;
 344
 345 /*
 346  * Shifts M left by "nbits", filling with 0s.
 347  * Returns the leftmost "nbits" of M in a digit.
 348  * Assumes 0 <= nbits < BN_DIGIT_SIZE.
 349  */
 350
 351 {
 352     register BigNumDigit res = 0, save;
 353              int         rnbits;
 354
 355
 356     if (nbits != 0)
 357     {
 358         rnbits = BN_DIGIT_SIZE - nbits;
 359
 360         while (--ml >= 0)
 361         {
 362             save = *mm;
 363             *mm++ = (save << nbits) | res;
 364             res = save >> rnbits;
 365         }
 366     }
 367
 368     return (res);
 369 }
 370
 371                 /***************************************/
 372 /*\f*/
 373
 374
 375 BigNumDigit BnnShiftRight (mm, ml, nbits)
 376
 377 register BigNum mm;
 378 register int    ml;
 379          int    nbits;
 380
 381 /*
 382  * Shifts M right by "nbits", filling with 0s.
 383  * Returns the rightmost "nbits" of M in a digit.
 384  * Assumes 0 <= nbits < BN_DIGIT_SIZE.
 385  */
 386
 387 {
 388     register BigNumDigit res = 0, save;
 389              int         lnbits;
 390
 391
 392     if (nbits != 0)
 393     {
 394         mm += ml;
 395         lnbits = BN_DIGIT_SIZE - nbits;
 396
 397         while (--ml >= 0)
 398         {
 399             save = *(--mm);
 400             *mm = (save >> nbits) | res;
 401             res = save << lnbits;
 402         }
 403     }
 404
 405     return (res);
 406 }
 407
 408                 /***************************************/
 409 /*\f*/
 410
 411
 412         /******************* Additions **************************/
 413
 414
 415 BigNumCarry BnnAddCarry (nn, nl, carryin)
 416
 417 register BigNum         nn;
 418 register int            nl;
 419          BigNumCarry    carryin;
 420
 421 /*
 422  * Performs the sum N + CarryIn => N.
 423  * Returns the CarryOut.
 424  */
 425
 426 {
 427     if (carryin == 0)
 428         return (0);
 429
 430     if (nl == 0)
 431         return (1);
 432
 433     while (--nl >= 0 && !(++(*nn++)))
 434         ;
 435
 436     return (nl >= 0 ? 0 : 1);
 437 }
 438
 439                 /***************************************/
 440 /*\f*/
 441
 442
 443 BigNumCarry BnnAdd (mm, ml, nn, nl, carryin)
 444
 445 register BigNum         mm, nn;
 446          int            ml;
 447 register int            nl;
 448          BigNumCarry    carryin;
 449
 450 /*
 451  * Performs the sum M + N + CarryIn => M.
 452  * Returns the CarryOut.
 453  * Assumes Size(M) >= Size(N).
 454  */
 455
 456 {
 457     register BigNumProduct c = carryin;
 458
 459
 460     ml -= nl;
 461     /* test computed at compile time */
 462     if (sizeof (BigNumProduct) > sizeof (BigNumDigit))
 463     {
 464         while (--nl >= 0)
 465         {
 466             c += *mm + *(nn++);
 467             *(mm++) = c;
 468             c >>= BN_DIGIT_SIZE;
 469         }
 470     }
 471     else
 472     {
 473         register BigNumProduct save;
 474
 475         while (--nl >= 0)
 476         {
 477             save = *mm;
 478             c += save;
 479             if (c < save)
 480             {
 481                 *(mm++) = *(nn++);
 482                 c = 1;
 483             }
 484             else
 485             {
 486                 save = *(nn++);
 487                 c += save;
 488                 *(mm++) = c;
 489                 c = (c < save) ? 1 : 0;
 490             }
 491         }
 492     }
 493
 494     return (BnnAddCarry (mm, ml, (BigNumCarry) c));
 495 }
 496
 497                 /***************************************/
 498 /*\f*/
 499
 500         /****************** Subtraction *************************/
 501
 502
 503
 504 BigNumCarry BnnSubtractBorrow (nn, nl, carryin)
 505
 506 register BigNum         nn;
 507 register int            nl;
 508          BigNumCarry    carryin;
 509
 510 /*
 511  * Performs the difference N + CarryIn - 1 => N.
 512  * Returns the CarryOut.
 513  */
 514
 515 {
 516     if (carryin == 1)
 517         return (1);
 518     if (nl == 0)
 519         return (0);
 520
 521 #ifdef ibm032
 522     while (--nl >= 0 && !(*nn)--)  nn++;
 523 #else
 524     while (--nl >= 0 && !((*nn++)--))
 525 #endif
 526         ;
 527
 528     return (nl >= 0 ? 1 : 0);
 529 }
 530
 531                 /***************************************/
 532 /*\f*/
 533
 534
 535 BigNumCarry BnnSubtract (mm, ml, nn, nl, carryin)
 536
 537 register BigNum         mm, nn;
 538          int            ml;
 539 register int            nl;
 540          BigNumCarry    carryin;
 541
 542 /*
 543  * Performs the difference M - N + CarryIn - 1 => M.
 544  * Returns the CarryOut.
 545  * Assumes Size(M) >= Size(N).
 546  */
 547
 548 {
 549     register BigNumProduct      c = carryin;
 550     register BigNumDigit        invn;
 551
 552
 553     ml -= nl;
 554     /* test computed at compile time */
 555     if (sizeof (BigNumProduct) > sizeof (BigNumDigit))
 556     {
 557         while (--nl >= 0)
 558         {
 559             invn = *(nn++) ^ -1;
 560             c += *mm + invn;
 561             *(mm++) = c;
 562             c >>= BN_DIGIT_SIZE;
 563         }
 564     }
 565     else
 566     {
 567         register BigNumProduct save;
 568
 569         while (--nl >= 0)
 570         {
 571             save = *mm;
 572             invn = *(nn++) ^ -1;
 573             c += save;
 574
 575             if (c < save)
 576             {
 577                 *(mm++) = invn;
 578                 c = 1;
 579             }
 580             else
 581             {
 582                 c += invn;
 583                 *(mm++) = c;
 584                 c = (c < invn) ? 1 : 0;
 585             }
 586         }
 587     }
 588
 589     return (BnnSubtractBorrow (mm, ml, (BigNumCarry) c)); }
 590
 591
 592                 /***************************************/
 593 /*\f */
 594
 595         /***************** Multiplication ************************/
 596
 597
 598 BigNumCarry BnnMultiplyDigit (pp, pl, mm, ml, d)
 599
 600 register BigNum         pp, mm;
 601          int            pl, ml;
 602          BigNumDigit    d;
 603
 604 /*
 605  * Performs the product:
 606  * Q = P + M * d
 607  * BB = BBase(P)
 608  * Q mod BB => P
 609  * Q div BB => CarryOut
 610  * Returns the CarryOut.
 611  * Assumes Size(P) >= Size(M) + 1.
 612  */
 613
 614 {
 615     register BigNumProduct c = 0;
 616
 617
 618     if (d == 0)
 619         return (0);
 620
 621     if (d == 1)
 622         return (BnnAdd (pp, pl, mm, ml, (BigNumCarry) 0));
 623
 624     pl -= ml;
 625     /* test computed at compile time */
 626     if (sizeof (BigNumProduct) > sizeof (BigNumDigit))
 627     {
 628         while (ml != 0)
 629         {
 630             ml--;
 631             c += *pp + (d * (*(mm++)));
 632             *(pp++) = c;
 633             c >>= BN_DIGIT_SIZE;
 634         }
 635
 636         while (pl != 0) {
 637             pl--;
 638             c += *pp;
 639             *(pp++) = c;
 640             c >>= BN_DIGIT_SIZE;
 641         }
 642
 643         return (c);
 644     }
 645     else
 646     {
 647 /* help for stupid compilers--may actually be counter
 648    productive on pipelined machines with decent register allocation!! */
 649 #define m_digit X0
 650 #define X3 Lm
 651 #define X1 Hm
 652         register BigNumDigit Lm, Hm, Ld, Hd, X0, X2 /*, X1, X3 */;
 653
 654         Ld = d & ((1 << (BN_DIGIT_SIZE / 2)) -1);
 655         Hd = d >> (BN_DIGIT_SIZE / 2);
 656         while (ml != 0)
 657         {
 658             ml--;
 659             m_digit = *mm++;
 660             Lm = m_digit & ((1 << (BN_DIGIT_SIZE / 2)) -1);
 661             Hm = m_digit >> (BN_DIGIT_SIZE / 2);
 662             X0 = Ld * Lm;
 663             X2 = Hd * Lm;
 664             X3 = Hd * Hm;
 665             X1 = Ld * Hm;
 666
 667             if ((c += X0) < X0) X3++;
 668             if ((X1 += X2) < X2) X3 += (1<<(BN_DIGIT_SIZE / 2));
 669             X3 += (X1 >> (BN_DIGIT_SIZE / 2));
 670             X1 <<= (BN_DIGIT_SIZE / 2);
 671             if ((c += X1) < X1) X3++;
 672             if ((*pp += c) < c) X3++;
 673             pp++;
 674
 675             c = X3;
 676 #undef m_digit
 677 #undef X1
 678 #undef X3
 679         }
 680
 681         X0 = *pp;
 682         c += X0;
 683         *(pp++) = c;
 684
 685         if (c >= X0)
 686             return (0);
 687
 688         pl--;
 689         while (pl != 0 && !(++(*pp++)))
 690             pl--;
 691
 692         return (pl != 0 ? 0 : 1);
 693     }
 694 }
 695
 696 #ifdef mips
 697 BigNumCarry BnnMultiply2Digit (pp, pl, mm, ml, d0, d1)
 698
 699 register BigNum         pp, mm;
 700 register int            pl, ml;
 701          BigNumDigit    d0, d1;
 702
 703 /*
 704  * Provided for compatibility with mips assembler implementation.
 705  * Performs the product:
 706  * Q = P + M * d0_d1
 707  * BB = BBase(P)
 708  * Q mod BB => P
 709  * Q div BB => CarryOut
 710  * Returns the CarryOut.
 711  * Assumes Size(P) >= Size(M) + 1.
 712  */
 713
 714 {
 715     return
 716         BnnMultiplyDigit (pp, pl, mm, ml, d0)
 717         + BnnMultiplyDigit (pp+1, pl-1, mm, ml, d1);
 718 }
 719 #endif /* mips */
 720
 721                 /***************************************/
 722 /*\f*/
 723
 724         /********************** Division *************************/
 725
 726
 727                 /* xh:xl -= yh:yl */
 728 #define SUB(xh,xl,yh,yl)        if (yl > xl) {xl -= yl; xh -= yh + 1;}\
 729                                 else         {xl -= yl; xh -= yh;}
 730
 731 #define LOW(x)                  (x & ((1 << (BN_DIGIT_SIZE / 2)) -1))
 732 #define HIGH(x)                 (x >> (BN_DIGIT_SIZE / 2))
 733 #define L2H(x)                  (x << (BN_DIGIT_SIZE / 2))
 734
 735
 736 BigNumDigit BnnDivideDigit (qq, nn, nl, d)
 737
 738 register BigNum         qq, nn;
 739 register int            nl;
 740          BigNumDigit    d;
 741
 742 /* Performs the quotient: N div d => Q
 743  * Returns R = N mod d
 744  * Assumes leading digit of N < d, and d > 0.
 745  */
 746
 747 {
 748     /* test computed at compile time */
 749     if (sizeof (BigNumProduct) > sizeof (BigNumDigit))
 750     {
 751         register BigNumProduct quad;
 752
 753
 754         nn += nl;
 755         nl--;
 756         qq += nl;
 757         quad = *(--nn);
 758
 759         while (nl != 0)
 760         {
 761             nl--;
 762             quad = (quad << BN_DIGIT_SIZE) | *(--nn);
 763             *(--qq) = quad / d;
 764             quad = quad % d;
 765         }
 766
 767         return (quad);
 768     }
 769     else
 770     {
 771         int             k;
 772         int             orig_nl;
 773         BigNumDigit     rh;             /* Two halves of current remainder */
 774         BigNumDigit     rl;             /* Correspond to quad above */
 775         register BigNumDigit qa;        /* Current appr. to quotient */
 776         register BigNumDigit ph, pl;    /* product of c and qa */
 777         BigNumDigit     ch, cl, prev_qq;
 778
 779
 780         /* Normalize divisor */
 781         k = BnnNumLeadingZeroBitsInDigit (d);
 782         if (k != 0)
 783         {
 784             prev_qq = qq[-1];
 785             orig_nl = nl;
 786             d <<= k;
 787             BnnShiftLeft (nn, nl, k);
 788         }
 789
 790         nn += nl;
 791         nl--;
 792         qq += nl;
 793
 794         ch = HIGH (d);
 795         cl = LOW (d);
 796
 797         rl = *(--nn);
 798
 799         while (nl != 0)
 800         {
 801             nl--;
 802             rh = rl;
 803             rl = *(--nn);
 804             qa = rh / ch;       /* appr. quotient */
 805
 806             /* Compute ph, pl */
 807             pl = cl * qa;
 808             ph = ch * qa;
 809             ph += HIGH (pl);
 810             pl = L2H (pl);
 811
 812             /* While ph:pl > rh:rl, decrement qa, adjust qh:ql */
 813             while (ph > rh || ph == rh && pl > rl)
 814             {
 815                 qa--;
 816                 SUB (ph, pl, ch, L2H (cl));
 817             }
 818
 819             SUB (rh, rl, ph, pl);
 820
 821             /* Top half of quotient is correct; save it */
 822             *(--qq) = L2H (qa);
 823             qa = (L2H (rh) | HIGH (rl)) / ch;
 824
 825             /* Approx low half of q */
 826             /* Compute ph, pl, again */
 827             pl = cl * qa;
 828             ph = ch * qa;
 829             ph += HIGH (pl);
 830             pl = LOW (pl) | L2H (LOW (ph));
 831             ph = HIGH (ph);
 832
 833             /* While ph:pl > rh:rl, decrement qa, adjust qh:ql */
 834             while (ph > rh || ph == rh && pl > rl)
 835             {
 836                 qa--;
 837                 SUB (ph, pl, 0, d);
 838             }
 839
 840             /* Subtract ph:pl from rh:rl; we know rh will be 0 */
 841             rl -= pl;
 842             *qq |= qa;
 843         }
 844
 845         /* Denormalize dividend */
 846         if (k != 0) {
 847                 if((qq > nn) && (qq < &nn[orig_nl])) {
 848                         /* Overlap between qq and nn. Care of *qq! */
 849                         orig_nl = (qq - nn);
 850                         BnnShiftRight (nn, orig_nl, k);
 851                         nn[orig_nl - 1] = prev_qq;
 852                 } else if(qq == nn) {
 853                         BnnShiftRight(&nn[orig_nl - 1], 1, k);
 854                 } else {
 855                         BnnShiftRight (nn, orig_nl, k);
 856         }       }
 857         return (rl >> k);
 858     }
 859 }
 860
 861                 /***************************************/
 862
 863