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 ssh1_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 + (ssh1_bignum_bitcount(bn)+7)/8;
496 * Return a byte from a bignum; 0 is least significant, etc.
498 int bignum_byte(Bignum bn, int i) {
500 return 0; /* beyond the end */
502 return (bn[i/2+1] >> 8) & 0xFF;
504 return (bn[i/2+1] ) & 0xFF;
508 * Return a bit from a bignum; 0 is least significant, etc.
510 int bignum_bit(Bignum bn, int i) {
512 return 0; /* beyond the end */
514 return (bn[i/16+1] >> (i%16)) & 1;
518 * Set a bit in a bignum; 0 is least significant, etc.
520 void bignum_set_bit(Bignum bn, int bitnum, int value) {
521 if (bitnum >= 16*bn[0])
522 abort(); /* beyond the end */
525 int mask = 1 << (bitnum%16);
534 * Write a ssh1-format bignum into a buffer. It is assumed the
535 * buffer is big enough. Returns the number of bytes used.
537 int ssh1_write_bignum(void *data, Bignum bn) {
538 unsigned char *p = data;
539 int len = ssh1_bignum_length(bn);
541 int bitc = ssh1_bignum_bitcount(bn);
543 *p++ = (bitc >> 8) & 0xFF;
544 *p++ = (bitc ) & 0xFF;
545 for (i = len-2; i-- ;)
546 *p++ = bignum_byte(bn, i);
551 * Compare two bignums. Returns like strcmp.
553 int bignum_cmp(Bignum a, Bignum b) {
554 int amax = a[0], bmax = b[0];
555 int i = (amax > bmax ? amax : bmax);
557 unsigned short aval = (i > amax ? 0 : a[i]);
558 unsigned short bval = (i > bmax ? 0 : b[i]);
559 if (aval < bval) return -1;
560 if (aval > bval) return +1;
567 * Right-shift one bignum to form another.
569 Bignum bignum_rshift(Bignum a, int shift) {
571 int i, shiftw, shiftb, shiftbb, bits;
572 unsigned short ai, ai1;
574 bits = ssh1_bignum_bitcount(a) - shift;
575 ret = newbn((bits+15)/16);
580 shiftbb = 16 - shiftb;
583 for (i = 1; i <= ret[0]; i++) {
585 ai1 = (i+shiftw+1 <= a[0] ? a[i+shiftw+1] : 0);
586 ret[i] = ((ai >> shiftb) | (ai1 << shiftbb)) & 0xFFFF;
594 * Non-modular multiplication and addition.
596 Bignum bigmuladd(Bignum a, Bignum b, Bignum addend) {
597 int alen = a[0], blen = b[0];
598 int mlen = (alen > blen ? alen : blen);
599 int rlen, i, maxspot;
600 unsigned short *workspace;
603 /* mlen space for a, mlen space for b, 2*mlen for result */
604 workspace = smalloc(mlen * 4 * sizeof(unsigned short));
605 for (i = 0; i < mlen; i++) {
606 workspace[0*mlen + i] = (mlen-i <= a[0] ? a[mlen-i] : 0);
607 workspace[1*mlen + i] = (mlen-i <= b[0] ? b[mlen-i] : 0);
610 internal_mul(workspace+0*mlen, workspace+1*mlen, workspace+2*mlen, mlen);
612 /* now just copy the result back */
613 rlen = alen + blen + 1;
614 if (addend && rlen <= addend[0])
615 rlen = addend[0] + 1;
618 for (i = 1; i <= ret[0]; i++) {
619 ret[i] = (i <= 2*mlen ? workspace[4*mlen - i] : 0);
625 /* now add in the addend, if any */
627 unsigned long carry = 0;
628 for (i = 1; i <= rlen; i++) {
629 carry += (i <= ret[0] ? ret[i] : 0);
630 carry += (i <= addend[0] ? addend[i] : 0);
631 ret[i] = (unsigned short) carry & 0xFFFF;
633 if (ret[i] != 0 && i > maxspot)
643 * Non-modular multiplication.
645 Bignum bigmul(Bignum a, Bignum b) {
646 return bigmuladd(a, b, NULL);
650 * Create a bignum which is the bitmask covering another one. That
651 * is, the smallest integer which is >= N and is also one less than
654 Bignum bignum_bitmask(Bignum n) {
655 Bignum ret = copybn(n);
660 while (n[i] == 0 && i > 0)
663 return ret; /* input was zero */
674 * Convert a (max 16-bit) short into a bignum.
676 Bignum bignum_from_short(unsigned short n) {
681 ret[2] = (n >> 16) & 0xFFFF;
682 ret[0] = (ret[2] ? 2 : 1);
687 * Add a long to a bignum.
689 Bignum bignum_add_long(Bignum number, unsigned long addend) {
690 Bignum ret = newbn(number[0]+1);
692 unsigned long carry = 0;
694 for (i = 1; i <= ret[0]; i++) {
695 carry += addend & 0xFFFF;
696 carry += (i <= number[0] ? number[i] : 0);
698 ret[i] = (unsigned short) carry & 0xFFFF;
708 * Compute the residue of a bignum, modulo a (max 16-bit) short.
710 unsigned short bignum_mod_short(Bignum number, unsigned short modulus) {
711 unsigned long mod, r;
716 for (i = number[0]; i > 0; i--)
717 r = (r * 65536 + number[i]) % mod;
718 return (unsigned short) r;
721 void diagbn(char *prefix, Bignum md) {
722 int i, nibbles, morenibbles;
723 static const char hex[] = "0123456789ABCDEF";
725 debugprint(("%s0x", prefix ? prefix : ""));
727 nibbles = (3 + ssh1_bignum_bitcount(md))/4; if (nibbles<1) nibbles=1;
728 morenibbles = 4*md[0] - nibbles;
729 for (i=0; i<morenibbles; i++) debugprint(("-"));
730 for (i=nibbles; i-- ;)
731 debugprint(("%c",hex[(bignum_byte(md, i/2) >> (4*(i%2))) & 0xF]));
733 if (prefix) debugprint(("\n"));
737 * Greatest common divisor.
739 Bignum biggcd(Bignum av, Bignum bv) {
740 Bignum a = copybn(av);
741 Bignum b = copybn(bv);
745 while (bignum_cmp(b, Zero) != 0) {
746 Bignum t = newbn(b[0]);
747 bigmod(a, b, t, NULL);
749 while (t[0] > 1 && t[t[0]] == 0) t[0]--;
760 * Modular inverse, using Euclid's extended algorithm.
762 Bignum modinv(Bignum number, Bignum modulus) {
763 Bignum a = copybn(modulus);
764 Bignum b = copybn(number);
765 Bignum xp = copybn(Zero);
766 Bignum x = copybn(One);
769 while (bignum_cmp(b, One) != 0) {
770 Bignum t = newbn(b[0]);
771 Bignum q = newbn(a[0]);
773 while (t[0] > 1 && t[t[0]] == 0) t[0]--;
779 x = bigmuladd(q, xp, t);
788 /* now we know that sign * x == 1, and that x < modulus */
790 /* set a new x to be modulus - x */
791 Bignum newx = newbn(modulus[0]);
792 unsigned short carry = 0;
796 for (i = 1; i <= newx[0]; i++) {
797 unsigned short aword = (i <= modulus[0] ? modulus[i] : 0);
798 unsigned short bword = (i <= x[0] ? x[i] : 0);
799 newx[i] = aword - bword - carry;
801 carry = carry ? (newx[i] >= bword) : (newx[i] > bword);
815 * Render a bignum into decimal. Return a malloced string holding
816 * the decimal representation.
818 char *bignum_decimal(Bignum x) {
823 unsigned short *workspace;
826 * First, estimate the number of digits. Since log(10)/log(2)
827 * is just greater than 93/28 (the joys of continued fraction
828 * approximations...) we know that for every 93 bits, we need
829 * at most 28 digits. This will tell us how much to malloc.
831 * Formally: if x has i bits, that means x is strictly less
832 * than 2^i. Since 2 is less than 10^(28/93), this is less than
833 * 10^(28i/93). We need an integer power of ten, so we must
834 * round up (rounding down might make it less than x again).
835 * Therefore if we multiply the bit count by 28/93, rounding
836 * up, we will have enough digits.
838 i = ssh1_bignum_bitcount(x);
839 ndigits = (28*i + 92)/93; /* multiply by 28/93 and round up */
840 ndigits++; /* allow for trailing \0 */
841 ret = smalloc(ndigits);
844 * Now allocate some workspace to hold the binary form as we
845 * repeatedly divide it by ten. Initialise this to the
846 * big-endian form of the number.
848 workspace = smalloc(sizeof(unsigned short) * x[0]);
849 for (i = 0; i < x[0]; i++)
850 workspace[i] = x[x[0] - i];
853 * Next, write the decimal number starting with the last digit.
854 * We use ordinary short division, dividing 10 into the
862 for (i = 0; i < x[0]; i++) {
863 carry = (carry << 16) + workspace[i];
864 workspace[i] = (unsigned short) (carry / 10);
869 ret[--ndigit] = (char)(carry + '0');
873 * There's a chance we've fallen short of the start of the
874 * string. Correct if so.
877 memmove(ret, ret+ndigit, ndigits-ndigit);