2 * Bignum routines for RSA and DH and stuff.
9 #if 0 // use PuTTY main debugging for diagbn()
12 #define debugprint debug
14 #define debugprint(x) printf x
17 #define BIGNUM_INTERNAL
18 typedef unsigned short *Bignum;
22 unsigned short bnZero[1] = { 0 };
23 unsigned short bnOne[2] = { 1, 1 };
26 * The Bignum format is an array of `unsigned short'. The first
27 * element of the array counts the remaining elements. The
28 * remaining elements express the actual number, base 2^16, _least_
29 * significant digit first. (So it's trivial to extract the bit
30 * with value 2^n for any n.)
32 * All Bignums in this module are positive. Negative numbers must
33 * be dealt with outside it.
35 * INVARIANT: the most significant word of any Bignum must be
39 Bignum Zero = bnZero, One = bnOne;
41 static Bignum newbn(int length) {
42 Bignum b = smalloc((length+1)*sizeof(unsigned short));
45 memset(b, 0, (length+1)*sizeof(*b));
50 void bn_restore_invariant(Bignum b) {
51 while (b[0] > 1 && b[b[0]] == 0) b[0]--;
54 Bignum copybn(Bignum orig) {
55 Bignum b = smalloc((orig[0]+1)*sizeof(unsigned short));
58 memcpy(b, orig, (orig[0]+1)*sizeof(*b));
62 void freebn(Bignum b) {
64 * Burn the evidence, just in case.
66 memset(b, 0, sizeof(b[0]) * (b[0] + 1));
70 Bignum bn_power_2(int n) {
71 Bignum ret = newbn(n/16+1);
72 bignum_set_bit(ret, n, 1);
78 * Input is in the first len words of a and b.
79 * Result is returned in the first 2*len words of c.
81 static void internal_mul(unsigned short *a, unsigned short *b,
82 unsigned short *c, int len)
87 for (j = 0; j < 2*len; j++)
90 for (i = len - 1; i >= 0; i--) {
93 for (j = len - 1; j >= 0; j--) {
94 t += ai * (unsigned long) b[j];
95 t += (unsigned long) c[i+j+1];
96 c[i+j+1] = (unsigned short)t;
99 c[i] = (unsigned short)t;
103 static void internal_add_shifted(unsigned short *number,
104 unsigned n, int shift) {
105 int word = 1 + (shift / 16);
106 int bshift = shift % 16;
107 unsigned long addend;
109 addend = n << bshift;
112 addend += number[word];
113 number[word] = (unsigned short) addend & 0xFFFF;
121 * Input in first alen words of a and first mlen words of m.
122 * Output in first alen words of a
123 * (of which first alen-mlen words will be zero).
124 * The MSW of m MUST have its high bit set.
125 * Quotient is accumulated in the `quotient' array, which is a Bignum
126 * rather than the internal bigendian format. Quotient parts are shifted
127 * left by `qshift' before adding into quot.
129 static void internal_mod(unsigned short *a, int alen,
130 unsigned short *m, int mlen,
131 unsigned short *quot, int qshift)
133 unsigned short m0, m1;
143 for (i = 0; i <= alen-mlen; i++) {
145 unsigned int q, r, c, ai1;
159 /* Find q = h:a[i] / m0 */
160 t = ((unsigned long) h << 16) + a[i];
164 /* Refine our estimate of q by looking at
165 h:a[i]:a[i+1] / m0:m1 */
166 t = (long) m1 * (long) q;
167 if (t > ((unsigned long) r << 16) + ai1) {
170 r = (r + m0) & 0xffff; /* overflow? */
171 if (r >= (unsigned long)m0 &&
172 t > ((unsigned long) r << 16) + ai1)
176 /* Subtract q * m from a[i...] */
178 for (k = mlen - 1; k >= 0; k--) {
179 t = (long) q * (long) m[k];
182 if ((unsigned short) t > a[i+k]) c++;
183 a[i+k] -= (unsigned short) t;
186 /* Add back m in case of borrow */
189 for (k = mlen - 1; k >= 0; k--) {
192 a[i+k] = (unsigned short)t;
198 internal_add_shifted(quot, q, qshift + 16 * (alen-mlen-i));
203 * Compute (base ^ exp) % mod.
204 * The base MUST be smaller than the modulus.
205 * The most significant word of mod MUST be non-zero.
206 * We assume that the result array is the same size as the mod array.
208 Bignum modpow(Bignum base, Bignum exp, Bignum mod)
210 unsigned short *a, *b, *n, *m;
215 /* Allocate m of size mlen, copy mod to m */
216 /* We use big endian internally */
218 m = smalloc(mlen * sizeof(unsigned short));
219 for (j = 0; j < mlen; j++) m[j] = mod[mod[0] - j];
221 /* Shift m left to make msb bit set */
222 for (mshift = 0; mshift < 15; mshift++)
223 if ((m[0] << mshift) & 0x8000) break;
225 for (i = 0; i < mlen - 1; i++)
226 m[i] = (m[i] << mshift) | (m[i+1] >> (16-mshift));
227 m[mlen-1] = m[mlen-1] << mshift;
230 /* Allocate n of size mlen, copy base to n */
231 n = smalloc(mlen * sizeof(unsigned short));
233 for (j = 0; j < i; j++) n[j] = 0;
234 for (j = 0; j < base[0]; j++) n[i+j] = base[base[0] - j];
236 /* Allocate a and b of size 2*mlen. Set a = 1 */
237 a = smalloc(2 * mlen * sizeof(unsigned short));
238 b = smalloc(2 * mlen * sizeof(unsigned short));
239 for (i = 0; i < 2*mlen; i++) a[i] = 0;
242 /* Skip leading zero bits of exp. */
244 while (i < exp[0] && (exp[exp[0] - i] & (1 << j)) == 0) {
246 if (j < 0) { i++; j = 15; }
249 /* Main computation */
252 internal_mul(a + mlen, a + mlen, b, mlen);
253 internal_mod(b, mlen*2, m, mlen, NULL, 0);
254 if ((exp[exp[0] - i] & (1 << j)) != 0) {
255 internal_mul(b + mlen, n, a, mlen);
256 internal_mod(a, mlen*2, m, mlen, NULL, 0);
266 /* Fixup result in case the modulus was shifted */
268 for (i = mlen - 1; i < 2*mlen - 1; i++)
269 a[i] = (a[i] << mshift) | (a[i+1] >> (16-mshift));
270 a[2*mlen-1] = a[2*mlen-1] << mshift;
271 internal_mod(a, mlen*2, m, mlen, NULL, 0);
272 for (i = 2*mlen - 1; i >= mlen; i--)
273 a[i] = (a[i] >> mshift) | (a[i-1] << (16-mshift));
276 /* Copy result to buffer */
277 result = newbn(mod[0]);
278 for (i = 0; i < mlen; i++)
279 result[result[0] - i] = a[i+mlen];
280 while (result[0] > 1 && result[result[0]] == 0) result[0]--;
282 /* Free temporary arrays */
283 for (i = 0; i < 2*mlen; i++) a[i] = 0; sfree(a);
284 for (i = 0; i < 2*mlen; i++) b[i] = 0; sfree(b);
285 for (i = 0; i < mlen; i++) m[i] = 0; sfree(m);
286 for (i = 0; i < mlen; i++) n[i] = 0; sfree(n);
292 * Compute (p * q) % mod.
293 * The most significant word of mod MUST be non-zero.
294 * We assume that the result array is the same size as the mod array.
296 Bignum modmul(Bignum p, Bignum q, Bignum mod)
298 unsigned short *a, *n, *m, *o;
300 int pqlen, mlen, rlen, i, j;
303 /* Allocate m of size mlen, copy mod to m */
304 /* We use big endian internally */
306 m = smalloc(mlen * sizeof(unsigned short));
307 for (j = 0; j < mlen; j++) m[j] = mod[mod[0] - j];
309 /* Shift m left to make msb bit set */
310 for (mshift = 0; mshift < 15; mshift++)
311 if ((m[0] << mshift) & 0x8000) break;
313 for (i = 0; i < mlen - 1; i++)
314 m[i] = (m[i] << mshift) | (m[i+1] >> (16-mshift));
315 m[mlen-1] = m[mlen-1] << mshift;
318 pqlen = (p[0] > q[0] ? p[0] : q[0]);
320 /* Allocate n of size pqlen, copy p to n */
321 n = smalloc(pqlen * sizeof(unsigned short));
323 for (j = 0; j < i; j++) n[j] = 0;
324 for (j = 0; j < p[0]; j++) n[i+j] = p[p[0] - j];
326 /* Allocate o of size pqlen, copy q to o */
327 o = smalloc(pqlen * sizeof(unsigned short));
329 for (j = 0; j < i; j++) o[j] = 0;
330 for (j = 0; j < q[0]; j++) o[i+j] = q[q[0] - j];
332 /* Allocate a of size 2*pqlen for result */
333 a = smalloc(2 * pqlen * sizeof(unsigned short));
335 /* Main computation */
336 internal_mul(n, o, a, pqlen);
337 internal_mod(a, pqlen*2, m, mlen, NULL, 0);
339 /* Fixup result in case the modulus was shifted */
341 for (i = 2*pqlen - mlen - 1; i < 2*pqlen - 1; i++)
342 a[i] = (a[i] << mshift) | (a[i+1] >> (16-mshift));
343 a[2*pqlen-1] = a[2*pqlen-1] << mshift;
344 internal_mod(a, pqlen*2, m, mlen, NULL, 0);
345 for (i = 2*pqlen - 1; i >= 2*pqlen - mlen; i--)
346 a[i] = (a[i] >> mshift) | (a[i-1] << (16-mshift));
349 /* Copy result to buffer */
350 rlen = (mlen < pqlen*2 ? mlen : pqlen*2);
351 result = newbn(rlen);
352 for (i = 0; i < rlen; i++)
353 result[result[0] - i] = a[i+2*pqlen-rlen];
354 while (result[0] > 1 && result[result[0]] == 0) result[0]--;
356 /* Free temporary arrays */
357 for (i = 0; i < 2*pqlen; i++) a[i] = 0; sfree(a);
358 for (i = 0; i < mlen; i++) m[i] = 0; sfree(m);
359 for (i = 0; i < pqlen; i++) n[i] = 0; sfree(n);
360 for (i = 0; i < pqlen; i++) o[i] = 0; sfree(o);
367 * The most significant word of mod MUST be non-zero.
368 * We assume that the result array is the same size as the mod array.
369 * We optionally write out a quotient.
371 void bigmod(Bignum p, Bignum mod, Bignum result, Bignum quotient)
373 unsigned short *n, *m;
375 int plen, mlen, i, j;
377 /* Allocate m of size mlen, copy mod to m */
378 /* We use big endian internally */
380 m = smalloc(mlen * sizeof(unsigned short));
381 for (j = 0; j < mlen; j++) m[j] = mod[mod[0] - j];
383 /* Shift m left to make msb bit set */
384 for (mshift = 0; mshift < 15; mshift++)
385 if ((m[0] << mshift) & 0x8000) break;
387 for (i = 0; i < mlen - 1; i++)
388 m[i] = (m[i] << mshift) | (m[i+1] >> (16-mshift));
389 m[mlen-1] = m[mlen-1] << mshift;
393 /* Ensure plen > mlen */
394 if (plen <= mlen) plen = mlen+1;
396 /* Allocate n of size plen, copy p to n */
397 n = smalloc(plen * sizeof(unsigned short));
398 for (j = 0; j < plen; j++) n[j] = 0;
399 for (j = 1; j <= p[0]; j++) n[plen-j] = p[j];
401 /* Main computation */
402 internal_mod(n, plen, m, mlen, quotient, mshift);
404 /* Fixup result in case the modulus was shifted */
406 for (i = plen - mlen - 1; i < plen - 1; i++)
407 n[i] = (n[i] << mshift) | (n[i+1] >> (16-mshift));
408 n[plen-1] = n[plen-1] << mshift;
409 internal_mod(n, plen, m, mlen, quotient, 0);
410 for (i = plen - 1; i >= plen - mlen; i--)
411 n[i] = (n[i] >> mshift) | (n[i-1] << (16-mshift));
414 /* Copy result to buffer */
415 for (i = 1; i <= result[0]; i++) {
417 result[i] = j>=0 ? n[j] : 0;
420 /* Free temporary arrays */
421 for (i = 0; i < mlen; i++) m[i] = 0; sfree(m);
422 for (i = 0; i < plen; i++) n[i] = 0; sfree(n);
426 * Decrement a number.
428 void decbn(Bignum bn) {
430 while (i < bn[0] && bn[i] == 0)
435 Bignum bignum_from_bytes(unsigned char *data, int nbytes) {
439 w = (nbytes+1)/2; /* bytes -> words */
444 for (i=nbytes; i-- ;) {
445 unsigned char byte = *data++;
447 result[1+i/2] |= byte<<8;
449 result[1+i/2] |= byte;
452 while (result[0] > 1 && result[result[0]] == 0) result[0]--;
457 * Read an ssh1-format bignum from a data buffer. Return the number
460 int ssh1_read_bignum(unsigned char *data, Bignum *result) {
461 unsigned char *p = data;
468 b = (w+7)/8; /* bits -> bytes */
470 if (!result) /* just return length */
473 *result = bignum_from_bytes(p, b);
479 * Return the bit count of a bignum, for ssh1 encoding.
481 int bignum_bitcount(Bignum bn) {
482 int bitcount = bn[0] * 16 - 1;
483 while (bitcount >= 0 && (bn[bitcount/16+1] >> (bitcount % 16)) == 0)
489 * Return the byte length of a bignum when ssh1 encoded.
491 int ssh1_bignum_length(Bignum bn) {
492 return 2 + (bignum_bitcount(bn)+7)/8;
496 * Return the byte length of a bignum when ssh2 encoded.
498 int ssh2_bignum_length(Bignum bn) {
499 return 4 + (bignum_bitcount(bn)+8)/8;
503 * Return a byte from a bignum; 0 is least significant, etc.
505 int bignum_byte(Bignum bn, int i) {
507 return 0; /* beyond the end */
509 return (bn[i/2+1] >> 8) & 0xFF;
511 return (bn[i/2+1] ) & 0xFF;
515 * Return a bit from a bignum; 0 is least significant, etc.
517 int bignum_bit(Bignum bn, int i) {
519 return 0; /* beyond the end */
521 return (bn[i/16+1] >> (i%16)) & 1;
525 * Set a bit in a bignum; 0 is least significant, etc.
527 void bignum_set_bit(Bignum bn, int bitnum, int value) {
528 if (bitnum >= 16*bn[0])
529 abort(); /* beyond the end */
532 int mask = 1 << (bitnum%16);
541 * Write a ssh1-format bignum into a buffer. It is assumed the
542 * buffer is big enough. Returns the number of bytes used.
544 int ssh1_write_bignum(void *data, Bignum bn) {
545 unsigned char *p = data;
546 int len = ssh1_bignum_length(bn);
548 int bitc = bignum_bitcount(bn);
550 *p++ = (bitc >> 8) & 0xFF;
551 *p++ = (bitc ) & 0xFF;
552 for (i = len-2; i-- ;)
553 *p++ = bignum_byte(bn, i);
558 * Compare two bignums. Returns like strcmp.
560 int bignum_cmp(Bignum a, Bignum b) {
561 int amax = a[0], bmax = b[0];
562 int i = (amax > bmax ? amax : bmax);
564 unsigned short aval = (i > amax ? 0 : a[i]);
565 unsigned short bval = (i > bmax ? 0 : b[i]);
566 if (aval < bval) return -1;
567 if (aval > bval) return +1;
574 * Right-shift one bignum to form another.
576 Bignum bignum_rshift(Bignum a, int shift) {
578 int i, shiftw, shiftb, shiftbb, bits;
579 unsigned short ai, ai1;
581 bits = bignum_bitcount(a) - shift;
582 ret = newbn((bits+15)/16);
587 shiftbb = 16 - shiftb;
590 for (i = 1; i <= ret[0]; i++) {
592 ai1 = (i+shiftw+1 <= a[0] ? a[i+shiftw+1] : 0);
593 ret[i] = ((ai >> shiftb) | (ai1 << shiftbb)) & 0xFFFF;
601 * Non-modular multiplication and addition.
603 Bignum bigmuladd(Bignum a, Bignum b, Bignum addend) {
604 int alen = a[0], blen = b[0];
605 int mlen = (alen > blen ? alen : blen);
606 int rlen, i, maxspot;
607 unsigned short *workspace;
610 /* mlen space for a, mlen space for b, 2*mlen for result */
611 workspace = smalloc(mlen * 4 * sizeof(unsigned short));
612 for (i = 0; i < mlen; i++) {
613 workspace[0*mlen + i] = (mlen-i <= a[0] ? a[mlen-i] : 0);
614 workspace[1*mlen + i] = (mlen-i <= b[0] ? b[mlen-i] : 0);
617 internal_mul(workspace+0*mlen, workspace+1*mlen, workspace+2*mlen, mlen);
619 /* now just copy the result back */
620 rlen = alen + blen + 1;
621 if (addend && rlen <= addend[0])
622 rlen = addend[0] + 1;
625 for (i = 1; i <= ret[0]; i++) {
626 ret[i] = (i <= 2*mlen ? workspace[4*mlen - i] : 0);
632 /* now add in the addend, if any */
634 unsigned long carry = 0;
635 for (i = 1; i <= rlen; i++) {
636 carry += (i <= ret[0] ? ret[i] : 0);
637 carry += (i <= addend[0] ? addend[i] : 0);
638 ret[i] = (unsigned short) carry & 0xFFFF;
640 if (ret[i] != 0 && i > maxspot)
650 * Non-modular multiplication.
652 Bignum bigmul(Bignum a, Bignum b) {
653 return bigmuladd(a, b, NULL);
657 * Create a bignum which is the bitmask covering another one. That
658 * is, the smallest integer which is >= N and is also one less than
661 Bignum bignum_bitmask(Bignum n) {
662 Bignum ret = copybn(n);
667 while (n[i] == 0 && i > 0)
670 return ret; /* input was zero */
681 * Convert a (max 16-bit) short into a bignum.
683 Bignum bignum_from_short(unsigned short n) {
688 ret[2] = (n >> 16) & 0xFFFF;
689 ret[0] = (ret[2] ? 2 : 1);
694 * Add a long to a bignum.
696 Bignum bignum_add_long(Bignum number, unsigned long addend) {
697 Bignum ret = newbn(number[0]+1);
699 unsigned long carry = 0;
701 for (i = 1; i <= ret[0]; i++) {
702 carry += addend & 0xFFFF;
703 carry += (i <= number[0] ? number[i] : 0);
705 ret[i] = (unsigned short) carry & 0xFFFF;
715 * Compute the residue of a bignum, modulo a (max 16-bit) short.
717 unsigned short bignum_mod_short(Bignum number, unsigned short modulus) {
718 unsigned long mod, r;
723 for (i = number[0]; i > 0; i--)
724 r = (r * 65536 + number[i]) % mod;
725 return (unsigned short) r;
728 void diagbn(char *prefix, Bignum md) {
729 int i, nibbles, morenibbles;
730 static const char hex[] = "0123456789ABCDEF";
732 debugprint(("%s0x", prefix ? prefix : ""));
734 nibbles = (3 + bignum_bitcount(md))/4; if (nibbles<1) nibbles=1;
735 morenibbles = 4*md[0] - nibbles;
736 for (i=0; i<morenibbles; i++) debugprint(("-"));
737 for (i=nibbles; i-- ;)
738 debugprint(("%c",hex[(bignum_byte(md, i/2) >> (4*(i%2))) & 0xF]));
740 if (prefix) debugprint(("\n"));
744 * Greatest common divisor.
746 Bignum biggcd(Bignum av, Bignum bv) {
747 Bignum a = copybn(av);
748 Bignum b = copybn(bv);
752 while (bignum_cmp(b, Zero) != 0) {
753 Bignum t = newbn(b[0]);
754 bigmod(a, b, t, NULL);
756 while (t[0] > 1 && t[t[0]] == 0) t[0]--;
767 * Modular inverse, using Euclid's extended algorithm.
769 Bignum modinv(Bignum number, Bignum modulus) {
770 Bignum a = copybn(modulus);
771 Bignum b = copybn(number);
772 Bignum xp = copybn(Zero);
773 Bignum x = copybn(One);
776 while (bignum_cmp(b, One) != 0) {
777 Bignum t = newbn(b[0]);
778 Bignum q = newbn(a[0]);
780 while (t[0] > 1 && t[t[0]] == 0) t[0]--;
786 x = bigmuladd(q, xp, t);
795 /* now we know that sign * x == 1, and that x < modulus */
797 /* set a new x to be modulus - x */
798 Bignum newx = newbn(modulus[0]);
799 unsigned short carry = 0;
803 for (i = 1; i <= newx[0]; i++) {
804 unsigned short aword = (i <= modulus[0] ? modulus[i] : 0);
805 unsigned short bword = (i <= x[0] ? x[i] : 0);
806 newx[i] = aword - bword - carry;
808 carry = carry ? (newx[i] >= bword) : (newx[i] > bword);
822 * Render a bignum into decimal. Return a malloced string holding
823 * the decimal representation.
825 char *bignum_decimal(Bignum x) {
830 unsigned short *workspace;
833 * First, estimate the number of digits. Since log(10)/log(2)
834 * is just greater than 93/28 (the joys of continued fraction
835 * approximations...) we know that for every 93 bits, we need
836 * at most 28 digits. This will tell us how much to malloc.
838 * Formally: if x has i bits, that means x is strictly less
839 * than 2^i. Since 2 is less than 10^(28/93), this is less than
840 * 10^(28i/93). We need an integer power of ten, so we must
841 * round up (rounding down might make it less than x again).
842 * Therefore if we multiply the bit count by 28/93, rounding
843 * up, we will have enough digits.
845 i = bignum_bitcount(x);
846 ndigits = (28*i + 92)/93; /* multiply by 28/93 and round up */
847 ndigits++; /* allow for trailing \0 */
848 ret = smalloc(ndigits);
851 * Now allocate some workspace to hold the binary form as we
852 * repeatedly divide it by ten. Initialise this to the
853 * big-endian form of the number.
855 workspace = smalloc(sizeof(unsigned short) * x[0]);
856 for (i = 0; i < x[0]; i++)
857 workspace[i] = x[x[0] - i];
860 * Next, write the decimal number starting with the last digit.
861 * We use ordinary short division, dividing 10 into the
869 for (i = 0; i < x[0]; i++) {
870 carry = (carry << 16) + workspace[i];
871 workspace[i] = (unsigned short) (carry / 10);
876 ret[--ndigit] = (char)(carry + '0');
880 * There's a chance we've fallen short of the start of the
881 * string. Correct if so.
884 memmove(ret, ret+ndigit, ndigits-ndigit);