source: filezilla/trunk/fuentes/src/putty/sshbn.c @ 130

Last change on this file since 130 was 130, checked in by jrpelegrina, 3 years ago

First release to xenial

File size: 54.0 KB
Line 
1/*
2 * Bignum routines for RSA and DH and stuff.
3 */
4
5#include <stdio.h>
6#include <assert.h>
7#include <stdlib.h>
8#include <string.h>
9#include <limits.h>
10#include <ctype.h>
11
12#include "misc.h"
13
14#include "sshbn.h"
15
16#define BIGNUM_INTERNAL
17typedef BignumInt *Bignum;
18
19#include "ssh.h"
20
21BignumInt bnZero[1] = { 0 };
22BignumInt bnOne[2] = { 1, 1 };
23BignumInt bnTen[2] = { 1, 10 };
24
25/*
26 * The Bignum format is an array of `BignumInt'. The first
27 * element of the array counts the remaining elements. The
28 * remaining elements express the actual number, base 2^BIGNUM_INT_BITS, _least_
29 * significant digit first. (So it's trivial to extract the bit
30 * with value 2^n for any n.)
31 *
32 * All Bignums in this module are positive. Negative numbers must
33 * be dealt with outside it.
34 *
35 * INVARIANT: the most significant word of any Bignum must be
36 * nonzero.
37 */
38
39Bignum Zero = bnZero, One = bnOne, Ten = bnTen;
40
41static Bignum newbn(int length)
42{
43    Bignum b;
44
45    assert(length >= 0 && length < INT_MAX / BIGNUM_INT_BITS);
46
47    b = snewn(length + 1, BignumInt);
48    memset(b, 0, (length + 1) * sizeof(*b));
49    b[0] = length;
50    return b;
51}
52
53void bn_restore_invariant(Bignum b)
54{
55    while (b[0] > 1 && b[b[0]] == 0)
56        b[0]--;
57}
58
59Bignum copybn(Bignum orig)
60{
61    Bignum b = snewn(orig[0] + 1, BignumInt);
62    if (!b)
63        abort();                       /* FIXME */
64    memcpy(b, orig, (orig[0] + 1) * sizeof(*b));
65    return b;
66}
67
68void freebn(Bignum b)
69{
70    /*
71     * Burn the evidence, just in case.
72     */
73    smemclr(b, sizeof(b[0]) * (b[0] + 1));
74    sfree(b);
75}
76
77Bignum bn_power_2(int n)
78{
79    Bignum ret;
80
81    assert(n >= 0);
82
83    ret = newbn(n / BIGNUM_INT_BITS + 1);
84    bignum_set_bit(ret, n, 1);
85    return ret;
86}
87
88/*
89 * Internal addition. Sets c = a - b, where 'a', 'b' and 'c' are all
90 * big-endian arrays of 'len' BignumInts. Returns a BignumInt carried
91 * off the top.
92 */
93static BignumInt internal_add(const BignumInt *a, const BignumInt *b,
94                              BignumInt *c, int len)
95{
96    int i;
97    BignumDblInt carry = 0;
98
99    for (i = len-1; i >= 0; i--) {
100        carry += (BignumDblInt)a[i] + b[i];
101        c[i] = (BignumInt)carry;
102        carry >>= BIGNUM_INT_BITS;
103    }
104
105    return (BignumInt)carry;
106}
107
108/*
109 * Internal subtraction. Sets c = a - b, where 'a', 'b' and 'c' are
110 * all big-endian arrays of 'len' BignumInts. Any borrow from the top
111 * is ignored.
112 */
113static void internal_sub(const BignumInt *a, const BignumInt *b,
114                         BignumInt *c, int len)
115{
116    int i;
117    BignumDblInt carry = 1;
118
119    for (i = len-1; i >= 0; i--) {
120        carry += (BignumDblInt)a[i] + (b[i] ^ BIGNUM_INT_MASK);
121        c[i] = (BignumInt)carry;
122        carry >>= BIGNUM_INT_BITS;
123    }
124}
125
126/*
127 * Compute c = a * b.
128 * Input is in the first len words of a and b.
129 * Result is returned in the first 2*len words of c.
130 *
131 * 'scratch' must point to an array of BignumInt of size at least
132 * mul_compute_scratch(len). (This covers the needs of internal_mul
133 * and all its recursive calls to itself.)
134 */
135#define KARATSUBA_THRESHOLD 50
136static int mul_compute_scratch(int len)
137{
138    int ret = 0;
139    while (len > KARATSUBA_THRESHOLD) {
140        int toplen = len/2, botlen = len - toplen; /* botlen is the bigger */
141        int midlen = botlen + 1;
142        ret += 4*midlen;
143        len = midlen;
144    }
145    return ret;
146}
147static void internal_mul(const BignumInt *a, const BignumInt *b,
148                         BignumInt *c, int len, BignumInt *scratch)
149{
150    if (len > KARATSUBA_THRESHOLD) {
151        int i;
152
153        /*
154         * Karatsuba divide-and-conquer algorithm. Cut each input in
155         * half, so that it's expressed as two big 'digits' in a giant
156         * base D:
157         *
158         *   a = a_1 D + a_0
159         *   b = b_1 D + b_0
160         *
161         * Then the product is of course
162         *
163         *  ab = a_1 b_1 D^2 + (a_1 b_0 + a_0 b_1) D + a_0 b_0
164         *
165         * and we compute the three coefficients by recursively
166         * calling ourself to do half-length multiplications.
167         *
168         * The clever bit that makes this worth doing is that we only
169         * need _one_ half-length multiplication for the central
170         * coefficient rather than the two that it obviouly looks
171         * like, because we can use a single multiplication to compute
172         *
173         *   (a_1 + a_0) (b_1 + b_0) = a_1 b_1 + a_1 b_0 + a_0 b_1 + a_0 b_0
174         *
175         * and then we subtract the other two coefficients (a_1 b_1
176         * and a_0 b_0) which we were computing anyway.
177         *
178         * Hence we get to multiply two numbers of length N in about
179         * three times as much work as it takes to multiply numbers of
180         * length N/2, which is obviously better than the four times
181         * as much work it would take if we just did a long
182         * conventional multiply.
183         */
184
185        int toplen = len/2, botlen = len - toplen; /* botlen is the bigger */
186        int midlen = botlen + 1;
187        BignumDblInt carry;
188#ifdef KARA_DEBUG
189        int i;
190#endif
191
192        /*
193         * The coefficients a_1 b_1 and a_0 b_0 just avoid overlapping
194         * in the output array, so we can compute them immediately in
195         * place.
196         */
197
198#ifdef KARA_DEBUG
199        printf("a1,a0 = 0x");
200        for (i = 0; i < len; i++) {
201            if (i == toplen) printf(", 0x");
202            printf("%0*x", BIGNUM_INT_BITS/4, a[i]);
203        }
204        printf("\n");
205        printf("b1,b0 = 0x");
206        for (i = 0; i < len; i++) {
207            if (i == toplen) printf(", 0x");
208            printf("%0*x", BIGNUM_INT_BITS/4, b[i]);
209        }
210        printf("\n");
211#endif
212
213        /* a_1 b_1 */
214        internal_mul(a, b, c, toplen, scratch);
215#ifdef KARA_DEBUG
216        printf("a1b1 = 0x");
217        for (i = 0; i < 2*toplen; i++) {
218            printf("%0*x", BIGNUM_INT_BITS/4, c[i]);
219        }
220        printf("\n");
221#endif
222
223        /* a_0 b_0 */
224        internal_mul(a + toplen, b + toplen, c + 2*toplen, botlen, scratch);
225#ifdef KARA_DEBUG
226        printf("a0b0 = 0x");
227        for (i = 0; i < 2*botlen; i++) {
228            printf("%0*x", BIGNUM_INT_BITS/4, c[2*toplen+i]);
229        }
230        printf("\n");
231#endif
232
233        /* Zero padding. midlen exceeds toplen by at most 2, so just
234         * zero the first two words of each input and the rest will be
235         * copied over. */
236        scratch[0] = scratch[1] = scratch[midlen] = scratch[midlen+1] = 0;
237
238        for (i = 0; i < toplen; i++) {
239            scratch[midlen - toplen + i] = a[i]; /* a_1 */
240            scratch[2*midlen - toplen + i] = b[i]; /* b_1 */
241        }
242
243        /* compute a_1 + a_0 */
244        scratch[0] = internal_add(scratch+1, a+toplen, scratch+1, botlen);
245#ifdef KARA_DEBUG
246        printf("a1plusa0 = 0x");
247        for (i = 0; i < midlen; i++) {
248            printf("%0*x", BIGNUM_INT_BITS/4, scratch[i]);
249        }
250        printf("\n");
251#endif
252        /* compute b_1 + b_0 */
253        scratch[midlen] = internal_add(scratch+midlen+1, b+toplen,
254                                       scratch+midlen+1, botlen);
255#ifdef KARA_DEBUG
256        printf("b1plusb0 = 0x");
257        for (i = 0; i < midlen; i++) {
258            printf("%0*x", BIGNUM_INT_BITS/4, scratch[midlen+i]);
259        }
260        printf("\n");
261#endif
262
263        /*
264         * Now we can do the third multiplication.
265         */
266        internal_mul(scratch, scratch + midlen, scratch + 2*midlen, midlen,
267                     scratch + 4*midlen);
268#ifdef KARA_DEBUG
269        printf("a1plusa0timesb1plusb0 = 0x");
270        for (i = 0; i < 2*midlen; i++) {
271            printf("%0*x", BIGNUM_INT_BITS/4, scratch[2*midlen+i]);
272        }
273        printf("\n");
274#endif
275
276        /*
277         * Now we can reuse the first half of 'scratch' to compute the
278         * sum of the outer two coefficients, to subtract from that
279         * product to obtain the middle one.
280         */
281        scratch[0] = scratch[1] = scratch[2] = scratch[3] = 0;
282        for (i = 0; i < 2*toplen; i++)
283            scratch[2*midlen - 2*toplen + i] = c[i];
284        scratch[1] = internal_add(scratch+2, c + 2*toplen,
285                                  scratch+2, 2*botlen);
286#ifdef KARA_DEBUG
287        printf("a1b1plusa0b0 = 0x");
288        for (i = 0; i < 2*midlen; i++) {
289            printf("%0*x", BIGNUM_INT_BITS/4, scratch[i]);
290        }
291        printf("\n");
292#endif
293
294        internal_sub(scratch + 2*midlen, scratch,
295                     scratch + 2*midlen, 2*midlen);
296#ifdef KARA_DEBUG
297        printf("a1b0plusa0b1 = 0x");
298        for (i = 0; i < 2*midlen; i++) {
299            printf("%0*x", BIGNUM_INT_BITS/4, scratch[2*midlen+i]);
300        }
301        printf("\n");
302#endif
303
304        /*
305         * And now all we need to do is to add that middle coefficient
306         * back into the output. We may have to propagate a carry
307         * further up the output, but we can be sure it won't
308         * propagate right the way off the top.
309         */
310        carry = internal_add(c + 2*len - botlen - 2*midlen,
311                             scratch + 2*midlen,
312                             c + 2*len - botlen - 2*midlen, 2*midlen);
313        i = 2*len - botlen - 2*midlen - 1;
314        while (carry) {
315            assert(i >= 0);
316            carry += c[i];
317            c[i] = (BignumInt)carry;
318            carry >>= BIGNUM_INT_BITS;
319            i--;
320        }
321#ifdef KARA_DEBUG
322        printf("ab = 0x");
323        for (i = 0; i < 2*len; i++) {
324            printf("%0*x", BIGNUM_INT_BITS/4, c[i]);
325        }
326        printf("\n");
327#endif
328
329    } else {
330        int i;
331        BignumInt carry;
332        BignumDblInt t;
333        const BignumInt *ap, *bp;
334        BignumInt *cp, *cps;
335
336        /*
337         * Multiply in the ordinary O(N^2) way.
338         */
339
340        for (i = 0; i < 2 * len; i++)
341            c[i] = 0;
342
343        for (cps = c + 2*len, ap = a + len; ap-- > a; cps--) {
344            carry = 0;
345            for (cp = cps, bp = b + len; cp--, bp-- > b ;) {
346                t = (MUL_WORD(*ap, *bp) + carry) + *cp;
347                *cp = (BignumInt) t;
348                carry = (BignumInt)(t >> BIGNUM_INT_BITS);
349            }
350            *cp = carry;
351        }
352    }
353}
354
355/*
356 * Variant form of internal_mul used for the initial step of
357 * Montgomery reduction. Only bothers outputting 'len' words
358 * (everything above that is thrown away).
359 */
360static void internal_mul_low(const BignumInt *a, const BignumInt *b,
361                             BignumInt *c, int len, BignumInt *scratch)
362{
363    if (len > KARATSUBA_THRESHOLD) {
364        int i;
365
366        /*
367         * Karatsuba-aware version of internal_mul_low. As before, we
368         * express each input value as a shifted combination of two
369         * halves:
370         *
371         *   a = a_1 D + a_0
372         *   b = b_1 D + b_0
373         *
374         * Then the full product is, as before,
375         *
376         *  ab = a_1 b_1 D^2 + (a_1 b_0 + a_0 b_1) D + a_0 b_0
377         *
378         * Provided we choose D on the large side (so that a_0 and b_0
379         * are _at least_ as long as a_1 and b_1), we don't need the
380         * topmost term at all, and we only need half of the middle
381         * term. So there's no point in doing the proper Karatsuba
382         * optimisation which computes the middle term using the top
383         * one, because we'd take as long computing the top one as
384         * just computing the middle one directly.
385         *
386         * So instead, we do a much more obvious thing: we call the
387         * fully optimised internal_mul to compute a_0 b_0, and we
388         * recursively call ourself to compute the _bottom halves_ of
389         * a_1 b_0 and a_0 b_1, each of which we add into the result
390         * in the obvious way.
391         *
392         * In other words, there's no actual Karatsuba _optimisation_
393         * in this function; the only benefit in doing it this way is
394         * that we call internal_mul proper for a large part of the
395         * work, and _that_ can optimise its operation.
396         */
397
398        int toplen = len/2, botlen = len - toplen; /* botlen is the bigger */
399
400        /*
401         * Scratch space for the various bits and pieces we're going
402         * to be adding together: we need botlen*2 words for a_0 b_0
403         * (though we may end up throwing away its topmost word), and
404         * toplen words for each of a_1 b_0 and a_0 b_1. That adds up
405         * to exactly 2*len.
406         */
407
408        /* a_0 b_0 */
409        internal_mul(a + toplen, b + toplen, scratch + 2*toplen, botlen,
410                     scratch + 2*len);
411
412        /* a_1 b_0 */
413        internal_mul_low(a, b + len - toplen, scratch + toplen, toplen,
414                         scratch + 2*len);
415
416        /* a_0 b_1 */
417        internal_mul_low(a + len - toplen, b, scratch, toplen,
418                         scratch + 2*len);
419
420        /* Copy the bottom half of the big coefficient into place */
421        for (i = 0; i < botlen; i++)
422            c[toplen + i] = scratch[2*toplen + botlen + i];
423
424        /* Add the two small coefficients, throwing away the returned carry */
425        internal_add(scratch, scratch + toplen, scratch, toplen);
426
427        /* And add that to the large coefficient, leaving the result in c. */
428        internal_add(scratch, scratch + 2*toplen + botlen - toplen,
429                     c, toplen);
430
431    } else {
432        int i;
433        BignumInt carry;
434        BignumDblInt t;
435        const BignumInt *ap, *bp;
436        BignumInt *cp, *cps;
437
438        /*
439         * Multiply in the ordinary O(N^2) way.
440         */
441
442        for (i = 0; i < len; i++)
443            c[i] = 0;
444
445        for (cps = c + len, ap = a + len; ap-- > a; cps--) {
446            carry = 0;
447            for (cp = cps, bp = b + len; bp--, cp-- > c ;) {
448                t = (MUL_WORD(*ap, *bp) + carry) + *cp;
449                *cp = (BignumInt) t;
450                carry = (BignumInt)(t >> BIGNUM_INT_BITS);
451            }
452        }
453    }
454}
455
456/*
457 * Montgomery reduction. Expects x to be a big-endian array of 2*len
458 * BignumInts whose value satisfies 0 <= x < rn (where r = 2^(len *
459 * BIGNUM_INT_BITS) is the Montgomery base). Returns in the same array
460 * a value x' which is congruent to xr^{-1} mod n, and satisfies 0 <=
461 * x' < n.
462 *
463 * 'n' and 'mninv' should be big-endian arrays of 'len' BignumInts
464 * each, containing respectively n and the multiplicative inverse of
465 * -n mod r.
466 *
467 * 'tmp' is an array of BignumInt used as scratch space, of length at
468 * least 3*len + mul_compute_scratch(len).
469 */
470static void monty_reduce(BignumInt *x, const BignumInt *n,
471                         const BignumInt *mninv, BignumInt *tmp, int len)
472{
473    int i;
474    BignumInt carry;
475
476    /*
477     * Multiply x by (-n)^{-1} mod r. This gives us a value m such
478     * that mn is congruent to -x mod r. Hence, mn+x is an exact
479     * multiple of r, and is also (obviously) congruent to x mod n.
480     */
481    internal_mul_low(x + len, mninv, tmp, len, tmp + 3*len);
482
483    /*
484     * Compute t = (mn+x)/r in ordinary, non-modular, integer
485     * arithmetic. By construction this is exact, and is congruent mod
486     * n to x * r^{-1}, i.e. the answer we want.
487     *
488     * The following multiply leaves that answer in the _most_
489     * significant half of the 'x' array, so then we must shift it
490     * down.
491     */
492    internal_mul(tmp, n, tmp+len, len, tmp + 3*len);
493    carry = internal_add(x, tmp+len, x, 2*len);
494    for (i = 0; i < len; i++)
495        x[len + i] = x[i], x[i] = 0;
496
497    /*
498     * Reduce t mod n. This doesn't require a full-on division by n,
499     * but merely a test and single optional subtraction, since we can
500     * show that 0 <= t < 2n.
501     *
502     * Proof:
503     *  + we computed m mod r, so 0 <= m < r.
504     *  + so 0 <= mn < rn, obviously
505     *  + hence we only need 0 <= x < rn to guarantee that 0 <= mn+x < 2rn
506     *  + yielding 0 <= (mn+x)/r < 2n as required.
507     */
508    if (!carry) {
509        for (i = 0; i < len; i++)
510            if (x[len + i] != n[i])
511                break;
512    }
513    if (carry || i >= len || x[len + i] > n[i])
514        internal_sub(x+len, n, x+len, len);
515}
516
517static void internal_add_shifted(BignumInt *number,
518                                 BignumInt n, int shift)
519{
520    int word = 1 + (shift / BIGNUM_INT_BITS);
521    int bshift = shift % BIGNUM_INT_BITS;
522    BignumDblInt addend;
523
524    addend = (BignumDblInt)n << bshift;
525
526    while (addend) {
527        assert(word <= number[0]);
528        addend += number[word];
529        number[word] = (BignumInt) addend & BIGNUM_INT_MASK;
530        addend >>= BIGNUM_INT_BITS;
531        word++;
532    }
533}
534
535/*
536 * Compute a = a % m.
537 * Input in first alen words of a and first mlen words of m.
538 * Output in first alen words of a
539 * (of which first alen-mlen words will be zero).
540 * The MSW of m MUST have its high bit set.
541 * Quotient is accumulated in the `quotient' array, which is a Bignum
542 * rather than the internal bigendian format. Quotient parts are shifted
543 * left by `qshift' before adding into quot.
544 */
545static void internal_mod(BignumInt *a, int alen,
546                         BignumInt *m, int mlen,
547                         BignumInt *quot, int qshift)
548{
549    BignumInt m0, m1, h;
550    int i, k;
551
552    m0 = m[0];
553    assert(m0 >> (BIGNUM_INT_BITS-1) == 1);
554    if (mlen > 1)
555        m1 = m[1];
556    else
557        m1 = 0;
558
559    for (i = 0; i <= alen - mlen; i++) {
560        BignumDblInt t;
561        BignumInt q, r, c, ai1;
562
563        if (i == 0) {
564            h = 0;
565        } else {
566            h = a[i - 1];
567            a[i - 1] = 0;
568        }
569
570        if (i == alen - 1)
571            ai1 = 0;
572        else
573            ai1 = a[i + 1];
574
575        /* Find q = h:a[i] / m0 */
576        if (h >= m0) {
577            /*
578             * Special case.
579             *
580             * To illustrate it, suppose a BignumInt is 8 bits, and
581             * we are dividing (say) A1:23:45:67 by A1:B2:C3. Then
582             * our initial division will be 0xA123 / 0xA1, which
583             * will give a quotient of 0x100 and a divide overflow.
584             * However, the invariants in this division algorithm
585             * are not violated, since the full number A1:23:... is
586             * _less_ than the quotient prefix A1:B2:... and so the
587             * following correction loop would have sorted it out.
588             *
589             * In this situation we set q to be the largest
590             * quotient we _can_ stomach (0xFF, of course).
591             */
592            q = BIGNUM_INT_MASK;
593        } else {
594            /* Macro doesn't want an array subscript expression passed
595             * into it (see definition), so use a temporary. */
596            BignumInt tmplo = a[i];
597            DIVMOD_WORD(q, r, h, tmplo, m0);
598
599            /* Refine our estimate of q by looking at
600             h:a[i]:a[i+1] / m0:m1 */
601            t = MUL_WORD(m1, q);
602            if (t > ((BignumDblInt) r << BIGNUM_INT_BITS) + ai1) {
603                q--;
604                t -= m1;
605                r = (r + m0) & BIGNUM_INT_MASK;     /* overflow? */
606                if (r >= (BignumDblInt) m0 &&
607                    t > ((BignumDblInt) r << BIGNUM_INT_BITS) + ai1) q--;
608            }
609        }
610
611        /* Subtract q * m from a[i...] */
612        c = 0;
613        for (k = mlen - 1; k >= 0; k--) {
614            t = MUL_WORD(q, m[k]);
615            t += c;
616            c = (BignumInt)(t >> BIGNUM_INT_BITS);
617            if ((BignumInt) t > a[i + k])
618                c++;
619            a[i + k] -= (BignumInt) t;
620        }
621
622        /* Add back m in case of borrow */
623        if (c != h) {
624            t = 0;
625            for (k = mlen - 1; k >= 0; k--) {
626                t += m[k];
627                t += a[i + k];
628                a[i + k] = (BignumInt) t;
629                t = t >> BIGNUM_INT_BITS;
630            }
631            q--;
632        }
633        if (quot)
634            internal_add_shifted(quot, q, qshift + BIGNUM_INT_BITS * (alen - mlen - i));
635    }
636}
637
638/*
639 * Compute (base ^ exp) % mod, the pedestrian way.
640 */
641Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod)
642{
643    BignumInt *a, *b, *n, *m, *scratch;
644    int mshift;
645    int mlen, scratchlen, i, j;
646    Bignum base, result;
647
648    /*
649     * The most significant word of mod needs to be non-zero. It
650     * should already be, but let's make sure.
651     */
652    assert(mod[mod[0]] != 0);
653
654    /*
655     * Make sure the base is smaller than the modulus, by reducing
656     * it modulo the modulus if not.
657     */
658    base = bigmod(base_in, mod);
659
660    /* Allocate m of size mlen, copy mod to m */
661    /* We use big endian internally */
662    mlen = mod[0];
663    m = snewn(mlen, BignumInt);
664    for (j = 0; j < mlen; j++)
665        m[j] = mod[mod[0] - j];
666
667    /* Shift m left to make msb bit set */
668    for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)
669        if ((m[0] << mshift) & BIGNUM_TOP_BIT)
670            break;
671    if (mshift) {
672        for (i = 0; i < mlen - 1; i++)
673            m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));
674        m[mlen - 1] = m[mlen - 1] << mshift;
675    }
676
677    /* Allocate n of size mlen, copy base to n */
678    n = snewn(mlen, BignumInt);
679    i = mlen - base[0];
680    for (j = 0; j < i; j++)
681        n[j] = 0;
682    for (j = 0; j < (int)base[0]; j++)
683        n[i + j] = base[base[0] - j];
684
685    /* Allocate a and b of size 2*mlen. Set a = 1 */
686    a = snewn(2 * mlen, BignumInt);
687    b = snewn(2 * mlen, BignumInt);
688    for (i = 0; i < 2 * mlen; i++)
689        a[i] = 0;
690    a[2 * mlen - 1] = 1;
691
692    /* Scratch space for multiplies */
693    scratchlen = mul_compute_scratch(mlen);
694    scratch = snewn(scratchlen, BignumInt);
695
696    /* Skip leading zero bits of exp. */
697    i = 0;
698    j = BIGNUM_INT_BITS-1;
699    while (i < (int)exp[0] && (exp[exp[0] - i] & ((BignumInt)1 << j)) == 0) {
700        j--;
701        if (j < 0) {
702            i++;
703            j = BIGNUM_INT_BITS-1;
704        }
705    }
706
707    /* Main computation */
708    while (i < (int)exp[0]) {
709        while (j >= 0) {
710            internal_mul(a + mlen, a + mlen, b, mlen, scratch);
711            internal_mod(b, mlen * 2, m, mlen, NULL, 0);
712            if ((exp[exp[0] - i] & ((BignumInt)1 << j)) != 0) {
713                internal_mul(b + mlen, n, a, mlen, scratch);
714                internal_mod(a, mlen * 2, m, mlen, NULL, 0);
715            } else {
716                BignumInt *t;
717                t = a;
718                a = b;
719                b = t;
720            }
721            j--;
722        }
723        i++;
724        j = BIGNUM_INT_BITS-1;
725    }
726
727    /* Fixup result in case the modulus was shifted */
728    if (mshift) {
729        for (i = mlen - 1; i < 2 * mlen - 1; i++)
730            a[i] = (a[i] << mshift) | (a[i + 1] >> (BIGNUM_INT_BITS - mshift));
731        a[2 * mlen - 1] = a[2 * mlen - 1] << mshift;
732        internal_mod(a, mlen * 2, m, mlen, NULL, 0);
733        for (i = 2 * mlen - 1; i >= mlen; i--)
734            a[i] = (a[i] >> mshift) | (a[i - 1] << (BIGNUM_INT_BITS - mshift));
735    }
736
737    /* Copy result to buffer */
738    result = newbn(mod[0]);
739    for (i = 0; i < mlen; i++)
740        result[result[0] - i] = a[i + mlen];
741    while (result[0] > 1 && result[result[0]] == 0)
742        result[0]--;
743
744    /* Free temporary arrays */
745    smemclr(a, 2 * mlen * sizeof(*a));
746    sfree(a);
747    smemclr(scratch, scratchlen * sizeof(*scratch));
748    sfree(scratch);
749    smemclr(b, 2 * mlen * sizeof(*b));
750    sfree(b);
751    smemclr(m, mlen * sizeof(*m));
752    sfree(m);
753    smemclr(n, mlen * sizeof(*n));
754    sfree(n);
755
756    freebn(base);
757
758    return result;
759}
760
761/*
762 * Compute (base ^ exp) % mod. Uses the Montgomery multiplication
763 * technique where possible, falling back to modpow_simple otherwise.
764 */
765Bignum modpow(Bignum base_in, Bignum exp, Bignum mod)
766{
767    BignumInt *a, *b, *x, *n, *mninv, *scratch;
768    int len, scratchlen, i, j;
769    Bignum base, base2, r, rn, inv, result;
770
771    /*
772     * The most significant word of mod needs to be non-zero. It
773     * should already be, but let's make sure.
774     */
775    assert(mod[mod[0]] != 0);
776
777    /*
778     * mod had better be odd, or we can't do Montgomery multiplication
779     * using a power of two at all.
780     */
781    if (!(mod[1] & 1))
782        return modpow_simple(base_in, exp, mod);
783
784    /*
785     * Make sure the base is smaller than the modulus, by reducing
786     * it modulo the modulus if not.
787     */
788    base = bigmod(base_in, mod);
789
790    /*
791     * Compute the inverse of n mod r, for monty_reduce. (In fact we
792     * want the inverse of _minus_ n mod r, but we'll sort that out
793     * below.)
794     */
795    len = mod[0];
796    r = bn_power_2(BIGNUM_INT_BITS * len);
797    inv = modinv(mod, r);
798    assert(inv); /* cannot fail, since mod is odd and r is a power of 2 */
799
800    /*
801     * Multiply the base by r mod n, to get it into Montgomery
802     * representation.
803     */
804    base2 = modmul(base, r, mod);
805    freebn(base);
806    base = base2;
807
808    rn = bigmod(r, mod);               /* r mod n, i.e. Montgomerified 1 */
809
810    freebn(r);                         /* won't need this any more */
811
812    /*
813     * Set up internal arrays of the right lengths, in big-endian
814     * format, containing the base, the modulus, and the modulus's
815     * inverse.
816     */
817    n = snewn(len, BignumInt);
818    for (j = 0; j < len; j++)
819        n[len - 1 - j] = mod[j + 1];
820
821    mninv = snewn(len, BignumInt);
822    for (j = 0; j < len; j++)
823        mninv[len - 1 - j] = (j < (int)inv[0] ? inv[j + 1] : 0);
824    freebn(inv);         /* we don't need this copy of it any more */
825    /* Now negate mninv mod r, so it's the inverse of -n rather than +n. */
826    x = snewn(len, BignumInt);
827    for (j = 0; j < len; j++)
828        x[j] = 0;
829    internal_sub(x, mninv, mninv, len);
830
831    /* x = snewn(len, BignumInt); */ /* already done above */
832    for (j = 0; j < len; j++)
833        x[len - 1 - j] = (j < (int)base[0] ? base[j + 1] : 0);
834    freebn(base);        /* we don't need this copy of it any more */
835
836    a = snewn(2*len, BignumInt);
837    b = snewn(2*len, BignumInt);
838    for (j = 0; j < len; j++)
839        a[2*len - 1 - j] = (j < (int)rn[0] ? rn[j + 1] : 0);
840    freebn(rn);
841
842    /* Scratch space for multiplies */
843    scratchlen = 3*len + mul_compute_scratch(len);
844    scratch = snewn(scratchlen, BignumInt);
845
846    /* Skip leading zero bits of exp. */
847    i = 0;
848    j = BIGNUM_INT_BITS-1;
849    while (i < (int)exp[0] && (exp[exp[0] - i] & ((BignumInt)1 << j)) == 0) {
850        j--;
851        if (j < 0) {
852            i++;
853            j = BIGNUM_INT_BITS-1;
854        }
855    }
856
857    /* Main computation */
858    while (i < (int)exp[0]) {
859        while (j >= 0) {
860            internal_mul(a + len, a + len, b, len, scratch);
861            monty_reduce(b, n, mninv, scratch, len);
862            if ((exp[exp[0] - i] & ((BignumInt)1 << j)) != 0) {
863                internal_mul(b + len, x, a, len,  scratch);
864                monty_reduce(a, n, mninv, scratch, len);
865            } else {
866                BignumInt *t;
867                t = a;
868                a = b;
869                b = t;
870            }
871            j--;
872        }
873        i++;
874        j = BIGNUM_INT_BITS-1;
875    }
876
877    /*
878     * Final monty_reduce to get back from the adjusted Montgomery
879     * representation.
880     */
881    monty_reduce(a, n, mninv, scratch, len);
882
883    /* Copy result to buffer */
884    result = newbn(mod[0]);
885    for (i = 0; i < len; i++)
886        result[result[0] - i] = a[i + len];
887    while (result[0] > 1 && result[result[0]] == 0)
888        result[0]--;
889
890    /* Free temporary arrays */
891    smemclr(scratch, scratchlen * sizeof(*scratch));
892    sfree(scratch);
893    smemclr(a, 2 * len * sizeof(*a));
894    sfree(a);
895    smemclr(b, 2 * len * sizeof(*b));
896    sfree(b);
897    smemclr(mninv, len * sizeof(*mninv));
898    sfree(mninv);
899    smemclr(n, len * sizeof(*n));
900    sfree(n);
901    smemclr(x, len * sizeof(*x));
902    sfree(x);
903
904    return result;
905}
906
907/*
908 * Compute (p * q) % mod.
909 * The most significant word of mod MUST be non-zero.
910 * We assume that the result array is the same size as the mod array.
911 */
912Bignum modmul(Bignum p, Bignum q, Bignum mod)
913{
914    BignumInt *a, *n, *m, *o, *scratch;
915    int mshift, scratchlen;
916    int pqlen, mlen, rlen, i, j;
917    Bignum result;
918
919    /*
920     * The most significant word of mod needs to be non-zero. It
921     * should already be, but let's make sure.
922     */
923    assert(mod[mod[0]] != 0);
924
925    /* Allocate m of size mlen, copy mod to m */
926    /* We use big endian internally */
927    mlen = mod[0];
928    m = snewn(mlen, BignumInt);
929    for (j = 0; j < mlen; j++)
930        m[j] = mod[mod[0] - j];
931
932    /* Shift m left to make msb bit set */
933    for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)
934        if ((m[0] << mshift) & BIGNUM_TOP_BIT)
935            break;
936    if (mshift) {
937        for (i = 0; i < mlen - 1; i++)
938            m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));
939        m[mlen - 1] = m[mlen - 1] << mshift;
940    }
941
942    pqlen = (p[0] > q[0] ? p[0] : q[0]);
943
944    /*
945     * Make sure that we're allowing enough space. The shifting below
946     * will underflow the vectors we allocate if pqlen is too small.
947     */
948    if (2*pqlen <= mlen)
949        pqlen = mlen/2 + 1;
950
951    /* Allocate n of size pqlen, copy p to n */
952    n = snewn(pqlen, BignumInt);
953    i = pqlen - p[0];
954    for (j = 0; j < i; j++)
955        n[j] = 0;
956    for (j = 0; j < (int)p[0]; j++)
957        n[i + j] = p[p[0] - j];
958
959    /* Allocate o of size pqlen, copy q to o */
960    o = snewn(pqlen, BignumInt);
961    i = pqlen - q[0];
962    for (j = 0; j < i; j++)
963        o[j] = 0;
964    for (j = 0; j < (int)q[0]; j++)
965        o[i + j] = q[q[0] - j];
966
967    /* Allocate a of size 2*pqlen for result */
968    a = snewn(2 * pqlen, BignumInt);
969
970    /* Scratch space for multiplies */
971    scratchlen = mul_compute_scratch(pqlen);
972    scratch = snewn(scratchlen, BignumInt);
973
974    /* Main computation */
975    internal_mul(n, o, a, pqlen, scratch);
976    internal_mod(a, pqlen * 2, m, mlen, NULL, 0);
977
978    /* Fixup result in case the modulus was shifted */
979    if (mshift) {
980        for (i = 2 * pqlen - mlen - 1; i < 2 * pqlen - 1; i++)
981            a[i] = (a[i] << mshift) | (a[i + 1] >> (BIGNUM_INT_BITS - mshift));
982        a[2 * pqlen - 1] = a[2 * pqlen - 1] << mshift;
983        internal_mod(a, pqlen * 2, m, mlen, NULL, 0);
984        for (i = 2 * pqlen - 1; i >= 2 * pqlen - mlen; i--)
985            a[i] = (a[i] >> mshift) | (a[i - 1] << (BIGNUM_INT_BITS - mshift));
986    }
987
988    /* Copy result to buffer */
989    rlen = (mlen < pqlen * 2 ? mlen : pqlen * 2);
990    result = newbn(rlen);
991    for (i = 0; i < rlen; i++)
992        result[result[0] - i] = a[i + 2 * pqlen - rlen];
993    while (result[0] > 1 && result[result[0]] == 0)
994        result[0]--;
995
996    /* Free temporary arrays */
997    smemclr(scratch, scratchlen * sizeof(*scratch));
998    sfree(scratch);
999    smemclr(a, 2 * pqlen * sizeof(*a));
1000    sfree(a);
1001    smemclr(m, mlen * sizeof(*m));
1002    sfree(m);
1003    smemclr(n, pqlen * sizeof(*n));
1004    sfree(n);
1005    smemclr(o, pqlen * sizeof(*o));
1006    sfree(o);
1007
1008    return result;
1009}
1010
1011Bignum modsub(const Bignum a, const Bignum b, const Bignum n)
1012{
1013    Bignum a1, b1, ret;
1014
1015    if (bignum_cmp(a, n) >= 0) a1 = bigmod(a, n);
1016    else a1 = a;
1017    if (bignum_cmp(b, n) >= 0) b1 = bigmod(b, n);
1018    else b1 = b;
1019
1020    if (bignum_cmp(a1, b1) >= 0) /* a >= b */
1021    {
1022        ret = bigsub(a1, b1);
1023    }
1024    else
1025    {
1026        /* Handle going round the corner of the modulus without having
1027         * negative support in Bignum */
1028        Bignum tmp = bigsub(n, b1);
1029        assert(tmp);
1030        ret = bigadd(tmp, a1);
1031        freebn(tmp);
1032    }
1033
1034    if (a != a1) freebn(a1);
1035    if (b != b1) freebn(b1);
1036
1037    return ret;
1038}
1039
1040/*
1041 * Compute p % mod.
1042 * The most significant word of mod MUST be non-zero.
1043 * We assume that the result array is the same size as the mod array.
1044 * We optionally write out a quotient if `quotient' is non-NULL.
1045 * We can avoid writing out the result if `result' is NULL.
1046 */
1047static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient)
1048{
1049    BignumInt *n, *m;
1050    int mshift;
1051    int plen, mlen, i, j;
1052
1053    /*
1054     * The most significant word of mod needs to be non-zero. It
1055     * should already be, but let's make sure.
1056     */
1057    assert(mod[mod[0]] != 0);
1058
1059    /* Allocate m of size mlen, copy mod to m */
1060    /* We use big endian internally */
1061    mlen = mod[0];
1062    m = snewn(mlen, BignumInt);
1063    for (j = 0; j < mlen; j++)
1064        m[j] = mod[mod[0] - j];
1065
1066    /* Shift m left to make msb bit set */
1067    for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)
1068        if ((m[0] << mshift) & BIGNUM_TOP_BIT)
1069            break;
1070    if (mshift) {
1071        for (i = 0; i < mlen - 1; i++)
1072            m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));
1073        m[mlen - 1] = m[mlen - 1] << mshift;
1074    }
1075
1076    plen = p[0];
1077    /* Ensure plen > mlen */
1078    if (plen <= mlen)
1079        plen = mlen + 1;
1080
1081    /* Allocate n of size plen, copy p to n */
1082    n = snewn(plen, BignumInt);
1083    for (j = 0; j < plen; j++)
1084        n[j] = 0;
1085    for (j = 1; j <= (int)p[0]; j++)
1086        n[plen - j] = p[j];
1087
1088    /* Main computation */
1089    internal_mod(n, plen, m, mlen, quotient, mshift);
1090
1091    /* Fixup result in case the modulus was shifted */
1092    if (mshift) {
1093        for (i = plen - mlen - 1; i < plen - 1; i++)
1094            n[i] = (n[i] << mshift) | (n[i + 1] >> (BIGNUM_INT_BITS - mshift));
1095        n[plen - 1] = n[plen - 1] << mshift;
1096        internal_mod(n, plen, m, mlen, quotient, 0);
1097        for (i = plen - 1; i >= plen - mlen; i--)
1098            n[i] = (n[i] >> mshift) | (n[i - 1] << (BIGNUM_INT_BITS - mshift));
1099    }
1100
1101    /* Copy result to buffer */
1102    if (result) {
1103        for (i = 1; i <= (int)result[0]; i++) {
1104            int j = plen - i;
1105            result[i] = j >= 0 ? n[j] : 0;
1106        }
1107    }
1108
1109    /* Free temporary arrays */
1110    smemclr(m, mlen * sizeof(*m));
1111    sfree(m);
1112    smemclr(n, plen * sizeof(*n));
1113    sfree(n);
1114}
1115
1116/*
1117 * Decrement a number.
1118 */
1119void decbn(Bignum bn)
1120{
1121    int i = 1;
1122    while (i < (int)bn[0] && bn[i] == 0)
1123        bn[i++] = BIGNUM_INT_MASK;
1124    bn[i]--;
1125}
1126
1127Bignum bignum_from_bytes(const unsigned char *data, int nbytes)
1128{
1129    Bignum result;
1130    int w, i;
1131
1132    assert(nbytes >= 0 && nbytes < INT_MAX/8);
1133
1134    w = (nbytes + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES; /* bytes->words */
1135
1136    result = newbn(w);
1137    for (i = 1; i <= w; i++)
1138        result[i] = 0;
1139    for (i = nbytes; i--;) {
1140        unsigned char byte = *data++;
1141        result[1 + i / BIGNUM_INT_BYTES] |=
1142            (BignumInt)byte << (8*i % BIGNUM_INT_BITS);
1143    }
1144
1145    while (result[0] > 1 && result[result[0]] == 0)
1146        result[0]--;
1147    return result;
1148}
1149
1150Bignum bignum_from_bytes_le(const unsigned char *data, int nbytes)
1151{
1152    Bignum result;
1153    int w, i;
1154
1155    assert(nbytes >= 0 && nbytes < INT_MAX/8);
1156
1157    w = (nbytes + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES; /* bytes->words */
1158
1159    result = newbn(w);
1160    for (i = 1; i <= w; i++)
1161        result[i] = 0;
1162    for (i = 0; i < nbytes; ++i) {
1163        unsigned char byte = *data++;
1164        result[1 + i / BIGNUM_INT_BYTES] |=
1165            (BignumInt)byte << (8*i % BIGNUM_INT_BITS);
1166    }
1167
1168    while (result[0] > 1 && result[result[0]] == 0)
1169        result[0]--;
1170    return result;
1171}
1172
1173Bignum bignum_from_decimal(const char *decimal)
1174{
1175    Bignum result = copybn(Zero);
1176
1177    while (*decimal) {
1178        Bignum tmp, tmp2;
1179
1180        if (!isdigit((unsigned char)*decimal)) {
1181            freebn(result);
1182            return 0;
1183        }
1184
1185        tmp = bigmul(result, Ten);
1186        tmp2 = bignum_from_long(*decimal - '0');
1187        result = bigadd(tmp, tmp2);
1188        freebn(tmp);
1189        freebn(tmp2);
1190
1191        decimal++;
1192    }
1193
1194    return result;
1195}
1196
1197Bignum bignum_random_in_range(const Bignum lower, const Bignum upper)
1198{
1199    Bignum ret = NULL;
1200    unsigned char *bytes;
1201    int upper_len = bignum_bitcount(upper);
1202    int upper_bytes = upper_len / 8;
1203    int upper_bits = upper_len % 8;
1204    if (upper_bits) ++upper_bytes;
1205
1206    bytes = snewn(upper_bytes, unsigned char);
1207    do {
1208        int i;
1209
1210        if (ret) freebn(ret);
1211
1212        for (i = 0; i < upper_bytes; ++i)
1213        {
1214            bytes[i] = (unsigned char)random_byte();
1215        }
1216        /* Mask the top to reduce failure rate to 50/50 */
1217        if (upper_bits)
1218        {
1219            bytes[i - 1] &= 0xFF >> (8 - upper_bits);
1220        }
1221
1222        ret = bignum_from_bytes(bytes, upper_bytes);
1223    } while (bignum_cmp(ret, lower) < 0 || bignum_cmp(ret, upper) > 0);
1224    smemclr(bytes, upper_bytes);
1225    sfree(bytes);
1226
1227    return ret;
1228}
1229
1230/*
1231 * Read an SSH-1-format bignum from a data buffer. Return the number
1232 * of bytes consumed, or -1 if there wasn't enough data.
1233 */
1234int ssh1_read_bignum(const unsigned char *data, int len, Bignum * result)
1235{
1236    const unsigned char *p = data;
1237    int i;
1238    int w, b;
1239
1240    if (len < 2)
1241        return -1;
1242
1243    w = 0;
1244    for (i = 0; i < 2; i++)
1245        w = (w << 8) + *p++;
1246    b = (w + 7) / 8;                   /* bits -> bytes */
1247
1248    if (len < b+2)
1249        return -1;
1250
1251    if (!result)                       /* just return length */
1252        return b + 2;
1253
1254    *result = bignum_from_bytes(p, b);
1255
1256    return p + b - data;
1257}
1258
1259/*
1260 * Return the bit count of a bignum, for SSH-1 encoding.
1261 */
1262int bignum_bitcount(Bignum bn)
1263{
1264    int bitcount = bn[0] * BIGNUM_INT_BITS - 1;
1265    while (bitcount >= 0
1266           && (bn[bitcount / BIGNUM_INT_BITS + 1] >> (bitcount % BIGNUM_INT_BITS)) == 0) bitcount--;
1267    return bitcount + 1;
1268}
1269
1270/*
1271 * Return the byte length of a bignum when SSH-1 encoded.
1272 */
1273int ssh1_bignum_length(Bignum bn)
1274{
1275    return 2 + (bignum_bitcount(bn) + 7) / 8;
1276}
1277
1278/*
1279 * Return the byte length of a bignum when SSH-2 encoded.
1280 */
1281int ssh2_bignum_length(Bignum bn)
1282{
1283    return 4 + (bignum_bitcount(bn) + 8) / 8;
1284}
1285
1286/*
1287 * Return a byte from a bignum; 0 is least significant, etc.
1288 */
1289int bignum_byte(Bignum bn, int i)
1290{
1291    if (i < 0 || i >= (int)(BIGNUM_INT_BYTES * bn[0]))
1292        return 0;                      /* beyond the end */
1293    else
1294        return (bn[i / BIGNUM_INT_BYTES + 1] >>
1295                ((i % BIGNUM_INT_BYTES)*8)) & 0xFF;
1296}
1297
1298/*
1299 * Return a bit from a bignum; 0 is least significant, etc.
1300 */
1301int bignum_bit(Bignum bn, int i)
1302{
1303    if (i < 0 || i >= (int)(BIGNUM_INT_BITS * bn[0]))
1304        return 0;                      /* beyond the end */
1305    else
1306        return (bn[i / BIGNUM_INT_BITS + 1] >> (i % BIGNUM_INT_BITS)) & 1;
1307}
1308
1309/*
1310 * Set a bit in a bignum; 0 is least significant, etc.
1311 */
1312void bignum_set_bit(Bignum bn, int bitnum, int value)
1313{
1314    if (bitnum < 0 || bitnum >= (int)(BIGNUM_INT_BITS * bn[0]))
1315        abort();                       /* beyond the end */
1316    else {
1317        int v = bitnum / BIGNUM_INT_BITS + 1;
1318        BignumInt mask = (BignumInt)1 << (bitnum % BIGNUM_INT_BITS);
1319        if (value)
1320            bn[v] |= mask;
1321        else
1322            bn[v] &= ~mask;
1323    }
1324}
1325
1326/*
1327 * Write a SSH-1-format bignum into a buffer. It is assumed the
1328 * buffer is big enough. Returns the number of bytes used.
1329 */
1330int ssh1_write_bignum(void *data, Bignum bn)
1331{
1332    unsigned char *p = data;
1333    int len = ssh1_bignum_length(bn);
1334    int i;
1335    int bitc = bignum_bitcount(bn);
1336
1337    *p++ = (bitc >> 8) & 0xFF;
1338    *p++ = (bitc) & 0xFF;
1339    for (i = len - 2; i--;)
1340        *p++ = bignum_byte(bn, i);
1341    return len;
1342}
1343
1344/*
1345 * Compare two bignums. Returns like strcmp.
1346 */
1347int bignum_cmp(Bignum a, Bignum b)
1348{
1349    int amax = a[0], bmax = b[0];
1350    int i;
1351
1352    /* Annoyingly we have two representations of zero */
1353    if (amax == 1 && a[amax] == 0)
1354        amax = 0;
1355    if (bmax == 1 && b[bmax] == 0)
1356        bmax = 0;
1357
1358    assert(amax == 0 || a[amax] != 0);
1359    assert(bmax == 0 || b[bmax] != 0);
1360
1361    i = (amax > bmax ? amax : bmax);
1362    while (i) {
1363        BignumInt aval = (i > amax ? 0 : a[i]);
1364        BignumInt bval = (i > bmax ? 0 : b[i]);
1365        if (aval < bval)
1366            return -1;
1367        if (aval > bval)
1368            return +1;
1369        i--;
1370    }
1371    return 0;
1372}
1373
1374/*
1375 * Right-shift one bignum to form another.
1376 */
1377Bignum bignum_rshift(Bignum a, int shift)
1378{
1379    Bignum ret;
1380    int i, shiftw, shiftb, shiftbb, bits;
1381    BignumInt ai, ai1;
1382
1383    assert(shift >= 0);
1384
1385    bits = bignum_bitcount(a) - shift;
1386    ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);
1387
1388    if (ret) {
1389        shiftw = shift / BIGNUM_INT_BITS;
1390        shiftb = shift % BIGNUM_INT_BITS;
1391        shiftbb = BIGNUM_INT_BITS - shiftb;
1392
1393        ai1 = a[shiftw + 1];
1394        for (i = 1; i <= (int)ret[0]; i++) {
1395            ai = ai1;
1396            ai1 = (i + shiftw + 1 <= (int)a[0] ? a[i + shiftw + 1] : 0);
1397            ret[i] = ((ai >> shiftb) | (ai1 << shiftbb)) & BIGNUM_INT_MASK;
1398        }
1399    }
1400
1401    return ret;
1402}
1403
1404/*
1405 * Left-shift one bignum to form another.
1406 */
1407Bignum bignum_lshift(Bignum a, int shift)
1408{
1409    Bignum ret;
1410    int bits, shiftWords, shiftBits;
1411
1412    assert(shift >= 0);
1413
1414    bits = bignum_bitcount(a) + shift;
1415    ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);
1416
1417    shiftWords = shift / BIGNUM_INT_BITS;
1418    shiftBits = shift % BIGNUM_INT_BITS;
1419
1420    if (shiftBits == 0)
1421    {
1422        memcpy(&ret[1 + shiftWords], &a[1], sizeof(BignumInt) * a[0]);
1423    }
1424    else
1425    {
1426        int i;
1427        BignumInt carry = 0;
1428
1429        /* Remember that Bignum[0] is length, so add 1 */
1430        for (i = shiftWords + 1; i < ((int)a[0]) + shiftWords + 1; ++i)
1431        {
1432            BignumInt from = a[i - shiftWords];
1433            ret[i] = (from << shiftBits) | carry;
1434            carry = from >> (BIGNUM_INT_BITS - shiftBits);
1435        }
1436        if (carry) ret[i] = carry;
1437    }
1438
1439    return ret;
1440}
1441
1442/*
1443 * Non-modular multiplication and addition.
1444 */
1445Bignum bigmuladd(Bignum a, Bignum b, Bignum addend)
1446{
1447    int alen = a[0], blen = b[0];
1448    int mlen = (alen > blen ? alen : blen);
1449    int rlen, i, maxspot;
1450    int wslen;
1451    BignumInt *workspace;
1452    Bignum ret;
1453
1454    /* mlen space for a, mlen space for b, 2*mlen for result,
1455     * plus scratch space for multiplication */
1456    wslen = mlen * 4 + mul_compute_scratch(mlen);
1457    workspace = snewn(wslen, BignumInt);
1458    for (i = 0; i < mlen; i++) {
1459        workspace[0 * mlen + i] = (mlen - i <= (int)a[0] ? a[mlen - i] : 0);
1460        workspace[1 * mlen + i] = (mlen - i <= (int)b[0] ? b[mlen - i] : 0);
1461    }
1462
1463    internal_mul(workspace + 0 * mlen, workspace + 1 * mlen,
1464                 workspace + 2 * mlen, mlen, workspace + 4 * mlen);
1465
1466    /* now just copy the result back */
1467    rlen = alen + blen + 1;
1468    if (addend && rlen <= (int)addend[0])
1469        rlen = addend[0] + 1;
1470    ret = newbn(rlen);
1471    maxspot = 0;
1472    for (i = 1; i <= (int)ret[0]; i++) {
1473        ret[i] = (i <= 2 * mlen ? workspace[4 * mlen - i] : 0);
1474        if (ret[i] != 0)
1475            maxspot = i;
1476    }
1477    ret[0] = maxspot;
1478
1479    /* now add in the addend, if any */
1480    if (addend) {
1481        BignumDblInt carry = 0;
1482        for (i = 1; i <= rlen; i++) {
1483            carry += (i <= (int)ret[0] ? ret[i] : 0);
1484            carry += (i <= (int)addend[0] ? addend[i] : 0);
1485            ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
1486            carry >>= BIGNUM_INT_BITS;
1487            if (ret[i] != 0 && i > maxspot)
1488                maxspot = i;
1489        }
1490    }
1491    ret[0] = maxspot;
1492
1493    smemclr(workspace, wslen * sizeof(*workspace));
1494    sfree(workspace);
1495    return ret;
1496}
1497
1498/*
1499 * Non-modular multiplication.
1500 */
1501Bignum bigmul(Bignum a, Bignum b)
1502{
1503    return bigmuladd(a, b, NULL);
1504}
1505
1506/*
1507 * Simple addition.
1508 */
1509Bignum bigadd(Bignum a, Bignum b)
1510{
1511    int alen = a[0], blen = b[0];
1512    int rlen = (alen > blen ? alen : blen) + 1;
1513    int i, maxspot;
1514    Bignum ret;
1515    BignumDblInt carry;
1516
1517    ret = newbn(rlen);
1518
1519    carry = 0;
1520    maxspot = 0;
1521    for (i = 1; i <= rlen; i++) {
1522        carry += (i <= (int)a[0] ? a[i] : 0);
1523        carry += (i <= (int)b[0] ? b[i] : 0);
1524        ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
1525        carry >>= BIGNUM_INT_BITS;
1526        if (ret[i] != 0 && i > maxspot)
1527            maxspot = i;
1528    }
1529    ret[0] = maxspot;
1530
1531    return ret;
1532}
1533
1534/*
1535 * Subtraction. Returns a-b, or NULL if the result would come out
1536 * negative (recall that this entire bignum module only handles
1537 * positive numbers).
1538 */
1539Bignum bigsub(Bignum a, Bignum b)
1540{
1541    int alen = a[0], blen = b[0];
1542    int rlen = (alen > blen ? alen : blen);
1543    int i, maxspot;
1544    Bignum ret;
1545    BignumDblInt carry;
1546
1547    ret = newbn(rlen);
1548
1549    carry = 1;
1550    maxspot = 0;
1551    for (i = 1; i <= rlen; i++) {
1552        carry += (i <= (int)a[0] ? a[i] : 0);
1553        carry += (i <= (int)b[0] ? b[i] ^ BIGNUM_INT_MASK : BIGNUM_INT_MASK);
1554        ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
1555        carry >>= BIGNUM_INT_BITS;
1556        if (ret[i] != 0 && i > maxspot)
1557            maxspot = i;
1558    }
1559    ret[0] = maxspot;
1560
1561    if (!carry) {
1562        freebn(ret);
1563        return NULL;
1564    }
1565
1566    return ret;
1567}
1568
1569/*
1570 * Create a bignum which is the bitmask covering another one. That
1571 * is, the smallest integer which is >= N and is also one less than
1572 * a power of two.
1573 */
1574Bignum bignum_bitmask(Bignum n)
1575{
1576    Bignum ret = copybn(n);
1577    int i;
1578    BignumInt j;
1579
1580    i = ret[0];
1581    while (n[i] == 0 && i > 0)
1582        i--;
1583    if (i <= 0)
1584        return ret;                    /* input was zero */
1585    j = 1;
1586    while (j < n[i])
1587        j = 2 * j + 1;
1588    ret[i] = j;
1589    while (--i > 0)
1590        ret[i] = BIGNUM_INT_MASK;
1591    return ret;
1592}
1593
1594/*
1595 * Convert a (max 32-bit) long into a bignum.
1596 */
1597Bignum bignum_from_long(unsigned long nn)
1598{
1599    Bignum ret;
1600    BignumDblInt n = nn;
1601
1602    ret = newbn(3);
1603    ret[1] = (BignumInt)(n & BIGNUM_INT_MASK);
1604    ret[2] = (BignumInt)((n >> BIGNUM_INT_BITS) & BIGNUM_INT_MASK);
1605    ret[3] = 0;
1606    ret[0] = (ret[2]  ? 2 : 1);
1607    return ret;
1608}
1609
1610/*
1611 * Add a long to a bignum.
1612 */
1613Bignum bignum_add_long(Bignum number, unsigned long addendx)
1614{
1615    Bignum ret = newbn(number[0] + 1);
1616    int i, maxspot = 0;
1617    BignumDblInt carry = 0, addend = addendx;
1618
1619    for (i = 1; i <= (int)ret[0]; i++) {
1620        carry += addend & BIGNUM_INT_MASK;
1621        carry += (i <= (int)number[0] ? number[i] : 0);
1622        addend >>= BIGNUM_INT_BITS;
1623        ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
1624        carry >>= BIGNUM_INT_BITS;
1625        if (ret[i] != 0)
1626            maxspot = i;
1627    }
1628    ret[0] = maxspot;
1629    return ret;
1630}
1631
1632/*
1633 * Compute the residue of a bignum, modulo a (max 16-bit) short.
1634 */
1635unsigned short bignum_mod_short(Bignum number, unsigned short modulus)
1636{
1637    BignumDblInt mod, r;
1638    int i;
1639
1640    r = 0;
1641    mod = modulus;
1642    for (i = number[0]; i > 0; i--)
1643        r = (r * (BIGNUM_TOP_BIT % mod) * 2 + number[i] % mod) % mod;
1644    return (unsigned short) r;
1645}
1646
1647#ifdef DEBUG
1648void diagbn(char *prefix, Bignum md)
1649{
1650    int i, nibbles, morenibbles;
1651    static const char hex[] = "0123456789ABCDEF";
1652
1653    debug(("%s0x", prefix ? prefix : ""));
1654
1655    nibbles = (3 + bignum_bitcount(md)) / 4;
1656    if (nibbles < 1)
1657        nibbles = 1;
1658    morenibbles = 4 * md[0] - nibbles;
1659    for (i = 0; i < morenibbles; i++)
1660        debug(("-"));
1661    for (i = nibbles; i--;)
1662        debug(("%c",
1663               hex[(bignum_byte(md, i / 2) >> (4 * (i % 2))) & 0xF]));
1664
1665    if (prefix)
1666        debug(("\n"));
1667}
1668#endif
1669
1670/*
1671 * Simple division.
1672 */
1673Bignum bigdiv(Bignum a, Bignum b)
1674{
1675    Bignum q = newbn(a[0]);
1676    bigdivmod(a, b, NULL, q);
1677    while (q[0] > 1 && q[q[0]] == 0)
1678        q[0]--;
1679    return q;
1680}
1681
1682/*
1683 * Simple remainder.
1684 */
1685Bignum bigmod(Bignum a, Bignum b)
1686{
1687    Bignum r = newbn(b[0]);
1688    bigdivmod(a, b, r, NULL);
1689    while (r[0] > 1 && r[r[0]] == 0)
1690        r[0]--;
1691    return r;
1692}
1693
1694/*
1695 * Greatest common divisor.
1696 */
1697Bignum biggcd(Bignum av, Bignum bv)
1698{
1699    Bignum a = copybn(av);
1700    Bignum b = copybn(bv);
1701
1702    while (bignum_cmp(b, Zero) != 0) {
1703        Bignum t = newbn(b[0]);
1704        bigdivmod(a, b, t, NULL);
1705        while (t[0] > 1 && t[t[0]] == 0)
1706            t[0]--;
1707        freebn(a);
1708        a = b;
1709        b = t;
1710    }
1711
1712    freebn(b);
1713    return a;
1714}
1715
1716/*
1717 * Modular inverse, using Euclid's extended algorithm.
1718 */
1719Bignum modinv(Bignum number, Bignum modulus)
1720{
1721    Bignum a = copybn(modulus);
1722    Bignum b = copybn(number);
1723    Bignum xp = copybn(Zero);
1724    Bignum x = copybn(One);
1725    int sign = +1;
1726
1727    assert(number[number[0]] != 0);
1728    assert(modulus[modulus[0]] != 0);
1729
1730    while (bignum_cmp(b, One) != 0) {
1731        Bignum t, q;
1732
1733        if (bignum_cmp(b, Zero) == 0) {
1734            /*
1735             * Found a common factor between the inputs, so we cannot
1736             * return a modular inverse at all.
1737             */
1738            freebn(b);
1739            freebn(a);
1740            freebn(xp);
1741            freebn(x);
1742            return NULL;
1743        }
1744
1745        t = newbn(b[0]);
1746        q = newbn(a[0]);
1747        bigdivmod(a, b, t, q);
1748        while (t[0] > 1 && t[t[0]] == 0)
1749            t[0]--;
1750        while (q[0] > 1 && q[q[0]] == 0)
1751            q[0]--;
1752        freebn(a);
1753        a = b;
1754        b = t;
1755        t = xp;
1756        xp = x;
1757        x = bigmuladd(q, xp, t);
1758        sign = -sign;
1759        freebn(t);
1760        freebn(q);
1761    }
1762
1763    freebn(b);
1764    freebn(a);
1765    freebn(xp);
1766
1767    /* now we know that sign * x == 1, and that x < modulus */
1768    if (sign < 0) {
1769        /* set a new x to be modulus - x */
1770        Bignum newx = newbn(modulus[0]);
1771        BignumInt carry = 0;
1772        int maxspot = 1;
1773        int i;
1774
1775        for (i = 1; i <= (int)newx[0]; i++) {
1776            BignumInt aword = (i <= (int)modulus[0] ? modulus[i] : 0);
1777            BignumInt bword = (i <= (int)x[0] ? x[i] : 0);
1778            newx[i] = aword - bword - carry;
1779            bword = ~bword;
1780            carry = carry ? (newx[i] >= bword) : (newx[i] > bword);
1781            if (newx[i] != 0)
1782                maxspot = i;
1783        }
1784        newx[0] = maxspot;
1785        freebn(x);
1786        x = newx;
1787    }
1788
1789    /* and return. */
1790    return x;
1791}
1792
1793/*
1794 * Render a bignum into decimal. Return a malloced string holding
1795 * the decimal representation.
1796 */
1797char *bignum_decimal(Bignum x)
1798{
1799    int ndigits, ndigit;
1800    int i, iszero;
1801    BignumDblInt carry;
1802    char *ret;
1803    BignumInt *workspace;
1804
1805    /*
1806     * First, estimate the number of digits. Since log(10)/log(2)
1807     * is just greater than 93/28 (the joys of continued fraction
1808     * approximations...) we know that for every 93 bits, we need
1809     * at most 28 digits. This will tell us how much to malloc.
1810     *
1811     * Formally: if x has i bits, that means x is strictly less
1812     * than 2^i. Since 2 is less than 10^(28/93), this is less than
1813     * 10^(28i/93). We need an integer power of ten, so we must
1814     * round up (rounding down might make it less than x again).
1815     * Therefore if we multiply the bit count by 28/93, rounding
1816     * up, we will have enough digits.
1817     *
1818     * i=0 (i.e., x=0) is an irritating special case.
1819     */
1820    i = bignum_bitcount(x);
1821    if (!i)
1822        ndigits = 1;                   /* x = 0 */
1823    else
1824        ndigits = (28 * i + 92) / 93;  /* multiply by 28/93 and round up */
1825    ndigits++;                         /* allow for trailing \0 */
1826    ret = snewn(ndigits, char);
1827
1828    /*
1829     * Now allocate some workspace to hold the binary form as we
1830     * repeatedly divide it by ten. Initialise this to the
1831     * big-endian form of the number.
1832     */
1833    workspace = snewn(x[0], BignumInt);
1834    for (i = 0; i < (int)x[0]; i++)
1835        workspace[i] = x[x[0] - i];
1836
1837    /*
1838     * Next, write the decimal number starting with the last digit.
1839     * We use ordinary short division, dividing 10 into the
1840     * workspace.
1841     */
1842    ndigit = ndigits - 1;
1843    ret[ndigit] = '\0';
1844    do {
1845        iszero = 1;
1846        carry = 0;
1847        for (i = 0; i < (int)x[0]; i++) {
1848            carry = (carry << BIGNUM_INT_BITS) + workspace[i];
1849            workspace[i] = (BignumInt) (carry / 10);
1850            if (workspace[i])
1851                iszero = 0;
1852            carry %= 10;
1853        }
1854        ret[--ndigit] = (char) (carry + '0');
1855    } while (!iszero);
1856
1857    /*
1858     * There's a chance we've fallen short of the start of the
1859     * string. Correct if so.
1860     */
1861    if (ndigit > 0)
1862        memmove(ret, ret + ndigit, ndigits - ndigit);
1863
1864    /*
1865     * Done.
1866     */
1867    smemclr(workspace, x[0] * sizeof(*workspace));
1868    sfree(workspace);
1869    return ret;
1870}
1871
1872#ifdef TESTBN
1873
1874#include <stdio.h>
1875#include <stdlib.h>
1876#include <ctype.h>
1877
1878/*
1879 * gcc -Wall -g -O0 -DTESTBN -o testbn sshbn.c misc.c conf.c tree234.c unix/uxmisc.c -I. -I unix -I charset
1880 *
1881 * Then feed to this program's standard input the output of
1882 * testdata/bignum.py .
1883 */
1884
1885void modalfatalbox(const char *p, ...)
1886{
1887    va_list ap;
1888    fprintf(stderr, "FATAL ERROR: ");
1889    va_start(ap, p);
1890    vfprintf(stderr, p, ap);
1891    va_end(ap);
1892    fputc('\n', stderr);
1893    exit(1);
1894}
1895
1896int random_byte(void)
1897{
1898    modalfatalbox("random_byte called in testbn");
1899    return 0;
1900}
1901
1902#define fromxdigit(c) ( (c)>'9' ? ((c)&0xDF) - 'A' + 10 : (c) - '0' )
1903
1904int main(int argc, char **argv)
1905{
1906    char *buf;
1907    int line = 0;
1908    int passes = 0, fails = 0;
1909
1910    while ((buf = fgetline(stdin)) != NULL) {
1911        int maxlen = strlen(buf);
1912        unsigned char *data = snewn(maxlen, unsigned char);
1913        unsigned char *ptrs[5], *q;
1914        int ptrnum;
1915        char *bufp = buf;
1916
1917        line++;
1918
1919        q = data;
1920        ptrnum = 0;
1921
1922        while (*bufp && !isspace((unsigned char)*bufp))
1923            bufp++;
1924        if (bufp)
1925            *bufp++ = '\0';
1926
1927        while (*bufp) {
1928            char *start, *end;
1929            int i;
1930
1931            while (*bufp && !isxdigit((unsigned char)*bufp))
1932                bufp++;
1933            start = bufp;
1934
1935            if (!*bufp)
1936                break;
1937
1938            while (*bufp && isxdigit((unsigned char)*bufp))
1939                bufp++;
1940            end = bufp;
1941
1942            if (ptrnum >= lenof(ptrs))
1943                break;
1944            ptrs[ptrnum++] = q;
1945           
1946            for (i = -((end - start) & 1); i < end-start; i += 2) {
1947                unsigned char val = (i < 0 ? 0 : fromxdigit(start[i]));
1948                val = val * 16 + fromxdigit(start[i+1]);
1949                *q++ = val;
1950            }
1951
1952            ptrs[ptrnum] = q;
1953        }
1954
1955        if (!strcmp(buf, "mul")) {
1956            Bignum a, b, c, p;
1957
1958            if (ptrnum != 3) {
1959                printf("%d: mul with %d parameters, expected 3\n", line, ptrnum);
1960                exit(1);
1961            }
1962            a = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);
1963            b = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);
1964            c = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);
1965            p = bigmul(a, b);
1966
1967            if (bignum_cmp(c, p) == 0) {
1968                passes++;
1969            } else {
1970                char *as = bignum_decimal(a);
1971                char *bs = bignum_decimal(b);
1972                char *cs = bignum_decimal(c);
1973                char *ps = bignum_decimal(p);
1974               
1975                printf("%d: fail: %s * %s gave %s expected %s\n",
1976                       line, as, bs, ps, cs);
1977                fails++;
1978
1979                sfree(as);
1980                sfree(bs);
1981                sfree(cs);
1982                sfree(ps);
1983            }
1984            freebn(a);
1985            freebn(b);
1986            freebn(c);
1987            freebn(p);
1988        } else if (!strcmp(buf, "modmul")) {
1989            Bignum a, b, m, c, p;
1990
1991            if (ptrnum != 4) {
1992                printf("%d: modmul with %d parameters, expected 4\n",
1993                       line, ptrnum);
1994                exit(1);
1995            }
1996            a = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);
1997            b = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);
1998            m = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);
1999            c = bignum_from_bytes(ptrs[3], ptrs[4]-ptrs[3]);
2000            p = modmul(a, b, m);
2001
2002            if (bignum_cmp(c, p) == 0) {
2003                passes++;
2004            } else {
2005                char *as = bignum_decimal(a);
2006                char *bs = bignum_decimal(b);
2007                char *ms = bignum_decimal(m);
2008                char *cs = bignum_decimal(c);
2009                char *ps = bignum_decimal(p);
2010               
2011                printf("%d: fail: %s * %s mod %s gave %s expected %s\n",
2012                       line, as, bs, ms, ps, cs);
2013                fails++;
2014
2015                sfree(as);
2016                sfree(bs);
2017                sfree(ms);
2018                sfree(cs);
2019                sfree(ps);
2020            }
2021            freebn(a);
2022            freebn(b);
2023            freebn(m);
2024            freebn(c);
2025            freebn(p);
2026        } else if (!strcmp(buf, "pow")) {
2027            Bignum base, expt, modulus, expected, answer;
2028
2029            if (ptrnum != 4) {
2030                printf("%d: mul with %d parameters, expected 4\n", line, ptrnum);
2031                exit(1);
2032            }
2033
2034            base = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);
2035            expt = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);
2036            modulus = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);
2037            expected = bignum_from_bytes(ptrs[3], ptrs[4]-ptrs[3]);
2038            answer = modpow(base, expt, modulus);
2039
2040            if (bignum_cmp(expected, answer) == 0) {
2041                passes++;
2042            } else {
2043                char *as = bignum_decimal(base);
2044                char *bs = bignum_decimal(expt);
2045                char *cs = bignum_decimal(modulus);
2046                char *ds = bignum_decimal(answer);
2047                char *ps = bignum_decimal(expected);
2048               
2049                printf("%d: fail: %s ^ %s mod %s gave %s expected %s\n",
2050                       line, as, bs, cs, ds, ps);
2051                fails++;
2052
2053                sfree(as);
2054                sfree(bs);
2055                sfree(cs);
2056                sfree(ds);
2057                sfree(ps);
2058            }
2059            freebn(base);
2060            freebn(expt);
2061            freebn(modulus);
2062            freebn(expected);
2063            freebn(answer);
2064        } else {
2065            printf("%d: unrecognised test keyword: '%s'\n", line, buf);
2066            exit(1);
2067        }
2068
2069        sfree(buf);
2070        sfree(data);
2071    }
2072
2073    printf("passed %d failed %d total %d\n", passes, fails, passes+fails);
2074    return fails != 0;
2075}
2076
2077#endif
Note: See TracBrowser for help on using the repository browser.