Upgrade AXTLS import to version 1.1.5-a
[people/adir/gpxe.git] / src / crypto / axtls / bigint.c
1 /*
2  *  Copyright(C) 2006 Cameron Rich
3  *
4  *  This library is free software; you can redistribute it and/or modify
5  *  it under the terms of the GNU Lesser General Public License as published by
6  *  the Free Software Foundation; either version 2.1 of the License, or
7  *  (at your option) any later version.
8  *
9  *  This library is distributed in the hope that it will be useful,
10  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  *  GNU Lesser General Public License for more details.
13  *
14  *  You should have received a copy of the GNU Lesser General Public License
15  *  along with this library; if not, write to the Free Software
16  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17  */
18
19 /**
20  * @defgroup bigint_api Big Integer API
21  * @brief The bigint implementation as used by the axTLS project.
22  *
23  * The bigint library is for RSA encryption/decryption as well as signing.
24  * This code tries to minimise use of malloc/free by maintaining a small 
25  * cache. A bigint context may maintain state by being made "permanent". 
26  * It be be later released with a bi_depermanent() and bi_free() call.
27  *
28  * It supports the following reduction techniques:
29  * - Classical
30  * - Barrett
31  * - Montgomery
32  *
33  * It also implements the following:
34  * - Karatsuba multiplication
35  * - Squaring
36  * - Sliding window exponentiation
37  * - Chinese Remainder Theorem (implemented in rsa.c).
38  *
39  * All the algorithms used are pretty standard, and designed for different
40  * data bus sizes. Negative numbers are not dealt with at all, so a subtraction
41  * may need to be tested for negativity.
42  *
43  * This library steals some ideas from Jef Poskanzer
44  * <http://cs.marlboro.edu/term/cs-fall02/algorithms/crypto/RSA/bigint>
45  * and GMP <http://www.swox.com/gmp>. It gets most of its implementation
46  * detail from "The Handbook of Applied Cryptography"
47  * <http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf>
48  * @{
49  */
50
51 #include <stdlib.h>
52 #include <limits.h>
53 #include <string.h>
54 #include <stdio.h>
55 #include <time.h>
56 #include "bigint.h"
57 #include "crypto.h"
58
59 static bigint *bi_int_multiply(BI_CTX *ctx, bigint *bi, comp i);
60 static bigint *bi_int_divide(BI_CTX *ctx, bigint *biR, comp denom);
61 static bigint *alloc(BI_CTX *ctx, int size);
62 static bigint *trim(bigint *bi);
63 static void more_comps(bigint *bi, int n);
64 #if defined(CONFIG_BIGINT_KARATSUBA) || defined(CONFIG_BIGINT_BARRETT) || \
65     defined(CONFIG_BIGINT_MONTGOMERY)
66 static bigint *comp_right_shift(bigint *biR, int num_shifts);
67 static bigint *comp_left_shift(bigint *biR, int num_shifts);
68 #endif
69
70 #ifdef CONFIG_BIGINT_CHECK_ON
71 static void check(const bigint *bi);
72 #endif
73
74 /**
75  * @brief Start a new bigint context.
76  * @return A bigint context.
77  */
78 BI_CTX *bi_initialize(void)
79 {
80     /* calloc() sets everything to zero */
81     BI_CTX *ctx = (BI_CTX *)calloc(1, sizeof(BI_CTX));
82    
83     /* the radix */
84     ctx->bi_radix = alloc(ctx, 2); 
85     ctx->bi_radix->comps[0] = 0;
86     ctx->bi_radix->comps[1] = 1;
87     bi_permanent(ctx->bi_radix);
88     return ctx;
89 }
90
91 /**
92  * @brief Close the bigint context and free any resources.
93  *
94  * Free up any used memory - a check is done if all objects were not 
95  * properly freed.
96  * @param ctx [in]   The bigint session context.
97  */
98 void bi_terminate(BI_CTX *ctx)
99 {
100     bigint *p, *pn;
101
102     bi_depermanent(ctx->bi_radix); 
103     bi_free(ctx, ctx->bi_radix);
104
105     if (ctx->active_count != 0)
106     {
107 #ifdef CONFIG_SSL_FULL_MODE
108         printf("bi_terminate: there were %d un-freed bigints\n",
109                        ctx->active_count);
110 #endif
111         abort();
112     }
113
114     for (p = ctx->free_list; p != NULL; p = pn)
115     {
116         pn = p->next;
117         free(p->comps);
118         free(p);
119     }
120
121     free(ctx);
122 }
123
124 /**
125  * @brief Increment the number of references to this object. 
126  * It does not do a full copy.
127  * @param bi [in]   The bigint to copy.
128  * @return A reference to the same bigint.
129  */
130 bigint *bi_copy(bigint *bi)
131 {
132     check(bi);
133     if (bi->refs != PERMANENT)
134         bi->refs++;
135     return bi;
136 }
137
138 /**
139  * @brief Simply make a bigint object "unfreeable" if bi_free() is called on it.
140  *
141  * For this object to be freed, bi_depermanent() must be called.
142  * @param bi [in]   The bigint to be made permanent.
143  */
144 void bi_permanent(bigint *bi)
145 {
146     check(bi);
147     if (bi->refs != 1)
148     {
149 #ifdef CONFIG_SSL_FULL_MODE
150         printf("bi_permanent: refs was not 1\n");
151 #endif
152         abort();
153     }
154
155     bi->refs = PERMANENT;
156 }
157
158 /**
159  * @brief Take a permanent object and make it eligible for freedom.
160  * @param bi [in]   The bigint to be made back to temporary.
161  */
162 void bi_depermanent(bigint *bi)
163 {
164     check(bi);
165     if (bi->refs != PERMANENT)
166     {
167 #ifdef CONFIG_SSL_FULL_MODE
168         printf("bi_depermanent: bigint was not permanent\n");
169 #endif
170         abort();
171     }
172
173     bi->refs = 1;
174 }
175
176 /**
177  * @brief Free a bigint object so it can be used again. 
178  *
179  * The memory itself it not actually freed, just tagged as being available 
180  * @param ctx [in]   The bigint session context.
181  * @param bi [in]    The bigint to be freed.
182  */
183 void bi_free(BI_CTX *ctx, bigint *bi)
184 {
185     check(bi);
186     if (bi->refs == PERMANENT)
187     {
188         return;
189     }
190
191     if (--bi->refs > 0)
192     {
193         return;
194     }
195
196     bi->next = ctx->free_list;
197     ctx->free_list = bi;
198     ctx->free_count++;
199
200     if (--ctx->active_count < 0)
201     {
202 #ifdef CONFIG_SSL_FULL_MODE
203         printf("bi_free: active_count went negative "
204                 "- double-freed bigint?\n");
205 #endif
206         abort();
207     }
208 }
209
210 /**
211  * @brief Convert an (unsigned) integer into a bigint.
212  * @param ctx [in]   The bigint session context.
213  * @param i [in]     The (unsigned) integer to be converted.
214  * 
215  */
216 bigint *int_to_bi(BI_CTX *ctx, comp i)
217 {
218     bigint *biR = alloc(ctx, 1);
219     biR->comps[0] = i;
220     return biR;
221 }
222
223 /**
224  * @brief Do a full copy of the bigint object.
225  * @param ctx [in]   The bigint session context.
226  * @param bi  [in]   The bigint object to be copied.
227  */
228 bigint *bi_clone(BI_CTX *ctx, const bigint *bi)
229 {
230     bigint *biR = alloc(ctx, bi->size);
231     check(bi);
232     memcpy(biR->comps, bi->comps, bi->size*COMP_BYTE_SIZE);
233     return biR;
234 }
235
236 /**
237  * @brief Perform an addition operation between two bigints.
238  * @param ctx [in]  The bigint session context.
239  * @param bia [in]  A bigint.
240  * @param bib [in]  Another bigint.
241  * @return The result of the addition.
242  */
243 bigint *bi_add(BI_CTX *ctx, bigint *bia, bigint *bib)
244 {
245     int n;
246     comp carry = 0;
247     comp *pa, *pb;
248
249     check(bia);
250     check(bib);
251
252     n = max(bia->size, bib->size);
253     more_comps(bia, n+1);
254     more_comps(bib, n);
255     pa = bia->comps;
256     pb = bib->comps;
257
258     do
259     {
260         comp  sl, rl, cy1;
261         sl = *pa + *pb++;
262         rl = sl + carry;
263         cy1 = sl < *pa;
264         carry = cy1 | (rl < sl);
265         *pa++ = rl;
266     } while (--n != 0);
267
268     *pa = carry;                  /* do overflow */
269     bi_free(ctx, bib);
270     return trim(bia);
271 }
272
273 /**
274  * @brief Perform a subtraction operation between two bigints.
275  * @param ctx [in]  The bigint session context.
276  * @param bia [in]  A bigint.
277  * @param bib [in]  Another bigint.
278  * @param is_negative [out] If defined, indicates that the result was negative.
279  * is_negative may be null.
280  * @return The result of the subtraction. The result is always positive.
281  */
282 bigint *bi_subtract(BI_CTX *ctx, 
283         bigint *bia, bigint *bib, int *is_negative)
284 {
285     int n = bia->size;
286     comp *pa, *pb, carry = 0;
287
288     check(bia);
289     check(bib);
290
291     more_comps(bib, n);
292     pa = bia->comps;
293     pb = bib->comps;
294
295     do 
296     {
297         comp sl, rl, cy1;
298         sl = *pa - *pb++;
299         rl = sl - carry;
300         cy1 = sl > *pa;
301         carry = cy1 | (rl > sl);
302         *pa++ = rl;
303     } while (--n != 0);
304
305     if (is_negative)    /* indicate a negative result */
306     {
307         *is_negative = carry;
308     }
309
310     bi_free(ctx, trim(bib));    /* put bib back to the way it was */
311     return trim(bia);
312 }
313
314 /**
315  * Perform a multiply between a bigint an an (unsigned) integer
316  */
317 static bigint *bi_int_multiply(BI_CTX *ctx, bigint *bia, comp b)
318 {
319     int j = 0, n = bia->size;
320     bigint *biR = alloc(ctx, n + 1);
321     comp carry = 0;
322     comp *r = biR->comps;
323     comp *a = bia->comps;
324
325     check(bia);
326
327     /* clear things to start with */
328     memset(r, 0, ((n+1)*COMP_BYTE_SIZE));
329
330     do
331     {
332         long_comp tmp = *r + (long_comp)a[j]*b + carry;
333         *r++ = (comp)tmp;              /* downsize */
334         carry = (comp)(tmp >> COMP_BIT_SIZE);
335     } while (++j < n);
336
337     *r = carry;
338     bi_free(ctx, bia);
339     return trim(biR);
340 }
341
342 /**
343  * @brief Does both division and modulo calculations. 
344  *
345  * Used extensively when doing classical reduction.
346  * @param ctx [in]  The bigint session context.
347  * @param u [in]    A bigint which is the numerator.
348  * @param v [in]    Either the denominator or the modulus depending on the mode.
349  * @param is_mod [n] Determines if this is a normal division (0) or a reduction
350  * (1).
351  * @return  The result of the division/reduction.
352  */
353 bigint *bi_divide(BI_CTX *ctx, bigint *u, bigint *v, int is_mod)
354 {
355     int n = v->size, m = u->size-n;
356     int j = 0, orig_u_size = u->size;
357     uint8_t mod_offset = ctx->mod_offset;
358     comp d;
359     bigint *quotient, *tmp_u;
360     comp q_dash;
361
362     check(u);
363     check(v);
364
365     /* if doing reduction and we are < mod, then return mod */
366     if (is_mod && bi_compare(v, u) > 0)
367     {
368         bi_free(ctx, v);
369         return u;
370     }
371
372     quotient = alloc(ctx, m+1);
373     tmp_u = alloc(ctx, n+1);
374     v = trim(v);        /* make sure we have no leading 0's */
375     d = (comp)((long_comp)COMP_RADIX/(V1+1));
376
377     /* clear things to start with */
378     memset(quotient->comps, 0, ((quotient->size)*COMP_BYTE_SIZE));
379
380     /* normalise */
381     if (d > 1)
382     {
383         u = bi_int_multiply(ctx, u, d);
384
385         if (is_mod)
386         {
387             v = ctx->bi_normalised_mod[mod_offset];
388         }
389         else
390         {
391             v = bi_int_multiply(ctx, v, d);
392         }
393     }
394
395     if (orig_u_size == u->size)  /* new digit position u0 */
396     {
397         more_comps(u, orig_u_size + 1);
398     }
399
400     do
401     {
402         /* get a temporary short version of u */
403         memcpy(tmp_u->comps, &u->comps[u->size-n-1-j], (n+1)*COMP_BYTE_SIZE);
404
405         /* calculate q' */
406         if (U(0) == V1)
407         {
408             q_dash = COMP_RADIX-1;
409         }
410         else
411         {
412             q_dash = (comp)(((long_comp)U(0)*COMP_RADIX + U(1))/V1);
413         }
414
415         if (v->size > 1 && V2)
416         {
417             /* we are implementing the following:
418             if (V2*q_dash > (((U(0)*COMP_RADIX + U(1) - 
419                     q_dash*V1)*COMP_RADIX) + U(2))) ... */
420             comp inner = (comp)((long_comp)COMP_RADIX*U(0) + U(1) - 
421                                         (long_comp)q_dash*V1);
422             if ((long_comp)V2*q_dash > (long_comp)inner*COMP_RADIX + U(2))
423             {
424                 q_dash--;
425             }
426         }
427
428         /* multiply and subtract */
429         if (q_dash)
430         {
431             int is_negative;
432             tmp_u = bi_subtract(ctx, tmp_u, 
433                     bi_int_multiply(ctx, bi_copy(v), q_dash), &is_negative);
434             more_comps(tmp_u, n+1);
435
436             Q(j) = q_dash; 
437
438             /* add back */
439             if (is_negative)
440             {
441                 Q(j)--;
442                 tmp_u = bi_add(ctx, tmp_u, bi_copy(v));
443
444                 /* lop off the carry */
445                 tmp_u->size--;
446                 v->size--;
447             }
448         }
449         else
450         {
451             Q(j) = 0; 
452         }
453
454         /* copy back to u */
455         memcpy(&u->comps[u->size-n-1-j], tmp_u->comps, (n+1)*COMP_BYTE_SIZE);
456     } while (++j <= m);
457
458     bi_free(ctx, tmp_u);
459     bi_free(ctx, v);
460
461     if (is_mod)     /* get the remainder */
462     {
463         bi_free(ctx, quotient);
464         return bi_int_divide(ctx, trim(u), d);
465     }
466     else            /* get the quotient */
467     {
468         bi_free(ctx, u);
469         return trim(quotient);
470     }
471 }
472
473 /*
474  * Perform an integer divide on a bigint.
475  */
476 static bigint *bi_int_divide(BI_CTX *ctx, bigint *biR, comp denom)
477 {
478     int i = biR->size - 1;
479     long_comp r = 0;
480
481     check(biR);
482
483     do
484     {
485         r = (r<<COMP_BIT_SIZE) + biR->comps[i];
486         biR->comps[i] = (comp)(r / denom);
487         r %= denom;
488     } while (--i != 0);
489
490     return trim(biR);
491 }
492
493 #ifdef CONFIG_BIGINT_MONTGOMERY
494 /**
495  * There is a need for the value of integer N' such that B^-1(B-1)-N^-1N'=1, 
496  * where B^-1(B-1) mod N=1. Actually, only the least significant part of 
497  * N' is needed, hence the definition N0'=N' mod b. We reproduce below the 
498  * simple algorithm from an article by Dusse and Kaliski to efficiently 
499  * find N0' from N0 and b */
500 static comp modular_inverse(bigint *bim)
501 {
502     int i;
503     comp t = 1;
504     comp two_2_i_minus_1 = 2;   /* 2^(i-1) */
505     long_comp two_2_i = 4;      /* 2^i */
506     comp N = bim->comps[0];
507
508     for (i = 2; i <= COMP_BIT_SIZE; i++)
509     {
510         if ((long_comp)N*t%two_2_i >= two_2_i_minus_1)
511         {
512             t += two_2_i_minus_1;
513         }
514
515         two_2_i_minus_1 <<= 1;
516         two_2_i <<= 1;
517     }
518
519     return (comp)(COMP_RADIX-t);
520 }
521 #endif
522
523 #if defined(CONFIG_BIGINT_KARATSUBA) || defined(CONFIG_BIGINT_BARRETT) || \
524     defined(CONFIG_BIGINT_MONTGOMERY)
525 /**
526  * Take each component and shift down (in terms of components) 
527  */
528 static bigint *comp_right_shift(bigint *biR, int num_shifts)
529 {
530     int i = biR->size-num_shifts;
531     comp *x = biR->comps;
532     comp *y = &biR->comps[num_shifts];
533
534     check(biR);
535
536     if (i <= 0)     /* have we completely right shifted? */
537     {
538         biR->comps[0] = 0;  /* return 0 */
539         biR->size = 1;
540         return biR;
541     }
542
543     do
544     {
545         *x++ = *y++;
546     } while (--i > 0);
547
548     biR->size -= num_shifts;
549     return biR;
550 }
551
552 /**
553  * Take each component and shift it up (in terms of components) 
554  */
555 static bigint *comp_left_shift(bigint *biR, int num_shifts)
556 {
557     int i = biR->size-1;
558     comp *x, *y;
559
560     check(biR);
561
562     if (num_shifts <= 0)
563     {
564         return biR;
565     }
566
567     more_comps(biR, biR->size + num_shifts);
568
569     x = &biR->comps[i+num_shifts];
570     y = &biR->comps[i];
571
572     do
573     {
574         *x-- = *y--;
575     } while (i--);
576
577     memset(biR->comps, 0, num_shifts*COMP_BYTE_SIZE); /* zero LS comps */
578     return biR;
579 }
580 #endif
581
582 /**
583  * @brief Allow a binary sequence to be imported as a bigint.
584  * @param ctx [in]  The bigint session context.
585  * @param data [in] The data to be converted.
586  * @param size [in] The number of bytes of data.
587  * @return A bigint representing this data.
588  */
589 bigint *bi_import(BI_CTX *ctx, const uint8_t *data, int size)
590 {
591     bigint *biR = alloc(ctx, (size+COMP_BYTE_SIZE-1)/COMP_BYTE_SIZE);
592     int i, j = 0, offset = 0;
593
594     memset(biR->comps, 0, biR->size*COMP_BYTE_SIZE);
595
596     for (i = size-1; i >= 0; i--)
597     {
598         biR->comps[offset] += data[i] << (j*8);
599
600         if (++j == COMP_BYTE_SIZE)
601         {
602             j = 0;
603             offset ++;
604         }
605     }
606
607     return trim(biR);
608 }
609
610 #ifdef CONFIG_SSL_FULL_MODE
611 /**
612  * @brief The testharness uses this code to import text hex-streams and 
613  * convert them into bigints.
614  * @param ctx [in]  The bigint session context.
615  * @param data [in] A string consisting of hex characters. The characters must
616  * be in upper case.
617  * @return A bigint representing this data.
618  */
619 bigint *bi_str_import(BI_CTX *ctx, const char *data)
620 {
621     int size = strlen(data);
622     bigint *biR = alloc(ctx, (size+COMP_NUM_NIBBLES-1)/COMP_NUM_NIBBLES);
623     int i, j = 0, offset = 0;
624     memset(biR->comps, 0, biR->size*COMP_BYTE_SIZE);
625
626     for (i = size-1; i >= 0; i--)
627     {
628         int num = (data[i] <= '9') ? (data[i] - '0') : (data[i] - 'A' + 10);
629         biR->comps[offset] += num << (j*4);
630
631         if (++j == COMP_NUM_NIBBLES)
632         {
633             j = 0;
634             offset ++;
635         }
636     }
637
638     return biR;
639 }
640
641 void bi_print(const char *label, bigint *x)
642 {
643     int i, j;
644
645     if (x == NULL)
646     {
647         printf("%s: (null)\n", label);
648         return;
649     }
650
651     printf("%s: (size %d)\n", label, x->size);
652     for (i = x->size-1; i >= 0; i--)
653     {
654         for (j = COMP_NUM_NIBBLES-1; j >= 0; j--)
655         {
656             comp mask = 0x0f << (j*4);
657             comp num = (x->comps[i] & mask) >> (j*4);
658             putc((num <= 9) ? (num + '0') : (num + 'A' - 10), stdout);
659         }
660     }  
661
662     printf("\n");
663 }
664 #endif
665
666 /**
667  * @brief Take a bigint and convert it into a byte sequence. 
668  *
669  * This is useful after a decrypt operation.
670  * @param ctx [in]  The bigint session context.
671  * @param x [in]  The bigint to be converted.
672  * @param data [out] The converted data as a byte stream.
673  * @param size [in] The maximum size of the byte stream. Unused bytes will be
674  * zeroed.
675  */
676 void bi_export(BI_CTX *ctx, bigint *x, uint8_t *data, int size)
677 {
678     int i, j, k = size-1;
679
680     check(x);
681     memset(data, 0, size);  /* ensure all leading 0's are cleared */
682
683     for (i = 0; i < x->size; i++)
684     {
685         for (j = 0; j < COMP_BYTE_SIZE; j++)
686         {
687             comp mask = 0xff << (j*8);
688             int num = (x->comps[i] & mask) >> (j*8);
689             data[k--] = num;
690
691             if (k < 0)
692             {
693                 break;
694             }
695         }
696     }
697
698     bi_free(ctx, x);
699 }
700
701 /**
702  * @brief Pre-calculate some of the expensive steps in reduction. 
703  *
704  * This function should only be called once (normally when a session starts).
705  * When the session is over, bi_free_mod() should be called. bi_mod_power()
706  * relies on this function being called.
707  * @param ctx [in]  The bigint session context.
708  * @param bim [in]  The bigint modulus that will be used.
709  * @param mod_offset [in] There are three moduluii that can be stored - the
710  * standard modulus, and its two primes p and q. This offset refers to which
711  * modulus we are referring to.
712  * @see bi_free_mod(), bi_mod_power().
713  */
714 void bi_set_mod(BI_CTX *ctx, bigint *bim, int mod_offset)
715 {
716     int k = bim->size;
717     comp d = (comp)((long_comp)COMP_RADIX/(bim->comps[k-1]+1));
718 #ifdef CONFIG_BIGINT_MONTGOMERY
719     bigint *R, *R2;
720 #endif
721
722     ctx->bi_mod[mod_offset] = bim;
723     bi_permanent(ctx->bi_mod[mod_offset]);
724     ctx->bi_normalised_mod[mod_offset] = bi_int_multiply(ctx, bim, d);
725     bi_permanent(ctx->bi_normalised_mod[mod_offset]);
726
727 #if defined(CONFIG_BIGINT_MONTGOMERY)
728     /* set montgomery variables */
729     R = comp_left_shift(bi_clone(ctx, ctx->bi_radix), k-1);     /* R */
730     R2 = comp_left_shift(bi_clone(ctx, ctx->bi_radix), k*2-1);  /* R^2 */
731     ctx->bi_RR_mod_m[mod_offset] = bi_mod(ctx, R2);             /* R^2 mod m */
732     ctx->bi_R_mod_m[mod_offset] = bi_mod(ctx, R);               /* R mod m */
733
734     bi_permanent(ctx->bi_RR_mod_m[mod_offset]);
735     bi_permanent(ctx->bi_R_mod_m[mod_offset]);
736
737     ctx->N0_dash[mod_offset] = modular_inverse(ctx->bi_mod[mod_offset]);
738
739 #elif defined (CONFIG_BIGINT_BARRETT)
740     ctx->bi_mu[mod_offset] = 
741         bi_divide(ctx, comp_left_shift(
742             bi_clone(ctx, ctx->bi_radix), k*2-1), ctx->bi_mod[mod_offset], 0);
743     bi_permanent(ctx->bi_mu[mod_offset]);
744 #endif
745 }
746
747 /**
748  * @brief Used when cleaning various bigints at the end of a session.
749  * @param ctx [in]  The bigint session context.
750  * @param mod_offset [in] The offset to use.
751  * @see bi_set_mod().
752  */
753 void bi_free_mod(BI_CTX *ctx, int mod_offset)
754 {
755     bi_depermanent(ctx->bi_mod[mod_offset]);
756     bi_free(ctx, ctx->bi_mod[mod_offset]);
757 #if defined (CONFIG_BIGINT_MONTGOMERY)
758     bi_depermanent(ctx->bi_RR_mod_m[mod_offset]);
759     bi_depermanent(ctx->bi_R_mod_m[mod_offset]);
760     bi_free(ctx, ctx->bi_RR_mod_m[mod_offset]);
761     bi_free(ctx, ctx->bi_R_mod_m[mod_offset]);
762 #elif defined(CONFIG_BIGINT_BARRETT)
763     bi_depermanent(ctx->bi_mu[mod_offset]); 
764     bi_free(ctx, ctx->bi_mu[mod_offset]);
765 #endif
766     bi_depermanent(ctx->bi_normalised_mod[mod_offset]); 
767     bi_free(ctx, ctx->bi_normalised_mod[mod_offset]);
768 }
769
770 /** 
771  * Perform a standard multiplication between two bigints.
772  */
773 static bigint *regular_multiply(BI_CTX *ctx, bigint *bia, bigint *bib)
774 {
775     int i, j, i_plus_j;
776     int n = bia->size; 
777     int t = bib->size;
778     bigint *biR = alloc(ctx, n + t);
779     comp *sr = biR->comps;
780     comp *sa = bia->comps;
781     comp *sb = bib->comps;
782
783     check(bia);
784     check(bib);
785
786     /* clear things to start with */
787     memset(biR->comps, 0, ((n+t)*COMP_BYTE_SIZE));
788     i = 0;
789
790     do 
791     {
792         comp carry = 0;
793         comp b = *sb++;
794         i_plus_j = i;
795         j = 0;
796
797         do
798         {
799             long_comp tmp = sr[i_plus_j] + (long_comp)sa[j]*b + carry;
800             sr[i_plus_j++] = (comp)tmp;              /* downsize */
801             carry = (comp)(tmp >> COMP_BIT_SIZE);
802         } while (++j < n);
803
804         sr[i_plus_j] = carry;
805     } while (++i < t);
806
807     bi_free(ctx, bia);
808     bi_free(ctx, bib);
809     return trim(biR);
810 }
811
812 #ifdef CONFIG_BIGINT_KARATSUBA
813 /*
814  * Karatsuba improves on regular multiplication due to only 3 multiplications 
815  * being done instead of 4. The additional additions/subtractions are O(N) 
816  * rather than O(N^2) and so for big numbers it saves on a few operations 
817  */
818 static bigint *karatsuba(BI_CTX *ctx, bigint *bia, bigint *bib, int is_square)
819 {
820     bigint *x0, *x1;
821     bigint *p0, *p1, *p2;
822     int m;
823
824     if (is_square)
825     {
826         m = (bia->size + 1)/2;
827     }
828     else
829     {
830         m = (max(bia->size, bib->size) + 1)/2;
831     }
832
833     x0 = bi_clone(ctx, bia);
834     x0->size = m;
835     x1 = bi_clone(ctx, bia);
836     comp_right_shift(x1, m);
837     bi_free(ctx, bia);
838
839     /* work out the 3 partial products */
840     if (is_square)
841     {
842         p0 = bi_square(ctx, bi_copy(x0));
843         p2 = bi_square(ctx, bi_copy(x1));
844         p1 = bi_square(ctx, bi_add(ctx, x0, x1));
845     }
846     else /* normal multiply */
847     {
848         bigint *y0, *y1;
849         y0 = bi_clone(ctx, bib);
850         y0->size = m;
851         y1 = bi_clone(ctx, bib);
852         comp_right_shift(y1, m);
853         bi_free(ctx, bib);
854
855         p0 = bi_multiply(ctx, bi_copy(x0), bi_copy(y0));
856         p2 = bi_multiply(ctx, bi_copy(x1), bi_copy(y1));
857         p1 = bi_multiply(ctx, bi_add(ctx, x0, x1), bi_add(ctx, y0, y1));
858     }
859
860     p1 = bi_subtract(ctx, 
861             bi_subtract(ctx, p1, bi_copy(p2), NULL), bi_copy(p0), NULL);
862
863     comp_left_shift(p1, m);
864     comp_left_shift(p2, 2*m);
865     return bi_add(ctx, p1, bi_add(ctx, p0, p2));
866 }
867 #endif
868
869 /**
870  * @brief Perform a multiplication operation between two bigints.
871  * @param ctx [in]  The bigint session context.
872  * @param bia [in]  A bigint.
873  * @param bib [in]  Another bigint.
874  * @return The result of the multiplication.
875  */
876 bigint *bi_multiply(BI_CTX *ctx, bigint *bia, bigint *bib)
877 {
878     check(bia);
879     check(bib);
880
881 #ifdef CONFIG_BIGINT_KARATSUBA
882     if (min(bia->size, bib->size) < MUL_KARATSUBA_THRESH)
883     {
884         return regular_multiply(ctx, bia, bib);
885     }
886
887     return karatsuba(ctx, bia, bib, 0);
888 #else
889     return regular_multiply(ctx, bia, bib);
890 #endif
891 }
892
893 #ifdef CONFIG_BIGINT_SQUARE
894 /*
895  * Perform the actual square operion. It takes into account overflow.
896  */
897 static bigint *regular_square(BI_CTX *ctx, bigint *bi)
898 {
899     int t = bi->size;
900     int i = 0, j;
901     bigint *biR = alloc(ctx, t*2);
902     comp *w = biR->comps;
903     comp *x = bi->comps;
904     comp carry;
905
906     memset(w, 0, biR->size*COMP_BYTE_SIZE);
907
908     do
909     {
910         long_comp tmp = w[2*i] + (long_comp)x[i]*x[i];
911         comp u = 0;
912         w[2*i] = (comp)tmp;
913         carry = (comp)(tmp >> COMP_BIT_SIZE);
914
915         for (j = i+1; j < t; j++)
916         {
917             long_comp xx = (long_comp)x[i]*x[j];
918             long_comp blob = (long_comp)w[i+j]+carry;
919
920             if (u)                  /* previous overflow */
921             {
922                 blob += COMP_RADIX;
923             }
924
925             u = 0;
926             if (xx & COMP_BIG_MSB)  /* check for overflow */
927             {
928                 u = 1;
929             }
930
931             tmp = 2*xx + blob;
932             w[i+j] = (comp)tmp;
933             carry = (comp)(tmp >> COMP_BIT_SIZE);
934         }
935
936         w[i+t] += carry;
937
938         if (u)
939         {
940             w[i+t+1] = 1;   /* add carry */
941         }
942     } while (++i < t);
943
944     bi_free(ctx, bi);
945     return trim(biR);
946 }
947
948 /**
949  * @brief Perform a square operation on a bigint.
950  * @param ctx [in]  The bigint session context.
951  * @param bia [in]  A bigint.
952  * @return The result of the multiplication.
953  */
954 bigint *bi_square(BI_CTX *ctx, bigint *bia)
955 {
956     check(bia);
957
958 #ifdef CONFIG_BIGINT_KARATSUBA
959     if (bia->size < SQU_KARATSUBA_THRESH) 
960     {
961         return regular_square(ctx, bia);
962     }
963
964     return karatsuba(ctx, bia, NULL, 1);
965 #else
966     return regular_square(ctx, bia);
967 #endif
968 }
969 #endif
970
971 /**
972  * @brief Compare two bigints.
973  * @param bia [in]  A bigint.
974  * @param bib [in]  Another bigint.
975  * @return -1 if smaller, 1 if larger and 0 if equal.
976  */
977 int bi_compare(bigint *bia, bigint *bib)
978 {
979     int r, i;
980
981     check(bia);
982     check(bib);
983
984     if (bia->size > bib->size)
985         r = 1;
986     else if (bia->size < bib->size)
987         r = -1;
988     else
989     {
990         comp *a = bia->comps; 
991         comp *b = bib->comps; 
992
993         /* Same number of components.  Compare starting from the high end
994          * and working down. */
995         r = 0;
996         i = bia->size - 1;
997
998         do 
999         {
1000             if (a[i] > b[i])
1001             { 
1002                 r = 1;
1003                 break; 
1004             }
1005             else if (a[i] < b[i])
1006             { 
1007                 r = -1;
1008                 break; 
1009             }
1010         } while (--i >= 0);
1011     }
1012
1013     return r;
1014 }
1015
1016 /*
1017  * Allocate and zero more components.  Does not consume bi. 
1018  */
1019 static void more_comps(bigint *bi, int n)
1020 {
1021     if (n > bi->max_comps)
1022     {
1023         bi->max_comps = max(bi->max_comps * 2, n);
1024         bi->comps = (comp*)realloc(bi->comps, bi->max_comps * COMP_BYTE_SIZE);
1025     }
1026
1027     if (n > bi->size)
1028     {
1029         memset(&bi->comps[bi->size], 0, (n-bi->size)*COMP_BYTE_SIZE);
1030     }
1031
1032     bi->size = n;
1033 }
1034
1035 /*
1036  * Make a new empty bigint. It may just use an old one if one is available.
1037  * Otherwise get one off the heap.
1038  */
1039 static bigint *alloc(BI_CTX *ctx, int size)
1040 {
1041     bigint *biR;
1042
1043     /* Can we recycle an old bigint? */
1044     if (ctx->free_list != NULL)
1045     {
1046         biR = ctx->free_list;
1047         ctx->free_list = biR->next;
1048         ctx->free_count--;
1049
1050         if (biR->refs != 0)
1051         {
1052 #ifdef CONFIG_SSL_FULL_MODE
1053             printf("alloc: refs was not 0\n");
1054 #endif
1055             abort();    /* create a stack trace from a core dump */
1056         }
1057
1058         more_comps(biR, size);
1059     }
1060     else
1061     {
1062         /* No free bigints available - create a new one. */
1063         biR = (bigint *)malloc(sizeof(bigint));
1064         biR->comps = (comp*)malloc(size * COMP_BYTE_SIZE);
1065         biR->max_comps = size;  /* give some space to spare */
1066     }
1067
1068     biR->size = size;
1069     biR->refs = 1;
1070     biR->next = NULL;
1071     ctx->active_count++;
1072     return biR;
1073 }
1074
1075 /*
1076  * Work out the highest '1' bit in an exponent. Used when doing sliding-window
1077  * exponentiation.
1078  */
1079 static int find_max_exp_index(bigint *biexp)
1080 {
1081     int i = COMP_BIT_SIZE-1;
1082     comp shift = COMP_RADIX/2;
1083     comp test = biexp->comps[biexp->size-1];    /* assume no leading zeroes */
1084
1085     check(biexp);
1086
1087     do
1088     {
1089         if (test & shift)
1090         {
1091             return i+(biexp->size-1)*COMP_BIT_SIZE;
1092         }
1093
1094         shift >>= 1;
1095     } while (--i != 0);
1096
1097     return -1;      /* error - must have been a leading 0 */
1098 }
1099
1100 /*
1101  * Is a particular bit is an exponent 1 or 0? Used when doing sliding-window
1102  * exponentiation.
1103  */
1104 static int exp_bit_is_one(bigint *biexp, int offset)
1105 {
1106     comp test = biexp->comps[offset / COMP_BIT_SIZE];
1107     int num_shifts = offset % COMP_BIT_SIZE;
1108     comp shift = 1;
1109     int i;
1110
1111     check(biexp);
1112
1113     for (i = 0; i < num_shifts; i++)
1114     {
1115         shift <<= 1;
1116     }
1117
1118     return test & shift;
1119 }
1120
1121 #ifdef CONFIG_BIGINT_CHECK_ON
1122 /*
1123  * Perform a sanity check on bi.
1124  */
1125 static void check(const bigint *bi)
1126 {
1127     if (bi->refs <= 0)
1128     {
1129         printf("check: zero or negative refs in bigint\n");
1130         abort();
1131     }
1132
1133     if (bi->next != NULL)
1134     {
1135         printf("check: attempt to use a bigint from "
1136                 "the free list\n");
1137         abort();
1138     }
1139 }
1140 #endif
1141
1142 /*
1143  * Delete any leading 0's (and allow for 0).
1144  */
1145 static bigint *trim(bigint *bi)
1146 {
1147     check(bi);
1148
1149     while (bi->comps[bi->size-1] == 0 && bi->size > 1)
1150     {
1151         bi->size--;
1152     }
1153
1154     return bi;
1155 }
1156
1157 #if defined(CONFIG_BIGINT_MONTGOMERY)
1158 /**
1159  * @brief Perform a single montgomery reduction.
1160  * @param ctx [in]  The bigint session context.
1161  * @param bixy [in]  A bigint.
1162  * @return The result of the montgomery reduction.
1163  */
1164 bigint *bi_mont(BI_CTX *ctx, bigint *bixy)
1165 {
1166     int i = 0, n;
1167     uint8_t mod_offset = ctx->mod_offset;
1168     bigint *bim = ctx->bi_mod[mod_offset];
1169     comp mod_inv = ctx->N0_dash[mod_offset];
1170
1171     check(bixy);
1172
1173     if (ctx->use_classical)     /* just use classical instead */
1174     {
1175         return bi_mod(ctx, bixy);
1176     }
1177
1178     n = bim->size;
1179
1180     do
1181     {
1182         bixy = bi_add(ctx, bixy, comp_left_shift(
1183                     bi_int_multiply(ctx, bim, bixy->comps[i]*mod_inv), i));
1184     } while (++i < n);
1185
1186     comp_right_shift(bixy, n);
1187
1188     if (bi_compare(bixy, bim) >= 0)
1189     {
1190         bixy = bi_subtract(ctx, bixy, bim, NULL);
1191     }
1192
1193     return bixy;
1194 }
1195
1196 #elif defined(CONFIG_BIGINT_BARRETT)
1197 /*
1198  * Stomp on the most significant components to give the illusion of a "mod base
1199  * radix" operation 
1200  */
1201 static bigint *comp_mod(bigint *bi, int mod)
1202 {
1203     check(bi);
1204
1205     if (bi->size > mod)
1206     {
1207         bi->size = mod;
1208     }
1209
1210     return bi;
1211 }
1212
1213 /*
1214  * Barrett reduction has no need for some parts of the product, so ignore bits
1215  * of the multiply. This routine gives Barrett its big performance
1216  * improvements over Classical/Montgomery reduction methods. 
1217  */
1218 static bigint *partial_multiply(BI_CTX *ctx, bigint *bia, bigint *bib, 
1219         int inner_partial, int outer_partial)
1220 {
1221     int i = 0, j, n = bia->size, t = bib->size;
1222     bigint *biR;
1223     comp carry;
1224     comp *sr, *sa, *sb;
1225
1226     check(bia);
1227     check(bib);
1228
1229     biR = alloc(ctx, n + t);
1230     sa = bia->comps;
1231     sb = bib->comps;
1232     sr = biR->comps;
1233
1234     if (inner_partial)
1235     {
1236         memset(sr, 0, inner_partial*COMP_BYTE_SIZE); 
1237     }
1238     else    /* outer partial */
1239     {
1240         if (n < outer_partial || t < outer_partial) /* should we bother? */
1241         {
1242             bi_free(ctx, bia);
1243             bi_free(ctx, bib);
1244             biR->comps[0] = 0;      /* return 0 */
1245             biR->size = 1;
1246             return biR;
1247         }
1248
1249         memset(&sr[outer_partial], 0, (n+t-outer_partial)*COMP_BYTE_SIZE);
1250     }
1251
1252     do 
1253     {
1254         comp *a = sa;
1255         comp b = *sb++;
1256         long_comp tmp;
1257         int i_plus_j = i;
1258         carry = 0;
1259         j = n;
1260
1261         if (outer_partial && i_plus_j < outer_partial)
1262         {
1263             i_plus_j = outer_partial;
1264             a = &sa[outer_partial-i];
1265             j = n-(outer_partial-i);
1266         }
1267
1268         do
1269         {
1270             if (inner_partial && i_plus_j >= inner_partial) 
1271             {
1272                 break;
1273             }
1274
1275             tmp = sr[i_plus_j] + ((long_comp)*a++)*b + carry;
1276             sr[i_plus_j++] = (comp)tmp;              /* downsize */
1277             carry = (comp)(tmp >> COMP_BIT_SIZE);
1278         } while (--j != 0);
1279
1280         sr[i_plus_j] = carry;
1281     } while (++i < t);
1282
1283     bi_free(ctx, bia);
1284     bi_free(ctx, bib);
1285     return trim(biR);
1286 }
1287
1288 /**
1289  * @brief Perform a single Barrett reduction.
1290  * @param ctx [in]  The bigint session context.
1291  * @param bi [in]  A bigint.
1292  * @return The result of the Barrett reduction.
1293  */
1294 bigint *bi_barrett(BI_CTX *ctx, bigint *bi)
1295 {
1296     bigint *q1, *q2, *q3, *r1, *r2, *r;
1297     uint8_t mod_offset = ctx->mod_offset;
1298     bigint *bim = ctx->bi_mod[mod_offset];
1299     int k = bim->size;
1300
1301     check(bi);
1302     check(bim);
1303
1304     /* use Classical method instead  - Barrett cannot help here */
1305     if (bi->size > k*2)
1306     {
1307         return bi_mod(ctx, bi);
1308     }
1309
1310     q1 = comp_right_shift(bi_clone(ctx, bi), k-1);
1311
1312     /* do outer partial multiply */
1313     q2 = partial_multiply(ctx, q1, ctx->bi_mu[mod_offset], 0, k-1); 
1314     q3 = comp_right_shift(q2, k+1);
1315     r1 = comp_mod(bi, k+1);
1316
1317     /* do inner partial multiply */
1318     r2 = comp_mod(partial_multiply(ctx, q3, bim, k+1, 0), k+1);
1319     r = bi_subtract(ctx, r1, r2, NULL);
1320
1321     /* if (r >= m) r = r - m; */
1322     if (bi_compare(r, bim) >= 0)
1323     {
1324         r = bi_subtract(ctx, r, bim, NULL);
1325     }
1326
1327     return r;
1328 }
1329 #endif /* CONFIG_BIGINT_BARRETT */
1330
1331 #ifdef CONFIG_BIGINT_SLIDING_WINDOW
1332 /*
1333  * Work out g1, g3, g5, g7... etc for the sliding-window algorithm 
1334  */
1335 static void precompute_slide_window(BI_CTX *ctx, int window, bigint *g1)
1336 {
1337     int k = 1, i;
1338     bigint *g2;
1339
1340     for (i = 0; i < window-1; i++)   /* compute 2^(window-1) */
1341     {
1342         k <<= 1;
1343     }
1344
1345     ctx->g = (bigint **)malloc(k*sizeof(bigint *));
1346     ctx->g[0] = bi_clone(ctx, g1);
1347     bi_permanent(ctx->g[0]);
1348     g2 = bi_residue(ctx, bi_square(ctx, ctx->g[0]));   /* g^2 */
1349
1350     for (i = 1; i < k; i++)
1351     {
1352         ctx->g[i] = bi_residue(ctx, bi_multiply(ctx, ctx->g[i-1], bi_copy(g2)));
1353         bi_permanent(ctx->g[i]);
1354     }
1355
1356     bi_free(ctx, g2);
1357     ctx->window = k;
1358 }
1359 #endif
1360
1361 /**
1362  * @brief Perform a modular exponentiation.
1363  *
1364  * This function requires bi_set_mod() to have been called previously. This is 
1365  * one of the optimisations used for performance.
1366  * @param ctx [in]  The bigint session context.
1367  * @param bi  [in]  The bigint on which to perform the mod power operation.
1368  * @param biexp [in] The bigint exponent.
1369  * @see bi_set_mod().
1370  */
1371 bigint *bi_mod_power(BI_CTX *ctx, bigint *bi, bigint *biexp)
1372 {
1373     int i = find_max_exp_index(biexp), j, window_size = 1;
1374     bigint *biR = int_to_bi(ctx, 1);
1375
1376 #if defined(CONFIG_BIGINT_MONTGOMERY)
1377     uint8_t mod_offset = ctx->mod_offset;
1378     if (!ctx->use_classical)
1379     {
1380         /* preconvert */
1381         bi = bi_mont(ctx, 
1382                 bi_multiply(ctx, bi, ctx->bi_RR_mod_m[mod_offset]));    /* x' */
1383         bi_free(ctx, biR);
1384         biR = ctx->bi_R_mod_m[mod_offset];                              /* A */
1385     }
1386 #endif
1387
1388     check(bi);
1389     check(biexp);
1390
1391 #ifdef CONFIG_BIGINT_SLIDING_WINDOW
1392     for (j = i; j > 32; j /= 5) /* work out an optimum size */
1393         window_size++;
1394
1395     /* work out the slide constants */
1396     precompute_slide_window(ctx, window_size, bi);
1397 #else   /* just one constant */
1398     ctx->g = (bigint **)malloc(sizeof(bigint *));
1399     ctx->g[0] = bi_clone(ctx, bi);
1400     ctx->window = 1;
1401     bi_permanent(ctx->g[0]);
1402 #endif
1403
1404     /* if sliding-window is off, then only one bit will be done at a time and
1405      * will reduce to standard left-to-right exponentiation */
1406     do
1407     {
1408         if (exp_bit_is_one(biexp, i))
1409         {
1410             int l = i-window_size+1;
1411             int part_exp = 0;
1412
1413             if (l < 0)  /* LSB of exponent will always be 1 */
1414                 l = 0;
1415             else
1416             {
1417                 while (exp_bit_is_one(biexp, l) == 0)
1418                     l++;    /* go back up */
1419             }
1420
1421             /* build up the section of the exponent */
1422             for (j = i; j >= l; j--)
1423             {
1424                 biR = bi_residue(ctx, bi_square(ctx, biR));
1425                 if (exp_bit_is_one(biexp, j))
1426                     part_exp++;
1427
1428                 if (j != l)
1429                     part_exp <<= 1;
1430             }
1431
1432             part_exp = (part_exp-1)/2;  /* adjust for array */
1433             biR = bi_residue(ctx, bi_multiply(ctx, biR, ctx->g[part_exp]));
1434             i = l-1;
1435         }
1436         else    /* square it */
1437         {
1438             biR = bi_residue(ctx, bi_square(ctx, biR));
1439             i--;
1440         }
1441     } while (i >= 0);
1442      
1443     /* cleanup */
1444     for (i = 0; i < ctx->window; i++)
1445     {
1446         bi_depermanent(ctx->g[i]);
1447         bi_free(ctx, ctx->g[i]);
1448     }
1449
1450     free(ctx->g);
1451     bi_free(ctx, bi);
1452     bi_free(ctx, biexp);
1453 #if defined CONFIG_BIGINT_MONTGOMERY
1454     return ctx->use_classical ? biR : bi_mont(ctx, biR); /* convert back */
1455 #else /* CONFIG_BIGINT_CLASSICAL or CONFIG_BIGINT_BARRETT */
1456     return biR;
1457 #endif
1458 }
1459
1460 #ifdef CONFIG_SSL_CERT_VERIFICATION
1461 /**
1462  * @brief Perform a modular exponentiation using a temporary modulus.
1463  *
1464  * We need this function to check the signatures of certificates. The modulus
1465  * of this function is temporary as it's just used for authentication.
1466  * @param ctx [in]  The bigint session context.
1467  * @param bi  [in]  The bigint to perform the exp/mod.
1468  * @param bim [in]  The temporary modulus.
1469  * @param biexp [in] The bigint exponent.
1470  * @see bi_set_mod().
1471  */
1472 bigint *bi_mod_power2(BI_CTX *ctx, bigint *bi, bigint *bim, bigint *biexp)
1473 {
1474     bigint *biR, *tmp_biR;
1475
1476     /* Set up a temporary bigint context and transfer what we need between
1477      * them. We need to do this since we want to keep the original modulus
1478      * which is already in this context. This operation is only called when
1479      * doing peer verification, and so is not expensive :-) */
1480     BI_CTX *tmp_ctx = bi_initialize();
1481     bi_set_mod(tmp_ctx, bi_clone(tmp_ctx, bim), BIGINT_M_OFFSET);
1482     tmp_biR = bi_mod_power(tmp_ctx, 
1483                 bi_clone(tmp_ctx, bi), 
1484                 bi_clone(tmp_ctx, biexp));
1485     biR = bi_clone(ctx, tmp_biR);
1486     bi_free(tmp_ctx, tmp_biR);
1487     bi_free_mod(tmp_ctx, BIGINT_M_OFFSET);
1488     bi_terminate(tmp_ctx);
1489
1490     bi_free(ctx, bi);
1491     bi_free(ctx, bim);
1492     bi_free(ctx, biexp);
1493     return biR;
1494 }
1495 #endif
1496 /** @} */