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