![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to mpilib.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [PGP] / pgp / src|
Return to mpilib.c CVS log | Up to [PGP] / pgp / src |
1.1 root 1: /* C source code for multiprecision arithmetic library routines.
2: Implemented Nov 86 by Philip Zimmermann
3: Last revised 27 Nov 91 by PRZ
4:
5: Boulder Software Engineering
6: 3021 Eleventh Street
7: Boulder, CO 80304
8: (303) 541-0140
9:
10: (c) Copyright 1986-92 by Philip Zimmermann. All rights reserved.
11: The author assumes no liability for damages resulting from the use
12: of this software, even if the damage results from defects in this
13: software. No warranty is expressed or implied. The use of this
14: cryptographic software for developing weapon systems is expressly
15: forbidden.
16:
17: These routines implement all of the multiprecision arithmetic
18: necessary for number-theoretic cryptographic algorithms such as
19: ElGamal, Diffie-Hellman, Rabin, or factoring studies for large
20: composite numbers, as well as Rivest-Shamir-Adleman (RSA) public
21: key cryptography.
22:
23: Although originally developed in Microsoft C for the IBM PC, this code
24: contains few machine dependancies. It assumes 2's complement
25: arithmetic. It can be adapted to 8-bit, 16-bit, or 32-bit machines,
26: lowbyte-highbyte order or highbyte-lowbyte order. This version
27: has been converted to ANSI C.
28:
29:
30: The internal representation for these extended precision integer
31: "registers" is an array of "units". A unit is a machine word, which
32: is either an 8-bit byte, a 16-bit unsigned integer, or a 32-bit
33: unsigned integer, depending on the machine's word size. For example,
34: an IBM PC or AT uses a unit size of 16 bits. To perform arithmetic
35: on these huge precision integers, we pass pointers to these unit
36: arrays to various subroutines. A pointer to an array of units is of
37: type unitptr. This is a pointer to a huge integer "register".
38:
39: When calling a subroutine, we always pass a pointer to the BEGINNING
40: of the array of units, regardless of the byte order of the machine.
41: On a lowbyte-first machine, such as the Intel 80x86, this unitptr
42: points to the LEAST significant unit, and the array of units increases
43: significance to the right. On a highbyte-first machine, such as the
44: Motorola 680x0, this unitptr points to the MOST significant unit, and
45: the array of units decreases significance to the right.
46:
47: Modified 8 Apr 92 - HAJK
48: Implement new VAX/VMS primitive support.
49:
50: */
51:
52: /* #define COUNTMULTS */ /* count modmults for performance studies */
53:
54: #ifdef DEBUG
55: #ifdef MSDOS
56: #include <conio.h>
57: #define poll_for_break() {while (kbhit()) getch();}
58: #endif /* MSDOS */
59: #endif /* DEBUG */
60:
61: #ifndef poll_for_break
62: #define poll_for_break() /* stub */
63: #endif
64:
65: #include "mpilib.h"
66:
67: short global_precision = 0; /* units of precision for all routines */
68: /* global_precision is the unit precision last set by set_precision.
69: Initially, set_precision() should be called to define global_precision
70: before using any of these other multiprecision library routines.
71: i.e.: set_precision(MAX_UNIT_PRECISION);
72: */
73:
74: #ifdef PORTABLE
75: /*************** multiprecision library primitives ****************/
76: /* The following portable C primitives should be recoded in assembly.
77: The equivalent assembly primitives are externally defined under
78: different names, unless PORTABLE is defined. See the header file
79: "mpilib.h" for further details.
80:
81: The carry bit mechanism we use for these portable primitives
82: will only work if the size of unit is smaller than the size of
83: a long integer. In most cases, this means UNITSIZE must be
84: less than 32 bits -- if you use the C portable primitives.
85: */
86:
87: typedef unsigned long int ulint;
88: #define carrybit ((ulint) 1 << UNITSIZE)
89: /* ...assumes sizeof(unit) < sizeof(unsigned long) */
90:
91: boolean mp_addc
92: (register unitptr r1,register unitptr r2,register boolean carry)
93: /* multiprecision add with carry r2 to r1, result in r1 */
94: /* carry is incoming carry flag-- value should be 0 or 1 */
95: { register ulint x; /* won't work if sizeof(unit)==sizeof(long) */
96: short precision; /* number of units to add */
97: precision = global_precision;
98: make_lsbptr(r1,precision);
99: make_lsbptr(r2,precision);
100: while (precision--)
101: { x = (ulint) *r1 + (ulint) *post_higherunit(r2) + (ulint) carry;
102: *post_higherunit(r1) = x;
103: carry = ((x & carrybit) != 0L);
104: }
105: return(carry); /* return the final carry flag bit */
106: } /* mp_addc */
107:
108:
109: boolean mp_subb
110: (register unitptr r1,register unitptr r2,register boolean borrow)
111: /* multiprecision subtract with borrow, r2 from r1, result in r1 */
112: /* borrow is incoming borrow flag-- value should be 0 or 1 */
113: { register ulint x; /* won't work if sizeof(unit)==sizeof(long) */
114: short precision; /* number of units to subtract */
115: precision = global_precision;
116: make_lsbptr(r1,precision);
117: make_lsbptr(r2,precision);
118: while (precision--)
119: { x = (ulint) *r1 - (ulint) *post_higherunit(r2) - (ulint) borrow;
120: *post_higherunit(r1) = x;
121: borrow = ((x & carrybit) != 0L);
122: }
123: return(borrow); /* return the final carry/borrow flag bit */
124: } /* mp_subb */
125:
126: #undef carrybit
127:
128:
129: boolean mp_rotate_left(register unitptr r1,register boolean carry)
130: /* multiprecision rotate left 1 bit with carry, result in r1. */
131: /* carry is incoming carry flag-- value should be 0 or 1 */
132: { register short precision; /* number of units to rotate */
133: register boolean nextcarry;
134: precision = global_precision;
135: make_lsbptr(r1,precision);
136: while (precision--)
137: {
138: nextcarry = (((signedunit) *r1) < 0);
139: /* nextcarry = ((*r1 & uppermostbit) != 0); */
140: *r1 <<= 1 ;
141: if (carry) *r1 |= 1;
142: carry = nextcarry;
143: pre_higherunit(r1);
144: }
145: return(nextcarry); /* return the final carry flag bit */
146: } /* mp_rotate_left */
147:
148: /************** end of primitives that should be in assembly *************/
149: #endif /* PORTABLE */
150:
151:
152: /* The mp_shift_right_bits function is not called in any time-critical
153: situations in public-key cryptographic functions, so it doesn't
154: need to be coded in assembly language.
155: */
156: void mp_shift_right_bits(register unitptr r1,register short bits)
157: /* multiprecision shift right bits, result in r1.
158: bits is how many bits to shift, must be < UNITSIZE.
159: */
160: { register short precision; /* number of units to shift */
161: register unit carry,nextcarry,bitmask;
162: register short unbits;
163: if (bits==0) return; /* shift zero bits is a no-op */
164: carry = 0;
165: bitmask = power_of_2(bits)-1;
166: unbits = UNITSIZE-bits; /* shift bits must be < UNITSIZE */
167: precision = global_precision;
168: make_msbptr(r1,precision);
169: while (precision--)
170: { nextcarry = *r1 & bitmask;
171: *r1 >>= bits;
172: *r1 |= carry << unbits;
173: carry = nextcarry;
174: pre_lowerunit(r1);
175: }
176: } /* mp_shift_right_bits */
177:
178:
179: #ifndef VMS
180:
181: short mp_compare(register unitptr r1,register unitptr r2)
182: /* Compares multiprecision integers *r1, *r2, and returns:
183: -1 iff *r1 < *r2
184: 0 iff *r1 == *r2
185: +1 iff *r1 > *r2
186: */
187: { register short precision; /* number of units to compare */
188:
189: precision = global_precision;
190: make_msbptr(r1,precision);
191: make_msbptr(r2,precision);
192: do
193: { if (*r1 < *r2)
194: return(-1);
195: if (*post_lowerunit(r1) > *post_lowerunit(r2))
196: return(1);
197: } while (--precision);
198: return(0); /* *r1 == *r2 */
199: } /* mp_compare */
200:
201: #endif /* VMS */
202:
203: boolean mp_inc(register unitptr r)
204: /* Increment multiprecision integer r. */
205: { register short precision;
206: precision = global_precision;
207: make_lsbptr(r,precision);
208: do
209: { if ( ++(*r) ) return(0); /* no carry */
210: post_higherunit(r);
211: } while (--precision);
212: return(1); /* return carry set */
213: } /* mp_inc */
214:
215:
216: boolean mp_dec(register unitptr r)
217: /* Decrement multiprecision integer r. */
218: { register short precision;
219: precision = global_precision;
220: make_lsbptr(r,precision);
221: do
222: { if ( (signedunit) (--(*r)) != (signedunit) -1 )
223: return(0); /* no borrow */
224: post_higherunit(r);
225: } while (--precision);
226: return(1); /* return borrow set */
227: } /* mp_dec */
228:
229:
230: void mp_neg(register unitptr r)
231: /* Compute 2's complement, the arithmetic negative, of r */
232: { register short precision; /* number of units to negate */
233: precision = global_precision;
234: mp_dec(r); /* 2's complement is 1's complement plus 1 */
235: do /* now do 1's complement */
236: { *r = ~(*r);
237: r++;
238: } while (--precision);
239: } /* mp_neg */
240:
241: #ifndef VAXC /* Replaced by a macro under VAX C */
242:
243: void mp_move(register unitptr dst,register unitptr src)
244: { register short precision; /* number of units to move */
245: precision = global_precision;
246: do { *dst++ = *src++; } while (--precision);
247: } /* mp_move */
248:
249: #endif /* VAXC */
250: void mp_init(register unitptr r, word16 value)
251: /* Init multiprecision register r with short value. */
252: { /* Note that mp_init doesn't extend sign bit for >32767 */
253: register short precision; /* number of units to init */
254: precision = global_precision;
255: make_lsbptr(r,precision);
256: *post_higherunit(r) = value; precision--;
257: #ifdef UNIT8
258: *post_higherunit(r) = value >> UNITSIZE; precision--;
259: #endif
260: #ifdef VAXC
261:
262: unitfill0( r, precision); /* Use character insts. on a VAX */
263:
264: #else /* VAXC */
265:
266: while (precision--)
267: *post_higherunit(r) = 0;
268:
269: #endif /* VAXC */
270:
271: } /* mp_init */
272:
273:
274: short significance(register unitptr r)
275: /* Returns number of significant units in r */
276: { register short precision;
277: precision = global_precision;
278: make_msbptr(r,precision);
279: do
280: { if (*post_lowerunit(r))
281: return(precision);
282: } while (--precision);
283: return(precision);
284: } /* significance */
285:
286:
287: #ifndef VAXC /* Replaced by a macro under VAX C */
288:
289: void unitfill0(unitptr r,word16 unitcount)
290: /* Zero-fill the unit buffer r. */
291: { while (unitcount--) *r++ = 0;
292: } /* unitfill0 */
293:
294: #endif /* VAXC */
295:
296: /* The macro normalize(r,precision) "normalizes" a multiprecision integer
297: by adjusting r and precision to new values. For Motorola-style processors
298: (MSB-first), r is a pointer to the MSB of the register, and must
299: be adjusted to point to the first nonzero unit. For Intel/VAX-style
300: (LSB-first) processors, r is a pointer to the LSB of the register,
301: and must be left unchanged. The precision counter is always adjusted,
302: regardless of processor type. In the case of precision = 0,
303: r becomes undefined.
304: */
305:
306:
307: /* The macro rescale(r,current_precision,new_precision) rescales
308: a multiprecision integer by adjusting r and its precision to new values.
309: It can be used to reverse the effects of the normalize
310: routine given above. See the comments in normalize concerning
311: Intel vs. Motorola LSB/MSB conventions.
312: WARNING: You can only safely call rescale on registers that
313: you have previously normalized with the above normalize routine,
314: or are known to be big enough for the new precision. You may
315: specify a new precision that is smaller than the current precision.
316: You must be careful not to specify a new_precision value that is
317: too big, or which adjusts the r pointer out of range.
318: */
319:
320:
321: /* The "bit sniffer" macros all begin sniffing at the MSB.
322: They are used internally by all the various multiply, divide,
323: modulo, exponentiation, and square root functions.
324: */
325:
326:
327: int mp_udiv(register unitptr remainder,register unitptr quotient,
328: register unitptr dividend,register unitptr divisor)
329: /* Unsigned divide, treats both operands as positive. */
330: { int bits;
331: short dprec;
332: register unit bitmask;
333: if (testeq(divisor,0))
334: return(-1); /* zero divisor means divide error */
335: mp_init0(remainder);
336: mp_init0(quotient);
337: /* normalize and compute number of bits in dividend first */
338: init_bitsniffer(dividend,bitmask,dprec,bits);
339: /* rescale quotient to same precision (dprec) as dividend */
340: rescale(quotient,global_precision,dprec);
341: make_msbptr(quotient,dprec);
342:
343: while (bits--)
344: { mp_rotate_left(remainder,(boolean)(sniff_bit(dividend,bitmask)!=0));
345: if (mp_compare(remainder,divisor) >= 0)
346: { mp_sub(remainder,divisor);
347: stuff_bit(quotient,bitmask);
348: }
349: bump_2bitsniffers(dividend,quotient,bitmask);
350: }
351: return(0);
352: } /* mp_udiv */
353:
354:
355: #define RECIPMARGIN 0 /* extra margin bits used by mp_recip() */
356:
357: int mp_recip(register unitptr quotient,register unitptr divisor)
358: /* Compute reciprocal (quotient) as 1/divisor. Used by faster modmult. */
359: { int bits;
360: short qprec;
361: register unit bitmask;
362: unit remainder[MAX_UNIT_PRECISION];
363: if (testeq(divisor,0))
364: return(-1); /* zero divisor means divide error */
365: mp_init0(remainder);
366: mp_init0(quotient);
367:
368: /* normalize and compute number of bits in quotient first */
369: bits = countbits(divisor) + RECIPMARGIN;
370: bitmask = bitmsk(bits); /* bitmask within a single unit */
371: qprec = bits2units(bits+1);
372: mp_setbit(remainder,(bits-RECIPMARGIN)-1);
373: /* rescale quotient to precision of divisor + RECIPMARGIN bits */
374: rescale(quotient,global_precision,qprec);
375: make_msbptr(quotient,qprec);
376:
377: while (bits--)
378: { mp_shift_left(remainder);
379: if (mp_compare(remainder,divisor) >= 0)
380: { mp_sub(remainder,divisor);
381: stuff_bit(quotient,bitmask);
382: }
383: bump_bitsniffer(quotient,bitmask);
384: }
385: mp_init0(remainder); /* burn sensitive data left on stack */
386: return(0);
387: } /* mp_recip */
388:
389:
390: int mp_div(register unitptr remainder,register unitptr quotient,
391: register unitptr dividend,register unitptr divisor)
392: /* Signed divide, either or both operands may be negative. */
393: { boolean dvdsign,dsign;
394: int status;
395: dvdsign = (mp_tstminus(dividend)!=0);
396: dsign = (mp_tstminus(divisor)!=0);
397: if (dvdsign) mp_neg(dividend);
398: if (dsign) mp_neg(divisor);
399: status = mp_udiv(remainder,quotient,dividend,divisor);
400: if (dvdsign) mp_neg(dividend); /* repair caller's dividend */
401: if (dsign) mp_neg(divisor); /* repair caller's divisor */
402: if (status<0) return(status); /* divide error? */
403: if (dvdsign) mp_neg(remainder);
404: if (dvdsign ^ dsign) mp_neg(quotient);
405: return(status);
406: } /* mp_div */
407:
408:
409: word16 mp_shortdiv(register unitptr quotient,
410: register unitptr dividend,register word16 divisor)
411: /* This function does a fast divide and mod on a multiprecision dividend
412: using a short integer divisor returning a short integer remainder.
413: This is an unsigned divide. It treats both operands as positive.
414: It is used mainly for faster printing of large numbers in base 10.
415: */
416: { int bits;
417: short dprec;
418: register unit bitmask;
419: register word16 remainder;
420: if (!divisor) /* if divisor == 0 */
421: return(-1); /* zero divisor means divide error */
422: remainder=0;
423: mp_init0(quotient);
424: /* normalize and compute number of bits in dividend first */
425: init_bitsniffer(dividend,bitmask,dprec,bits);
426: /* rescale quotient to same precision (dprec) as dividend */
427: rescale(quotient,global_precision,dprec);
428: make_msbptr(quotient,dprec);
429:
430: while (bits--)
431: { remainder <<= 1;
432: if (sniff_bit(dividend,bitmask))
433: remainder++;
434: if (remainder >= divisor)
435: { remainder -= divisor;
436: stuff_bit(quotient,bitmask);
437: }
438: bump_2bitsniffers(dividend,quotient,bitmask);
439: }
440: return(remainder);
441: } /* mp_shortdiv */
442:
443:
444: int mp_mod(register unitptr remainder,
445: register unitptr dividend,register unitptr divisor)
446: /* Unsigned divide, treats both operands as positive. */
447: { int bits;
448: short dprec;
449: register unit bitmask;
450: if (testeq(divisor,0))
451: return(-1); /* zero divisor means divide error */
452: mp_init0(remainder);
453: /* normalize and compute number of bits in dividend first */
454: init_bitsniffer(dividend,bitmask,dprec,bits);
455:
456: while (bits--)
457: { mp_rotate_left(remainder,(boolean)(sniff_bit(dividend,bitmask)!=0));
458: msub(remainder,divisor);
459: bump_bitsniffer(dividend,bitmask);
460: }
461: return(0);
462: } /* mp_mod */
463:
464:
465: word16 mp_shortmod(register unitptr dividend,register word16 divisor)
466: /* This function does a fast mod operation on a multprecision dividend
467: using a short integer modulus returning a short integer remainder.
468: This is an unsigned divide. It treats both operands as positive.
469: It is used mainly for fast sieve searches for large primes.
470: */
471: { int bits;
472: short dprec;
473: register unit bitmask;
474: register word16 remainder;
475: if (!divisor) /* if divisor == 0 */
476: return(-1); /* zero divisor means divide error */
477: remainder=0;
478: /* normalize and compute number of bits in dividend first */
479: init_bitsniffer(dividend,bitmask,dprec,bits);
480:
481: while (bits--)
482: { remainder <<= 1;
483: if (sniff_bit(dividend,bitmask))
484: remainder++;
485: if (remainder >= divisor) remainder -= divisor;
486: bump_bitsniffer(dividend,bitmask);
487: }
488: return(remainder);
489: } /* mp_shortmod */
490:
491:
492:
493: #ifdef COMB_MULT /* use faster "comb" multiply algorithm */
494: /* We are skipping this code because it has a bug... */
495:
496: int mp_mult(register unitptr prod,
497: register unitptr multiplicand, register unitptr multiplier)
498: /* Computes multiprecision prod = multiplicand * multiplier */
499: { /* Uses interleaved comb multiply algorithm.
500: This improved multiply more than twice as fast as a Russian
501: peasant multiply, because it does a lot fewer shifts.
502: Must have global_precision set to the size of the multiplicand
503: plus UNITSIZE-1 SLOP_BITS. Produces a product that is the sum
504: of the lengths of the multiplier and multiplicand.
505:
506: BUG ALERT: Unfortunately, this code has a bug. It fails for
507: some numbers. One such example:
508: x= 59DE 60CE 2345 8091 A02B 2A1C DBC3 8BE5
509: x*x= 59DE 60CE 2345 26B3 993B 67A5 2499 0B7D
510: 52C8 CDC7 AFB3 61C8 243C 741B
511: --which is obviously wrong. The answer should be:
512: x*x= 1F8C 607B 5EA6 C061 2714 04A9 A0C6 A17A
513: C9AB 6095 C62F 3756 3843 E4D0 3950 7AD9
514: We'll have to fix this some day. In the meantime, we'll
515: just have the compiler skip over this code.
516: */
517: int bits;
518: register unit bitmask;
519: unitptr product, mplier, temp;
520: short mprec,mprec2;
521: unit mplicand[MAX_UNIT_PRECISION];
522: mp_init(prod,0); /* better clear full width--double precision */
523: if (testeq(multiplicand,0))
524: return(0); /* zero multiplicand means zero product */
525:
526: mp_move(mplicand,multiplicand); /* save it from damage */
527:
528: normalize(multiplier,mprec);
529: if (!mprec)
530: return(0);
531:
532: make_lsbptr(multiplier,mprec);
533: bitmask = 1; /* start scan at LSB of multiplier */
534:
535: do /* UNITSIZE times */
536: { /* do for bits 0-15 */
537: product = prod;
538: mplier = multiplier;
539: mprec2 = mprec;
540: while (mprec2--) /* do for each word in multiplier */
541: {
542: if (sniff_bit(mplier,bitmask))
543: { if (mp_addc(product,multiplicand,0)) /* ripple carry */
544: { /* After 1st time thru, this is rarely encountered. */
545: temp = msbptr(product,global_precision);
546: pre_higherunit(temp);
547: /* temp now points to LSU of carry region. */
548: unmake_lsbptr(temp,global_precision);
549: mp_inc(temp);
550: } /* ripple carry */
551: }
552: pre_higherunit(mplier);
553: pre_higherunit(product);
554: }
555: if (!(bitmask <<= 1))
556: break;
557: mp_shift_left(multiplicand);
558:
559: } while (TRUE);
560:
561: mp_move(multiplicand,mplicand); /* recover */
562:
563: return(0); /* normal return */
564: } /* mp_mult */
565:
566: #endif /* COMB_MULT */
567:
568:
569: /* Because the "comb" multiply has a bug, we will use the slower
570: Russian peasant multiply instead. Fortunately, the mp_mult
571: function is not called from any time-critical code.
572: */
573:
574: int mp_mult(register unitptr prod,
575: register unitptr multiplicand,register unitptr multiplier)
576: /* Computes multiprecision prod = multiplicand * multiplier */
577: { /* Uses "Russian peasant" multiply algorithm. */
578: int bits;
579: register unit bitmask;
580: short mprec;
581: mp_init(prod,0);
582: if (testeq(multiplicand,0))
583: return(0); /* zero multiplicand means zero product */
584: /* normalize and compute number of bits in multiplier first */
585: init_bitsniffer(multiplier,bitmask,mprec,bits);
586:
587: while (bits--)
588: { mp_shift_left(prod);
589: if (sniff_bit(multiplier,bitmask))
590: mp_add(prod,multiplicand);
591: bump_bitsniffer(multiplier,bitmask);
592: }
593: return(0);
594: } /* mp_mult */
595:
596:
597:
598: /* mp_modmult computes a multiprecision multiply combined with a
599: modulo operation. This is the most time-critical function in
600: this multiprecision arithmetic library for performing modulo
601: exponentiation. We experimented with different versions of modmult,
602: depending on the machine architecture and performance requirements.
603: We will either use a Russian Peasant modmult (peasant_modmult),
604: Charlie Merritt's modmult (merritt_modmult) or Jimmy Upton's
605: modmult (upton_modmult). On machines with a hardware atomic
606: multiply instruction, Upton's modmult is fastest, which depends
607: on an assembly subroutine to speed up the hardware multiply logic.
608: If the machine lacks a fast hardware multiply, Merritt's modmult
609: is preferred, which doesn't call any assembly multiply routine.
610: We use the alias names mp_modmult, stage_modulus, and modmult_burn
611: for the corresponding true names, which depend on what flavor of
612: modmult we are using.
613:
614: Before making the first call to mp_modmult, you must set up the
615: modulus-dependant precomputated tables by calling stage_modulus.
616: After making all the calls to mp_modmult, you call modmult_burn to
617: erase the tables created by stage_modulus that were left in memory.
618: */
619:
620: #ifdef COUNTMULTS
621: /* "number of modmults" counters, used for performance studies. */
622: static unsigned int tally_modmults = 0;
623: static unsigned int tally_modsquares = 0;
624: #endif /* COUNTMULTS */
625:
626:
627: #ifdef PEASANT
628: /* Conventional Russian peasant multiply with modulo algorithm. */
629:
630: static unitptr modulus = 0; /* used only by mp_modmult */
631:
632: int stage_peasant_modulus(unitptr n)
633: /* Must pass modulus to stage_modulus before calling modmult.
634: Assumes that global_precision has already been adjusted to the
635: size of the modulus, plus SLOP_BITS.
636: */
637: { /* For this simple version of modmult, just copy unit pointer. */
638: modulus = n;
639: return(0); /* normal return */
640: } /* stage_peasant_modulus */
641:
642:
643: int peasant_modmult(register unitptr prod,
644: unitptr multiplicand,register unitptr multiplier)
645: { /* "Russian peasant" multiply algorithm, combined with a modulo
646: operation. This is a simple naive replacement for Merritt's
647: faster modmult algorithm. References global unitptr "modulus".
648: Computes: prod = (multiplicand*multiplier) mod modulus
649: WARNING: All the arguments must be less than the modulus!
650: */
651: int bits;
652: register unit bitmask;
653: short mprec;
654: mp_init(prod,0);
655: /* if (testeq(multiplicand,0))
656: return(0); */ /* zero multiplicand means zero product */
657: /* normalize and compute number of bits in multiplier first */
658: init_bitsniffer(multiplier,bitmask,mprec,bits);
659:
660: while (bits--)
661: { mp_shift_left(prod);
662: msub(prod,modulus); /* turns mult into modmult */
663: if (sniff_bit(multiplier,bitmask))
664: { mp_add(prod,multiplicand);
665: msub(prod,modulus); /* turns mult into modmult */
666: }
667: bump_bitsniffer(multiplier,bitmask);
668: }
669: return(0);
670: } /* peasant_modmult */
671:
672:
673: /* If we are using a version of mp_modmult that uses internal tables
674: in memory, we have to call modmult_burn() at the end of mp_modexp.
675: This is so that no sensitive data is left in memory after the program
676: exits. The Russian peasant method doesn't use any such tables.
677: */
678: static void peasant_burn(void)
679: /* Alias for modmult_burn, called only from mp_modexp(). Destroys
680: internal modmult tables. This version does nothing because no
681: tables are used by the Russian peasant modmult. */
682: { } /* peasant_burn */
683:
684: #endif /* PEASANT */
685:
686:
687: #ifdef MERRITT
688: /*=========================================================================*/
689: /*
690: This is Charlie Merritt's MODMULT algorithm, implemented in C by PRZ.
691: Also refined by Zhahai Stewart to reduce the number of subtracts.
692: It performs a multiply combined with a modulo operation, without
693: going into "double precision". It is faster than the Russian peasant
694: method, and still works well on machines that lack a fast hardware
695: multiply instruction.
696: */
697:
698: /* The following support functions, macros, and data structures
699: are used only by Merritt's modmult algorithm... */
700:
701: static void mp_lshift_unit(register unitptr r1)
702: /* Shift r1 1 whole unit to the left. Used only by modmult function. */
703: { register short precision;
704: register unitptr r2;
705: precision = global_precision;
706: make_msbptr(r1,precision);
707: r2 = r1;
708: while (--precision)
709: *post_lowerunit(r1) = *pre_lowerunit(r2);
710: *r1 = 0;
711: } /* mp_lshift_unit */
712:
713:
714: /* moduli_buf contains shifted images of the modulus, set by stage_modulus */
715: static unit moduli_buf[UNITSIZE][MAX_UNIT_PRECISION] = {0};
716: static unitptr moduli[UNITSIZE+1] = /* contains pointers into moduli_buf */
717: { 0
718: ,&moduli_buf[ 0][0], &moduli_buf[ 1][0], &moduli_buf[ 2][0], &moduli_buf[ 3][0],
719: &moduli_buf[ 4][0], &moduli_buf[ 5][0], &moduli_buf[ 6][0], &moduli_buf[ 7][0]
720: #ifndef UNIT8
721: ,&moduli_buf[ 8][0], &moduli_buf[ 9][0], &moduli_buf[10][0], &moduli_buf[11][0],
722: &moduli_buf[12][0], &moduli_buf[13][0], &moduli_buf[14][0], &moduli_buf[15][0]
723: #ifndef UNIT16 /* and not UNIT8 */
724: ,&moduli_buf[16][0], &moduli_buf[17][0], &moduli_buf[18][0], &moduli_buf[19][0],
725: &moduli_buf[20][0], &moduli_buf[21][0], &moduli_buf[22][0], &moduli_buf[23][0],
726: &moduli_buf[24][0], &moduli_buf[25][0], &moduli_buf[26][0], &moduli_buf[27][0],
727: &moduli_buf[28][0], &moduli_buf[29][0], &moduli_buf[30][0], &moduli_buf[31][0]
728: #endif /* UNIT16 and UNIT8 not defined */
729: #endif /* UNIT8 not defined */
730: };
731:
732: /* To optimize msubs, need following 2 unit arrays, each filled
733: with the most significant unit and next-to-most significant unit
734: of the preshifted images of the modulus. */
735: static unit msu_moduli[UNITSIZE+1] = {0}; /* most signif. unit */
736: static unit nmsu_moduli[UNITSIZE+1] = {0}; /* next-most signif. unit */
737:
738: /* mpdbuf contains preshifted images of the multiplicand, mod n.
739: It is used only by mp_modmult. It could be staticly declared
740: inside of mp_modmult, but we put it outside mp_modmult so that
741: it can be wiped clean by modmult_burn(), which is called at the
742: end of mp_modexp. This is so that no sensitive data is left in
743: memory after the program exits.
744: */
745: static unit mpdbuf[UNITSIZE-1][MAX_UNIT_PRECISION] = {0};
746:
747:
748: static void stage_mp_images(unitptr images[UNITSIZE],unitptr r)
749: /* Computes UNITSIZE images of r, each shifted left 1 more bit.
750: Used only by modmult function.
751: */
752: { short int i;
753: images[0] = r; /* no need to move the first image, just copy ptr */
754: for (i=1; i<UNITSIZE; i++)
755: { mp_move(images[i],images[i-1]);
756: mp_shift_left(images[i]);
757: }
758: } /* stage_mp_images */
759:
760:
761: int stage_merritt_modulus(unitptr n)
762: /* Computes UNITSIZE+1 images of modulus, each shifted left 1 more bit.
763: Before calling mp_modmult, you must first stage the moduli images by
764: calling stage_modulus. n is the pointer to the modulus.
765: Assumes that global_precision has already been adjusted to the
766: size of the modulus, plus SLOP_BITS.
767: */
768: { short int i;
769: unitptr msu; /* ptr to most significant unit, for faster msubs */
770: moduli[0] = n; /* no need to move the first image, just copy ptr */
771:
772: /* used by optimized msubs macro... */
773: msu = msbptr(n,global_precision); /* needed by msubs */
774: msu_moduli[0] = *post_lowerunit(msu);
775: nmsu_moduli[0] = *msu;
776:
777: for (i=1; i<UNITSIZE+1; i++)
778: { mp_move(moduli[i],moduli[i-1]);
779: mp_shift_left(moduli[i]);
780:
781: /* used by optimized msubs macro... */
782: msu = msbptr(moduli[i],global_precision); /* needed by msubs */
783: msu_moduli[i] = *post_lowerunit(msu); /* for faster msubs */
784: nmsu_moduli[i] = *msu;
785: }
786: return(0); /* normal return */
787: } /* stage_merritt_modulus */
788:
789:
790: /* The following macros, sniffadd and msubs, are used by modmult... */
791: #define sniffadd(i) if (*multiplier & power_of_2(i)) mp_add(prod,mpd[i])
792:
793: /* Unoptimized msubs macro (msubs0) follows... */
794: /* #define msubs0(i) msub(prod,moduli[i])
795: */
796:
797: /* To optimize msubs, compare the most significant units of the
798: product and the shifted modulus before deciding to call the full
799: mp_compare routine. Better still, compare the upper two units
800: before deciding to call mp_compare.
801: Optimization of msubs relies heavily on standard C left-to-right
802: parsimonious evaluation of logical expressions.
803: */
804:
805: /* Partially-optimized msubs macro (msubs1) follows... */
806: /* #define msubs1(i) if ( \
807: ((p_m = (*msu_prod-msu_moduli[i])) >= 0) && \
808: (p_m || (mp_compare(prod,moduli[i]) >= 0) ) \
809: ) mp_sub(prod,moduli[i])
810: */
811:
812: /* Fully-optimized msubs macro (msubs2) follows... */
813: #define msubs(i) if (((p_m = *msu_prod-msu_moduli[i])>0) || ( \
814: (p_m==0) && ( (*nmsu_prod>nmsu_moduli[i]) || ( \
815: (*nmsu_prod==nmsu_moduli[i]) && ((mp_compare(prod,moduli[i]) >= 0)) ))) ) \
816: mp_sub(prod,moduli[i])
817:
818:
819: int merritt_modmult(register unitptr prod,
820: unitptr multiplicand,register unitptr multiplier)
821: /* Performs combined multiply/modulo operation.
822: Computes: prod = (multiplicand*multiplier) mod modulus
823: WARNING: All the arguments must be less than the modulus!
824: Assumes the modulus has been predefined by first calling
825: stage_modulus.
826: */
827: {
828: /* p_m, msu_prod, and nmsu_prod are used by the optimized msubs macro...*/
829: register signedunit p_m;
830: register unitptr msu_prod; /* ptr to most significant unit of product */
831: register unitptr nmsu_prod; /* next-most signif. unit of product */
832: short mprec; /* precision of multiplier, in units */
833: /* Array mpd contains a list of pointers to preshifted images of
834: the multiplicand: */
835: static unitptr mpd[UNITSIZE] =
836: { 0, &mpdbuf[ 0][0], &mpdbuf[ 1][0], &mpdbuf[ 2][0],
837: &mpdbuf[ 3][0], &mpdbuf[ 4][0], &mpdbuf[ 5][0], &mpdbuf[ 6][0]
838: #ifndef UNIT8
839: ,&mpdbuf[ 7][0], &mpdbuf[ 8][0], &mpdbuf[ 9][0], &mpdbuf[10][0],
840: &mpdbuf[11][0], &mpdbuf[12][0], &mpdbuf[13][0], &mpdbuf[14][0]
841: #ifndef UNIT16 /* and not UNIT8 */
842: ,&mpdbuf[15][0], &mpdbuf[16][0], &mpdbuf[17][0], &mpdbuf[18][0],
843: &mpdbuf[19][0], &mpdbuf[20][0], &mpdbuf[21][0], &mpdbuf[22][0],
844: &mpdbuf[23][0], &mpdbuf[24][0], &mpdbuf[25][0], &mpdbuf[26][0],
845: &mpdbuf[27][0], &mpdbuf[28][0], &mpdbuf[29][0], &mpdbuf[30][0]
846: #endif /* UNIT16 and UNIT8 not defined */
847: #endif /* UNIT8 not defined */
848: };
849:
850: /* Compute preshifted images of multiplicand, mod n: */
851: stage_mp_images(mpd,multiplicand);
852:
853: /* To optimize msubs, set up msu_prod and nmsu_prod: */
854: msu_prod = msbptr(prod,global_precision); /* Get ptr to MSU of prod */
855: nmsu_prod = msu_prod;
856: post_lowerunit(nmsu_prod); /* Get next-MSU of prod */
857:
858: /* To understand this algorithm, it would be helpful to first
859: study the conventional Russian peasant modmult algorithm.
860: This one does about the same thing as Russian peasant, but
861: is organized differently to save some steps. It loops
862: through the multiplier a word (unit) at a time, instead of
863: a bit at a time. It word-shifts the product instead of
864: bit-shifting it, so it should be faster. It also does about
865: half as many subtracts as Russian peasant.
866: */
867:
868: mp_init(prod,0); /* Initialize product to 0. */
869:
870: /* The way mp_modmult is actually used in cryptographic
871: applications, there will NEVER be a zero multiplier or
872: multiplicand. So there is no need to optimize for that
873: condition.
874: */
875: /* if (testeq(multiplicand,0))
876: return(0); */ /* zero multiplicand means zero product */
877: /* Normalize and compute number of units in multiplier first: */
878: normalize(multiplier,mprec);
879: if (mprec==0) /* if precision of multiplier is 0 */
880: return(0); /* zero multiplier means zero product */
881: make_msbptr(multiplier,mprec); /* start at MSU of multiplier */
882:
883: while (mprec--) /* Loop for the number of units in the multiplier */
884: {
885: /* Shift the product one whole unit to the left.
886: This is part of the multiply phase of modmult.
887: */
888:
889: mp_lshift_unit(prod);
890:
891: /* Sniff each bit in the current unit of the multiplier,
892: and conditionally add the the corresponding preshifted
893: image of the multiplicand to the product.
894: This is also part of the multiply phase of modmult.
895:
896: The following loop is unrolled for speed, using macros...
897:
898: for (i=UNITSIZE-1; i>=0; i--)
899: if (*multiplier & power_of_2(i))
900: mp_add(prod,mpd[i]);
901: */
902: #ifndef UNIT8
903: #ifndef UNIT16 /* and not UNIT8 */
904: sniffadd(31);
905: sniffadd(30);
906: sniffadd(29);
907: sniffadd(28);
908: sniffadd(27);
909: sniffadd(26);
910: sniffadd(25);
911: sniffadd(24);
912: sniffadd(23);
913: sniffadd(22);
914: sniffadd(21);
915: sniffadd(20);
916: sniffadd(19);
917: sniffadd(18);
918: sniffadd(17);
919: sniffadd(16);
920: #endif /* not UNIT16 and not UNIT8 */
921: sniffadd(15);
922: sniffadd(14);
923: sniffadd(13);
924: sniffadd(12);
925: sniffadd(11);
926: sniffadd(10);
927: sniffadd(9);
928: sniffadd(8);
929: #endif /* not UNIT8 */
930: sniffadd(7);
931: sniffadd(6);
932: sniffadd(5);
933: sniffadd(4);
934: sniffadd(3);
935: sniffadd(2);
936: sniffadd(1);
937: sniffadd(0);
938:
939: /* The product may have grown by as many as UNITSIZE+1
940: bits. That's why we have global_precision set to the
941: size of the modulus plus UNITSIZE+1 slop bits.
942: Now reduce the product back down by conditionally
943: subtracting the UNITSIZE+1 preshifted images of the
944: modulus. This is the modulo reduction phase of modmult.
945:
946: The following loop is unrolled for speed, using macros...
947:
948: for (i=UNITSIZE; i>=0; i--)
949: if (mp_compare(prod,moduli[i]) >= 0)
950: mp_sub(prod,moduli[i]);
951: */
952: #ifndef UNIT8
953: #ifndef UNIT16 /* and not UNIT8 */
954: msubs(32);
955: msubs(31);
956: msubs(30);
957: msubs(29);
958: msubs(28);
959: msubs(27);
960: msubs(26);
961: msubs(25);
962: msubs(24);
963: msubs(23);
964: msubs(22);
965: msubs(21);
966: msubs(20);
967: msubs(19);
968: msubs(18);
969: msubs(17);
970: #endif /* not UNIT16 and not UNIT8 */
971: msubs(16);
972: msubs(15);
973: msubs(14);
974: msubs(13);
975: msubs(12);
976: msubs(11);
977: msubs(10);
978: msubs(9);
979: #endif /* not UNIT8 */
980: msubs(8);
981: msubs(7);
982: msubs(6);
983: msubs(5);
984: msubs(4);
985: msubs(3);
986: msubs(2);
987: msubs(1);
988: msubs(0);
989:
990: /* Bump pointer to next lower unit of multiplier: */
991: post_lowerunit(multiplier);
992:
993: } /* Loop for the number of units in the multiplier */
994:
995: return(0); /* normal return */
996:
997: } /* merritt_modmult */
998:
999:
1000: #undef msubs
1001: #undef sniffadd
1002:
1003:
1004: /* Merritt's mp_modmult function leaves some internal tables in memory,
1005: so we have to call modmult_burn() at the end of mp_modexp.
1006: This is so that no cryptographically sensitive data is left in memory
1007: after the program exits.
1008: */
1009: static void merritt_burn(void)
1010: /* Alias for modmult_burn, merritt_burn() is called only by mp_modexp. */
1011: { unitfill0(&(mpdbuf[0][0]),(UNITSIZE-1)*MAX_UNIT_PRECISION);
1012: unitfill0(&(moduli_buf[0][0]),(UNITSIZE)*MAX_UNIT_PRECISION);
1013: unitfill0(msu_moduli,UNITSIZE+1);
1014: unitfill0(nmsu_moduli,UNITSIZE+1);
1015: } /* merritt_burn() */
1016:
1017: /******* end of Merritt's MODMULT stuff. *******/
1018: /*=========================================================================*/
1019: #endif /* MERRITT */
1020:
1021:
1022: #ifdef UPTON /* using Jimmy Upton's modmult algorithm */
1023:
1024: /* Jimmy Upton's multiprecision modmult algorithm in C.
1025: Performs a multiply combined with a modulo operation.
1026:
1027: The following support functions and data structures
1028: are used only by Upton's modmult algorithm...
1029: */
1030:
1031: /* The MULTUNIT word is the biggest word size for the native atomic
1032: multiply instruction. It may or may not be smaller than UNITSIZE.
1033: Many CPU's have 16-by-16-bit multiplies, yielding a 32-bit product.
1034: In that case, MULTUINIT should be a 16-bit word, even if the rest of
1035: the multiprecision arithmetic functions assume a 32-bit UNIT word.
1036: */
1037: typedef unsigned short MULTUNIT;
1038: #define MULTUNITSIZE (8*sizeof(MULTUNIT)) /* size in bits */
1039:
1040: static short unit_prec; /* global_precision expressed in MULTUNITs */
1041:
1042:
1043: /* Define MPORTABLE if there is no assembly version of the mp_smul
1044: function available. An assembly version is much faster.
1045: Note that since the SPARC CPU has no hardware integer multiply
1046: instruction, there is not that much advantage in having an
1047: assembly version of mp_smul on that machine. It might be faster
1048: to use Merritt's modmult instead of Upton's modmult on the SPARC.
1049: */
1050: #ifdef MSDOS /* we do have an assembly version available on the 8086 */
1051: #undef MPORTABLE
1052: #endif
1053:
1054: #ifdef MPORTABLE /* use slow portable C version instead of assembly */
1055:
1056: /*
1057: Multiply the single-word multiplier times the multiprecision integer
1058: in multiplicand, accumulating result in prod. The resulting
1059: multiprecision prod will be 1 word longer than the multiplicand.
1060: multiplicand is unit_prec words long. We add into prod, so caller
1061: should zero it out first. For best results, this time-critical
1062: function should be implemented in assembly.
1063: NOTE: Unlike other functions in the multiprecision arithmetic
1064: library, both multiplicand and prod are pointing at the LSB,
1065: regardless of byte order of the machine. On an 80x86, this makes
1066: no difference. But if this assembly function is implemented
1067: on a 680x0, it becomes important.
1068: */
1069: void mp_smul (MULTUNIT *prod, MULTUNIT *multiplicand, MULTUNIT multiplier)
1070: {
1071: short i;
1072: unsigned long p, carry;
1073:
1074: carry = 0;
1075: for (i=0; i<unit_prec; ++i)
1076: { p = (unsigned long)multiplier * *post_higherunit(multiplicand);
1077: p += *prod + carry;
1078: *post_higherunit(prod) = (MULTUNIT)p;
1079: carry = p >> MULTUNITSIZE;
1080: }
1081: /* We know that the high-order word of prod will always be 0 */
1082: *prod = (MULTUNIT)carry; /* copy carry to empty high word of prod */
1083: } /* mp_smul */
1084:
1085: #else /* not MPORTABLE-- use assembly routine */
1086:
1087: #define mp_smul P_SMUL /* use external assembly routine */
1088:
1089: #endif /* MPORTABLE */
1090:
1091: #ifdef VMS
1092:
1093: #define mp_dmul p_dmul /* use external assembly routine */
1094:
1095: void mp_dmul (unitptr prod, unitptr multiplicand, unitptr multiplier);
1096:
1097: #else /* VMS */
1098:
1099: /* mp_dmul is a double-precision multiply multiplicand times multiplier,
1100: result into prod. prod must be pointing at a "double precision"
1101: buffer. E.g. If global_precision is 10 words, prod must be
1102: pointing at a 20-word buffer.
1103: */
1104: static void mp_dmul (unitptr prod, unitptr multiplicand, unitptr multiplier)
1105: {
1106: register int i;
1107: register MULTUNIT *p_multiplicand, *p_multiplier;
1108: register MULTUNIT *prodp;
1109:
1110: unitfill0(prod,global_precision*2); /* Pre-zero prod */
1111: /* Calculate precision in units of MULTUNIT */
1112: unit_prec = global_precision * UNITSIZE / MULTUNITSIZE;
1113: p_multiplicand = (MULTUNIT *)multiplicand;
1114: p_multiplier = (MULTUNIT *)multiplier;
1115: prodp = (MULTUNIT *)prod;
1116: make_lsbptr(p_multiplicand,unit_prec);
1117: make_lsbptr(p_multiplier,unit_prec);
1118: make_lsbptr(prodp,unit_prec*2);
1119: /* Multiply multiplicand by each word in multiplier, accumulating prod: */
1120: for (i=0; i<unit_prec; ++i)
1121: mp_smul (post_higherunit(prodp),
1122: p_multiplicand, *post_higherunit(p_multiplier));
1123: } /* mp_dmul */
1124:
1125: #endif /* VMS */
1126:
1127: /* These scratchpad arrays are used only by upton_modmult (mp_modmult).
1128: Some of them could be staticly declared inside of mp_modmult, but we
1129: put them outside mp_modmult so that they can be wiped clean by
1130: modmult_burn(), which is called at the end of mp_modexp. This is
1131: so that no sensitive data is left in memory after the program exits.
1132: */
1133:
1134: #ifdef VAXC
1135: /*
1136: * Operations on arrays are sped up on a VAX if the data is quadword (64-bit)
1137: * aligned. This is because many VAXen use a minimum of 64-bit data paths to
1138: * cache and memory. This technique is VAX C private, hence the conditional.
1139: */
1140: static unit _align( quadword ) modulus[MAX_UNIT_PRECISION];
1141: static unit _align( quadword ) reciprocal[MAX_UNIT_PRECISION];
1142: static unit _align( quadword ) dhi[MAX_UNIT_PRECISION];
1143: static unit _align( quadword ) d_data[MAX_UNIT_PRECISION*2]; /* double-wide register */
1144: static unit _align( quadword ) e_data[MAX_UNIT_PRECISION*2]; /* double-wide register */
1145: static unit _align( quadword ) f_data[MAX_UNIT_PRECISION*2]; /* double-wide register */
1146:
1147: #else /* VAXC */
1148:
1149: static unit modulus[MAX_UNIT_PRECISION];
1150: static unit reciprocal[MAX_UNIT_PRECISION];
1151: static unit dhi[MAX_UNIT_PRECISION];
1152: static unit d_data[MAX_UNIT_PRECISION*2]; /* double-wide register */
1153: static unit e_data[MAX_UNIT_PRECISION*2]; /* double-wide register */
1154: static unit f_data[MAX_UNIT_PRECISION*2]; /* double-wide register */
1155:
1156: #endif /* VAXC */
1157:
1158: static short nbits;
1159: static short nbitsDivUNITSIZE;
1160: static short nbitsModUNITSIZE;
1161:
1162: /* stage_upton_modulus() is aliased to stage_modulus().
1163: Prepare for an Upton modmult. Calculate the reciprocal of modulus,
1164: and save both. Note that reciprocal will have one more bit than
1165: modulus.
1166: Assumes that global_precision has already been adjusted to the
1167: size of the modulus, plus SLOP_BITS.
1168: */
1169: int stage_upton_modulus(unitptr n)
1170: {
1171: mp_move(modulus, n);
1172: mp_recip(reciprocal, modulus);
1173: nbits = countbits(modulus);
1174: nbitsDivUNITSIZE = nbits / UNITSIZE;
1175: nbitsModUNITSIZE = nbits % UNITSIZE;
1176: return(0); /* normal return */
1177: } /* stage_upton_modulus */
1178:
1179:
1180: /* Upton's algorithm performs a multiply combined with a modulo operation.
1181: Computes: prod = (multiplicand*multiplier) mod modulus
1182: WARNING: All the arguments must be less than the modulus!
1183: References global unitptr modulus and reciprocal.
1184: The reciprocal of modulus is 1 bit longer than the modulus.
1185: upton_modmult() is aliased to mp_modmult().
1186: */
1187: int upton_modmult (unitptr prod, unitptr multiplicand, unitptr multiplier)
1188: {
1189: unitptr d = d_data;
1190: unitptr d1 = d_data;
1191: unitptr e = e_data;
1192: unitptr f = f_data;
1193: short orig_precision, dprec;
1194: short i;
1195:
1196: orig_precision = global_precision;
1197: mp_dmul (d,multiplicand,multiplier);
1198:
1199: /* Throw off low nbits of d */
1200: #ifndef HIGHFIRST
1201: d1 = d + nbitsDivUNITSIZE;
1202: #else
1203: d1 = d + orig_precision - nbitsDivUNITSIZE;
1204: #endif
1205: mp_move(dhi, d1); /* Don't screw up d, we need it later */
1206: mp_shift_right_bits(dhi, nbitsModUNITSIZE);
1207:
1208: mp_dmul (e,dhi,reciprocal); /* Note - reciprocal has nbits+1 bits */
1209:
1210: #ifndef HIGHFIRST
1211: e += nbitsDivUNITSIZE;
1212: #else
1213: e += orig_precision - nbitsDivUNITSIZE;
1214: #endif
1215: mp_shift_right_bits(e, nbitsModUNITSIZE);
1216:
1217: mp_dmul (f,e,modulus);
1218:
1219: /* Now for the only double-precision call to mpilib: */
1220: set_precision (orig_precision * 2);
1221: mp_sub (d,f);
1222:
1223: /* d's precision should be <= orig_precision */
1224: rescale (d, orig_precision*2, orig_precision);
1225: set_precision (orig_precision);
1226:
1227: /* Should never have to do this final subtract more than twice: */
1228: while (mp_compare(d,modulus) > 0)
1229: mp_sub (d,modulus);
1230:
1231: mp_move(prod,d);
1232: return(0); /* normal return */
1233: } /* upton_modmult */
1234:
1235:
1236: /* Upton's mp_modmult function leaves some internal arrays in memory,
1237: so we have to call modmult_burn() at the end of mp_modexp.
1238: This is so that no cryptographically sensitive data is left in memory
1239: after the program exits.
1240: upton_burn() is aliased to modmult_burn().
1241: */
1242: void upton_burn(void)
1243: {
1244: unitfill0(modulus,MAX_UNIT_PRECISION);
1245: unitfill0(reciprocal,MAX_UNIT_PRECISION);
1246: unitfill0(dhi,MAX_UNIT_PRECISION);
1247: unitfill0(d_data,MAX_UNIT_PRECISION*2);
1248: unitfill0(e_data,MAX_UNIT_PRECISION*2);
1249: unitfill0(f_data,MAX_UNIT_PRECISION*2);
1250: nbits = nbitsDivUNITSIZE = nbitsModUNITSIZE = 0;
1251: } /* upton_burn */
1252:
1253: #endif /* UPTON */
1254:
1255:
1256: int countbits(unitptr r)
1257: /* Returns number of significant bits in r */
1258: { int bits;
1259: short prec;
1260: register unit bitmask;
1261: init_bitsniffer(r,bitmask,prec,bits);
1262: return(bits);
1263: } /* countbits */
1264:
1265:
1266: char *copyright_notice(void)
1267: /* force linker to include copyright notice in the executable object image. */
1268: { return ("(c)1986 Philip Zimmermann"); } /* copyright_notice */
1269:
1270:
1271: int mp_modexp(register unitptr expout,register unitptr expin,
1272: register unitptr exponent,register unitptr modulus)
1273: { /* Russian peasant combined exponentiation/modulo algorithm.
1274: Calls modmult instead of mult.
1275: Computes: expout = (expin**exponent) mod modulus
1276: WARNING: All the arguments must be less than the modulus!
1277: */
1278: int bits;
1279: short oldprecision;
1280: register unit bitmask;
1281: unit product[MAX_UNIT_PRECISION];
1282: short eprec;
1283:
1284: #ifdef COUNTMULTS
1285: tally_modmults = 0; /* clear "number of modmults" counter */
1286: tally_modsquares = 0; /* clear "number of modsquares" counter */
1287: #endif /* COUNTMULTS */
1288: mp_init(expout,1);
1289: if (testeq(exponent,0))
1290: { if (testeq(expin,0))
1291: return(-1); /* 0 to the 0th power means return error */
1292: return(0); /* otherwise, zero exponent means expout is 1 */
1293: }
1294: if (testeq(modulus,0))
1295: return(-2); /* zero modulus means error */
1296: #if SLOP_BITS > 0 /* if there's room for sign bits */
1297: if (mp_tstminus(modulus))
1298: return(-2); /* negative modulus means error */
1299: #endif /* SLOP_BITS > 0 */
1300: if (mp_compare(expin,modulus) >= 0)
1301: return(-3); /* if expin >= modulus, return error */
1302: if (mp_compare(exponent,modulus) >= 0)
1303: return(-4); /* if exponent >= modulus, return error */
1304:
1305: oldprecision = global_precision; /* save global_precision */
1306: /* set smallest optimum precision for this modulus */
1307: set_precision(bits2units(countbits(modulus)+SLOP_BITS));
1308: /* rescale all these registers to global_precision we just defined */
1309: rescale(modulus,oldprecision,global_precision);
1310: rescale(expin,oldprecision,global_precision);
1311: rescale(exponent,oldprecision,global_precision);
1312: rescale(expout,oldprecision,global_precision);
1313:
1314: if (stage_modulus(modulus))
1315: { set_precision(oldprecision); /* restore original precision */
1316: return(-5); /* unstageable modulus (STEWART algorithm) */
1317: }
1318:
1319: /* normalize and compute number of bits in exponent first */
1320: init_bitsniffer(exponent,bitmask,eprec,bits);
1321:
1322: /* We can "optimize out" the first modsquare and modmult: */
1323: bits--; /* We know for sure at this point that bits>0 */
1324: mp_move(expout,expin); /* expout = (1*1)*expin; */
1325: bump_bitsniffer(exponent,bitmask);
1326:
1327: while (bits--)
1328: {
1329: poll_for_break(); /* polls keyboard, allows ctrl-C to abort program */
1330: #ifdef COUNTMULTS
1331: tally_modsquares++; /* bump "number of modsquares" counter */
1332: #endif /* COUNTMULTS */
1333: mp_modsquare(product,expout);
1334: mp_move(expout,product);
1335: if (sniff_bit(exponent,bitmask))
1336: { mp_modmult(product,expout,expin);
1337: mp_move(expout,product);
1338: #ifdef COUNTMULTS
1339: tally_modmults++; /* bump "number of modmults" counter */
1340: #endif /* COUNTMULTS */
1341: }
1342: bump_bitsniffer(exponent,bitmask);
1343: } /* while bits-- */
1344: mp_burn(product); /* burn the evidence on the stack */
1345: modmult_burn(); /* ask mp_modmult to also burn its own evidence */
1346:
1347: #ifdef COUNTMULTS /* diagnostic analysis */
1348: { long atomic_mults;
1349: unsigned int unitcount,totalmults;
1350: unitcount = bits2units(countbits(modulus));
1351: /* calculation assumes modsquare takes as long as a modmult: */
1352: atomic_mults = (long) tally_modmults * (unitcount * unitcount);
1353: atomic_mults += (long) tally_modsquares * (unitcount * unitcount);
1354: printf("%ld atomic mults for ",atomic_mults);
1355: printf("%d+%d = %d modsqr+modmlt, at %d bits, %d words.\n",
1356: tally_modsquares,tally_modmults,
1357: tally_modsquares+tally_modmults,
1358: countbits(modulus), unitcount);
1359: }
1360: #endif /* COUNTMULTS */
1361:
1362: set_precision(oldprecision); /* restore original precision */
1363:
1364: /* Do an explicit reference to the copyright notice so that the linker
1365: will be forced to include it in the executable object image... */
1366: copyright_notice(); /* has no real effect at run time */
1367:
1368: return(0); /* normal return */
1369: } /* mp_modexp */
1370:
1371:
1372:
1373:
1374: /*********************************************************************
1375: RSA-specific routines follow. These are the only functions that
1376: are specific to the RSA public key cryptosystem. The other
1377: multiprecision integer math functions may be used for non-RSA
1378: applications.
1379:
1380: The RSA public key cryptosystem is patented by the Massachusetts
1381: Institute of Technology (U.S. patent #4,405,829). This patent does
1382: not apply outside the USA. Public Key Partners (PKP) holds the
1383: exclusive commercial license to sell and sub-license the RSA public
1384: key cryptosystem. The author of this software implementation of the
1385: RSA algorithm is providing this software for educational use only.
1386: Licensing this algorithm from PKP is the responsibility of you, the
1387: user, not Philip Zimmermann, the author of this software. The author
1388: assumes no liability for any breach of patent law resulting from the
1389: unlicensed use of this software by the user.
1390: *********************************************************************/
1391:
1392:
1393: int rsa_decrypt(unitptr M, unitptr C,
1394: unitptr d, unitptr p, unitptr q, unitptr u)
1395: /* Uses Quisquater's Chinese Remainder Theorem shortcut
1396: for RSA decryption.
1397: M is the output plaintext message.
1398: C is the input ciphertext message.
1399: d is the secret decryption exponent.
1400: p and q are the secret prime factors of n.
1401: u is the multiplicative inverse of p, mod q.
1402: Note that u is precomputed on the assumption that p<q.
1403: n, the common modulus, is not used here because of the
1404: Chinese Remainder Theorem shortcut.
1405: */
1406: {
1407: unit p2[MAX_UNIT_PRECISION];
1408: unit q2[MAX_UNIT_PRECISION];
1409: unit temp1[MAX_UNIT_PRECISION];
1410: unit temp2[MAX_UNIT_PRECISION];
1411: int status;
1412:
1413: mp_init(M,1); /* initialize result in case of error */
1414:
1415: if (mp_compare(p,q) >= 0) /* ensure that p<q */
1416: { /* swap the pointers p and q */
1417: unitptr t;
1418: t = p; p = q; q = t;
1419: }
1420:
1421: /* Rather than decrypting by computing modexp with full mod n
1422: precision, compute a shorter modexp with mod p precision... */
1423:
1424: /* p2 = [ (C mod p)**( d mod (p-1) ) ] mod p */
1425: mp_move(temp1,p);
1426: mp_dec(temp1); /* temp1 = p-1 */
1427: mp_mod(temp2,d,temp1); /* temp2 = d mod (p-1) ) */
1428: mp_mod(temp1,C,p); /* temp1 = C mod p */
1429: status = mp_modexp(p2,temp1,temp2,p);
1430: if (status < 0) /* mp_modexp returned an error. */
1431: return(status); /* error return */
1432:
1433: /* Now compute a short modexp with mod q precision... */
1434:
1435: /* q2 = [ (C mod q)**( d mod (q-1) ) ] mod q */
1436: mp_move(temp1,q);
1437: mp_dec(temp1); /* temp1 = q-1 */
1438: mp_mod(temp2,d,temp1); /* temp2 = d mod (q-1) */
1439: mp_mod(temp1,C,q); /* temp1 = C mod q */
1440: status = mp_modexp(q2,temp1,temp2,q);
1441: if (status < 0) /* mp_modexp returned an error. */
1442: return(status); /* error return */
1443:
1444: /* Now use the multiplicative inverse u to glue together the
1445: two halves, saving a lot of time by avoiding a full mod n
1446: exponentiation. At key generation time, u was computed
1447: with the secret key as the multiplicative inverse of
1448: p, mod q such that: (p*u) mod q = 1, assuming p<q.
1449: */
1450: if (mp_compare(p2,q2) == 0) /* only happens if C<p */
1451: mp_move(M,p2);
1452: else
1453: { /* q2 = q2 - p2; if q2<0 then q2 = q2 + q */
1454: if (mp_sub(q2,p2)) /* if the result went negative... */
1455: mp_add(q2,q); /* add q to q2 */
1456:
1457: /* M = p2 + ( p * [(q2*u) mod q] ) */
1458: mp_mult(temp1,q2,u); /* temp1 = q2*u */
1459: mp_mod(temp2,temp1,q); /* temp2 = temp1 mod q */
1460: mp_mult(temp1,p,temp2); /* temp1 = p * temp2 */
1461: mp_add(temp1,p2); /* temp1 = temp1 + p2 */
1462: mp_move(M,temp1); /* M = temp1 */
1463: }
1464:
1465: mp_burn(p2); /* burn the evidence on the stack...*/
1466: mp_burn(q2);
1467: mp_burn(temp1);
1468: mp_burn(temp2);
1469: return(0); /* normal return */
1470: } /* rsa_decrypt */
1471:
1472:
1473: /****************** end of MPI library ****************************/
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.