util/rsaref/nn.c

   1 /* NN.C - natural numbers routines
   2  */
   3
   4 /* Copyright (C) RSA Laboratories, a division of RSA Data Security,
   5      Inc., created 1991. All rights reserved.
   6  */
   7
   8 #include "global.h"
   9 #include "rsaref.h"
  10 #include "nn.h"
  11 #include "digit.h"
  12
  13 static NN_DIGIT NN_AddDigitMult PROTO_LIST
  14   ((NN_DIGIT *, NN_DIGIT *, NN_DIGIT, NN_DIGIT *, unsigned int));
  15 static NN_DIGIT NN_SubDigitMult PROTO_LIST
  16   ((NN_DIGIT *, NN_DIGIT *, NN_DIGIT, NN_DIGIT *, unsigned int));
  17
  18 static unsigned int NN_DigitBits PROTO_LIST ((NN_DIGIT));
  19
  20 /* Decodes character string b into a, where character string is ordered
  21    from most to least significant.
  22
  23    Lengths: a[digits], b[len].
  24    Assumes b[i] = 0 for i < len - digits * NN_DIGIT_LEN. (Otherwise most
  25    significant bytes are truncated.)
  26  */
  27 void NN_Decode (a, digits, b, len)
  28 NN_DIGIT *a;
  29 unsigned char *b;
  30 unsigned int digits, len;
  31 {
  32   NN_DIGIT t;
  33   int j;
  34   unsigned int i, u;
  35
  36   for (i = 0, j = len - 1; i < digits && j >= 0; i++) {
  37     t = 0;
  38     for (u = 0; j >= 0 && u < NN_DIGIT_BITS; j--, u += 8)
  39       t |= ((NN_DIGIT)b[j]) << u;
  40     a[i] = t;
  41   }
  42
  43   for (; i < digits; i++)
  44     a[i] = 0;
  45 }
  46
  47 /* Encodes b into character string a, where character string is ordered
  48    from most to least significant.
  49
  50    Lengths: a[len], b[digits].
  51    Assumes NN_Bits (b, digits) <= 8 * len. (Otherwise most significant
  52    digits are truncated.)
  53  */
  54 void NN_Encode (a, len, b, digits)
  55 NN_DIGIT *b;
  56 unsigned char *a;
  57 unsigned int digits, len;
  58 {
  59   NN_DIGIT t;
  60   int j;
  61   unsigned int i, u;
  62
  63   for (i = 0, j = len - 1; i < digits && j >= 0; i++) {
  64     t = b[i];
  65     for (u = 0; j >= 0 && u < NN_DIGIT_BITS; j--, u += 8)
  66       a[j] = (unsigned char)(t >> u);
  67   }
  68
  69   for (; j >= 0; j--)
  70     a[j] = 0;
  71 }
  72
  73 /* Assigns a = b.
  74
  75    Lengths: a[digits], b[digits].
  76  */
  77 void NN_Assign (a, b, digits)
  78 NN_DIGIT *a, *b;
  79 unsigned int digits;
  80 {
  81   unsigned int i;
  82
  83   for (i = 0; i < digits; i++)
  84     a[i] = b[i];
  85 }
  86
  87 /* Assigns a = 0.
  88
  89    Lengths: a[digits].
  90  */
  91 void NN_AssignZero (a, digits)
  92 NN_DIGIT *a;
  93 unsigned int digits;
  94 {
  95   unsigned int i;
  96
  97   for (i = 0; i < digits; i++)
  98     a[i] = 0;
  99 }
 100
 101 /* Assigns a = 2^b.
 102
 103    Lengths: a[digits].
 104    Requires b < digits * NN_DIGIT_BITS.
 105  */
 106 void NN_Assign2Exp (a, b, digits)
 107 NN_DIGIT *a;
 108 unsigned int b, digits;
 109 {
 110   NN_AssignZero (a, digits);
 111
 112   if (b >= digits * NN_DIGIT_BITS)
 113     return;
 114
 115   a[b / NN_DIGIT_BITS] = (NN_DIGIT)1 << (b % NN_DIGIT_BITS);
 116 }
 117
 118 /* Computes a = b + c. Returns carry.
 119
 120    Lengths: a[digits], b[digits], c[digits].
 121  */
 122 NN_DIGIT NN_Add (a, b, c, digits)
 123 NN_DIGIT *a, *b, *c;
 124 unsigned int digits;
 125 {
 126   NN_DIGIT ai, carry;
 127   unsigned int i;
 128
 129   carry = 0;
 130
 131   for (i = 0; i < digits; i++) {
 132     if ((ai = b[i] + carry) < carry)
 133       ai = c[i];
 134     else if ((ai += c[i]) < c[i])
 135       carry = 1;
 136     else
 137       carry = 0;
 138     a[i] = ai;
 139   }
 140
 141   return (carry);
 142 }
 143
 144 /* Computes a = b - c. Returns borrow.
 145
 146    Lengths: a[digits], b[digits], c[digits].
 147  */
 148 NN_DIGIT NN_Sub (a, b, c, digits)
 149 NN_DIGIT *a, *b, *c;
 150 unsigned int digits;
 151 {
 152   NN_DIGIT ai, borrow;
 153   unsigned int i;
 154
 155   borrow = 0;
 156
 157   for (i = 0; i < digits; i++) {
 158     if ((ai = b[i] - borrow) > (MAX_NN_DIGIT - borrow))
 159       ai = MAX_NN_DIGIT - c[i];
 160     else if ((ai -= c[i]) > (MAX_NN_DIGIT - c[i]))
 161       borrow = 1;
 162     else
 163       borrow = 0;
 164     a[i] = ai;
 165   }
 166
 167   return (borrow);
 168 }
 169
 170 /* Computes a = b * c.
 171
 172    Lengths: a[2*digits], b[digits], c[digits].
 173    Assumes digits < MAX_NN_DIGITS.
 174  */
 175 void NN_Mult (a, b, c, digits)
 176 NN_DIGIT *a, *b, *c;
 177 unsigned int digits;
 178 {
 179   NN_DIGIT t[2*MAX_NN_DIGITS];
 180   unsigned int bDigits, cDigits, i;
 181
 182   NN_AssignZero (t, 2 * digits);
 183
 184   bDigits = NN_Digits (b, digits);
 185   cDigits = NN_Digits (c, digits);
 186
 187   for (i = 0; i < bDigits; i++)
 188     t[i+cDigits] += NN_AddDigitMult (&t[i], &t[i], b[i], c, cDigits);
 189
 190   NN_Assign (a, t, 2 * digits);
 191
 192   /* Zeroize potentially sensitive information.
 193    */
 194   R_memset ((POINTER)t, 0, sizeof (t));
 195 }
 196
 197 /* Computes a = b * 2^c (i.e., shifts left c bits), returning carry.
 198
 199    Lengths: a[digits], b[digits].
 200    Requires c < NN_DIGIT_BITS.
 201  */
 202 NN_DIGIT NN_LShift (a, b, c, digits)
 203 NN_DIGIT *a, *b;
 204 unsigned int c, digits;
 205 {
 206   NN_DIGIT bi, carry;
 207   unsigned int i, t;
 208
 209   if (c >= NN_DIGIT_BITS)
 210     return (0);
 211
 212   t = NN_DIGIT_BITS - c;
 213
 214   carry = 0;
 215
 216   for (i = 0; i < digits; i++) {
 217     bi = b[i];
 218     a[i] = (bi << c) | carry;
 219     carry = c ? (bi >> t) : 0;
 220   }
 221
 222   return (carry);
 223 }
 224
 225 /* Computes a = c div 2^c (i.e., shifts right c bits), returning carry.
 226
 227    Lengths: a[digits], b[digits].
 228    Requires: c < NN_DIGIT_BITS.
 229  */
 230 NN_DIGIT NN_RShift (a, b, c, digits)
 231 NN_DIGIT *a, *b;
 232 unsigned int c, digits;
 233 {
 234   NN_DIGIT bi, carry;
 235   int i;
 236   unsigned int t;
 237
 238   if (c >= NN_DIGIT_BITS)
 239     return (0);
 240
 241   t = NN_DIGIT_BITS - c;
 242
 243   carry = 0;
 244
 245   for (i = digits - 1; i >= 0; i--) {
 246     bi = b[i];
 247     a[i] = (bi >> c) | carry;
 248     carry = c ? (bi << t) : 0;
 249   }
 250
 251   return (carry);
 252 }
 253
 254 /* Computes a = c div d and b = c mod d.
 255
 256    Lengths: a[cDigits], b[dDigits], c[cDigits], d[dDigits].
 257    Assumes d > 0, cDigits < 2 * MAX_NN_DIGITS,
 258            dDigits < MAX_NN_DIGITS.
 259  */
 260 void NN_Div (a, b, c, cDigits, d, dDigits)
 261 NN_DIGIT *a, *b, *c, *d;
 262 unsigned int cDigits, dDigits;
 263 {
 264   NN_DIGIT ai, cc[2*MAX_NN_DIGITS+1], dd[MAX_NN_DIGITS], t;
 265   int i;
 266   unsigned int ddDigits, shift;
 267
 268   ddDigits = NN_Digits (d, dDigits);
 269   if (ddDigits == 0)
 270     return;
 271
 272   /* Normalize operands.
 273    */
 274   shift = NN_DIGIT_BITS - NN_DigitBits (d[ddDigits-1]);
 275   NN_AssignZero (cc, ddDigits);
 276   cc[cDigits] = NN_LShift (cc, c, shift, cDigits);
 277   NN_LShift (dd, d, shift, ddDigits);
 278   t = dd[ddDigits-1];
 279
 280   NN_AssignZero (a, cDigits);
 281
 282   for (i = cDigits-ddDigits; i >= 0; i--) {
 283     /* Underestimate quotient digit and subtract.
 284      */
 285     if (t == MAX_NN_DIGIT)
 286       ai = cc[i+ddDigits];
 287     else
 288       NN_DigitDiv (&ai, &cc[i+ddDigits-1], t + 1);
 289     cc[i+ddDigits] -= NN_SubDigitMult (&cc[i], &cc[i], ai, dd, ddDigits);
 290
 291     /* Correct estimate.
 292      */
 293     while (cc[i+ddDigits] || (NN_Cmp (&cc[i], dd, ddDigits) >= 0)) {
 294       ai++;
 295       cc[i+ddDigits] -= NN_Sub (&cc[i], &cc[i], dd, ddDigits);
 296     }
 297
 298     a[i] = ai;
 299   }
 300
 301   /* Restore result.
 302    */
 303   NN_AssignZero (b, dDigits);
 304   NN_RShift (b, cc, shift, ddDigits);
 305
 306   /* Zeroize potentially sensitive information.
 307    */
 308   R_memset ((POINTER)cc, 0, sizeof (cc));
 309   R_memset ((POINTER)dd, 0, sizeof (dd));
 310 }
 311
 312 /* Computes a = b mod c.
 313
 314    Lengths: a[cDigits], b[bDigits], c[cDigits].
 315    Assumes c > 0, bDigits < 2 * MAX_NN_DIGITS, cDigits < MAX_NN_DIGITS.
 316  */
 317 void NN_Mod (a, b, bDigits, c, cDigits)
 318 NN_DIGIT *a, *b, *c;
 319 unsigned int bDigits, cDigits;
 320 {
 321   NN_DIGIT t[2 * MAX_NN_DIGITS];
 322
 323   NN_Div (t, a, b, bDigits, c, cDigits);
 324
 325   /* Zeroize potentially sensitive information.
 326    */
 327   R_memset ((POINTER)t, 0, sizeof (t));
 328 }
 329
 330 /* Computes a = b * c mod d.
 331
 332    Lengths: a[digits], b[digits], c[digits], d[digits].
 333    Assumes d > 0, digits < MAX_NN_DIGITS.
 334  */
 335 void NN_ModMult (a, b, c, d, digits)
 336 NN_DIGIT *a, *b, *c, *d;
 337 unsigned int digits;
 338 {
 339   NN_DIGIT t[2*MAX_NN_DIGITS];
 340
 341   NN_Mult (t, b, c, digits);
 342   NN_Mod (a, t, 2 * digits, d, digits);
 343
 344   /* Zeroize potentially sensitive information.
 345    */
 346   R_memset ((POINTER)t, 0, sizeof (t));
 347 }
 348
 349 /* Computes a = b^c mod d.
 350
 351    Lengths: a[dDigits], b[dDigits], c[cDigits], d[dDigits].
 352    Assumes d > 0, cDigits > 0, dDigits < MAX_NN_DIGITS.
 353  */
 354 void NN_ModExp (a, b, c, cDigits, d, dDigits)
 355 NN_DIGIT *a, *b, *c, *d;
 356 unsigned int cDigits, dDigits;
 357 {
 358   NN_DIGIT bPower[3][MAX_NN_DIGITS], ci, t[MAX_NN_DIGITS];
 359   int i;
 360   unsigned int ciBits, j, s;
 361
 362   /* Store b, b^2 mod d, and b^3 mod d.
 363    */
 364   NN_Assign (bPower[0], b, dDigits);
 365   NN_ModMult (bPower[1], bPower[0], b, d, dDigits);
 366   NN_ModMult (bPower[2], bPower[1], b, d, dDigits);
 367
 368   NN_ASSIGN_DIGIT (t, 1, dDigits);
 369
 370   cDigits = NN_Digits (c, cDigits);
 371   for (i = cDigits - 1; i >= 0; i--) {
 372     ci = c[i];
 373     ciBits = NN_DIGIT_BITS;
 374
 375     /* Scan past leading zero bits of most significant digit.
 376      */
 377     if (i == (int)(cDigits - 1)) {
 378       while (! DIGIT_2MSB (ci)) {
 379         ci <<= 2;
 380         ciBits -= 2;
 381       }
 382     }
 383
 384     for (j = 0; j < ciBits; j += 2, ci <<= 2) {
 385       /* Compute t = t^4 * b^s mod d, where s = two MSB's of ci.
 386        */
 387       NN_ModMult (t, t, t, d, dDigits);
 388       NN_ModMult (t, t, t, d, dDigits);
 389       if ((s = DIGIT_2MSB (ci)) != 0)
 390         NN_ModMult (t, t, bPower[s-1], d, dDigits);
 391     }
 392   }
 393
 394   NN_Assign (a, t, dDigits);
 395
 396   /* Zeroize potentially sensitive information.
 397    */
 398   R_memset ((POINTER)bPower, 0, sizeof (bPower));
 399   R_memset ((POINTER)t, 0, sizeof (t));
 400 }
 401
 402 /* Compute a = 1/b mod c, assuming inverse exists.
 403
 404    Lengths: a[digits], b[digits], c[digits].
 405    Assumes gcd (b, c) = 1, digits < MAX_NN_DIGITS.
 406  */
 407 void NN_ModInv (a, b, c, digits)
 408 NN_DIGIT *a, *b, *c;
 409 unsigned int digits;
 410 {
 411   NN_DIGIT q[MAX_NN_DIGITS], t1[MAX_NN_DIGITS], t3[MAX_NN_DIGITS],
 412     u1[MAX_NN_DIGITS], u3[MAX_NN_DIGITS], v1[MAX_NN_DIGITS],
 413     v3[MAX_NN_DIGITS], w[2*MAX_NN_DIGITS];
 414   int u1Sign;
 415
 416   /* Apply extended Euclidean algorithm, modified to avoid negative
 417      numbers.
 418    */
 419   NN_ASSIGN_DIGIT (u1, 1, digits);
 420   NN_AssignZero (v1, digits);
 421   NN_Assign (u3, b, digits);
 422   NN_Assign (v3, c, digits);
 423   u1Sign = 1;
 424
 425   while (! NN_Zero (v3, digits)) {
 426     NN_Div (q, t3, u3, digits, v3, digits);
 427     NN_Mult (w, q, v1, digits);
 428     NN_Add (t1, u1, w, digits);
 429     NN_Assign (u1, v1, digits);
 430     NN_Assign (v1, t1, digits);
 431     NN_Assign (u3, v3, digits);
 432     NN_Assign (v3, t3, digits);
 433     u1Sign = -u1Sign;
 434   }
 435
 436   /* Negate result if sign is negative.
 437     */
 438   if (u1Sign < 0)
 439     NN_Sub (a, c, u1, digits);
 440   else
 441     NN_Assign (a, u1, digits);
 442
 443   /* Zeroize potentially sensitive information.
 444    */
 445   R_memset ((POINTER)q, 0, sizeof (q));
 446   R_memset ((POINTER)t1, 0, sizeof (t1));
 447   R_memset ((POINTER)t3, 0, sizeof (t3));
 448   R_memset ((POINTER)u1, 0, sizeof (u1));
 449   R_memset ((POINTER)u3, 0, sizeof (u3));
 450   R_memset ((POINTER)v1, 0, sizeof (v1));
 451   R_memset ((POINTER)v3, 0, sizeof (v3));
 452   R_memset ((POINTER)w, 0, sizeof (w));
 453 }
 454
 455 /* Computes a = gcd(b, c).
 456
 457    Lengths: a[digits], b[digits], c[digits].
 458    Assumes b > c, digits < MAX_NN_DIGITS.
 459  */
 460 void NN_Gcd (a, b, c, digits)
 461 NN_DIGIT *a, *b, *c;
 462 unsigned int digits;
 463 {
 464   NN_DIGIT t[MAX_NN_DIGITS], u[MAX_NN_DIGITS], v[MAX_NN_DIGITS];
 465
 466   NN_Assign (u, b, digits);
 467   NN_Assign (v, c, digits);
 468
 469   while (! NN_Zero (v, digits)) {
 470     NN_Mod (t, u, digits, v, digits);
 471     NN_Assign (u, v, digits);
 472     NN_Assign (v, t, digits);
 473   }
 474
 475   NN_Assign (a, u, digits);
 476
 477   /* Zeroize potentially sensitive information.
 478    */
 479   R_memset ((POINTER)t, 0, sizeof (t));
 480   R_memset ((POINTER)u, 0, sizeof (u));
 481   R_memset ((POINTER)v, 0, sizeof (v));
 482 }
 483
 484 /* Returns sign of a - b.
 485
 486    Lengths: a[digits], b[digits].
 487  */
 488 int NN_Cmp (a, b, digits)
 489 NN_DIGIT *a, *b;
 490 unsigned int digits;
 491 {
 492   int i;
 493
 494   for (i = digits - 1; i >= 0; i--) {
 495     if (a[i] > b[i])
 496       return (1);
 497     if (a[i] < b[i])
 498       return (-1);
 499   }
 500
 501   return (0);
 502 }
 503
 504 /* Returns nonzero iff a is zero.
 505
 506    Lengths: a[digits].
 507  */
 508 int NN_Zero (a, digits)
 509 NN_DIGIT *a;
 510 unsigned int digits;
 511 {
 512   unsigned int i;
 513
 514   for (i = 0; i < digits; i++)
 515     if (a[i])
 516       return (0);
 517
 518   return (1);
 519 }
 520
 521 /* Returns the significant length of a in bits.
 522
 523    Lengths: a[digits].
 524  */
 525 unsigned int NN_Bits (a, digits)
 526 NN_DIGIT *a;
 527 unsigned int digits;
 528 {
 529   if ((digits = NN_Digits (a, digits)) == 0)
 530     return (0);
 531
 532   return ((digits - 1) * NN_DIGIT_BITS + NN_DigitBits (a[digits-1]));
 533 }
 534
 535 /* Returns the significant length of a in digits.
 536
 537    Lengths: a[digits].
 538  */
 539 unsigned int NN_Digits (a, digits)
 540 NN_DIGIT *a;
 541 unsigned int digits;
 542 {
 543   int i;
 544
 545   for (i = digits - 1; i >= 0; i--)
 546     if (a[i])
 547       break;
 548
 549   return (i + 1);
 550 }
 551
 552 /* Computes a = b + c*d, where c is a digit. Returns carry.
 553
 554    Lengths: a[digits], b[digits], d[digits].
 555  */
 556 static NN_DIGIT NN_AddDigitMult (a, b, c, d, digits)
 557 NN_DIGIT *a, *b, c, *d;
 558 unsigned int digits;
 559 {
 560   NN_DIGIT carry, t[2];
 561   unsigned int i;
 562
 563   if (c == 0)
 564     return (0);
 565
 566   carry = 0;
 567   for (i = 0; i < digits; i++) {
 568     NN_DigitMult (t, c, d[i]);
 569     if ((a[i] = b[i] + carry) < carry)
 570       carry = 1;
 571     else
 572       carry = 0;
 573     if ((a[i] += t[0]) < t[0])
 574       carry++;
 575     carry += t[1];
 576   }
 577
 578   return (carry);
 579 }
 580
 581 /* Computes a = b - c*d, where c is a digit. Returns borrow.
 582
 583    Lengths: a[digits], b[digits], d[digits].
 584  */
 585 static NN_DIGIT NN_SubDigitMult (a, b, c, d, digits)
 586 NN_DIGIT *a, *b, c, *d;
 587 unsigned int digits;
 588 {
 589   NN_DIGIT borrow, t[2];
 590   unsigned int i;
 591
 592   if (c == 0)
 593     return (0);
 594
 595   borrow = 0;
 596   for (i = 0; i < digits; i++) {
 597     NN_DigitMult (t, c, d[i]);
 598     if ((a[i] = b[i] - borrow) > (MAX_NN_DIGIT - borrow))
 599       borrow = 1;
 600     else
 601       borrow = 0;
 602     if ((a[i] -= t[0]) > (MAX_NN_DIGIT - t[0]))
 603       borrow++;
 604     borrow += t[1];
 605   }
 606
 607   return (borrow);
 608 }
 609
 610 /* Returns the significant length of a in bits, where a is a digit.
 611  */
 612 static unsigned int NN_DigitBits (a)
 613 NN_DIGIT a;
 614 {
 615   unsigned int i;
 616
 617   for (i = 0; i < NN_DIGIT_BITS; i++, a >>= 1)
 618     if (a == 0)
 619       break;
 620
 621   return (i);
 622 }