1 /* TomsFastMath, a fast ISO C bignum library.
3 * This project is meant to fill in where LibTomMath
4 * falls short. That is speed ;-)
6 * This project is public domain and free for all purposes.
8 * Tom St Denis, tomstdenis@iahu.ca
19 /* Assure these -Pekka */
23 #define MIN(x,y) ((x)<(y)?(x):(y))
25 #define MAX(x,y) ((x)>(y)?(x):(y))
27 /* do we want large code? */
30 /* do we want huge code (implies large)? The answer is, yes. */
33 /* imply TFM_LARGE as required */
35 #if !defined(TFM_LARGE)
40 /* Max size of any number in bits. Basically the largest size you will be multiplying
41 * should be half [or smaller] of FP_MAX_SIZE-four_digit
43 * You can externally define this or it defaults to 4096-bits.
47 #define FP_MAX_SIZE (8192+(4*DIGIT_BIT))
48 /* #define FP_MAX_SIZE (4096+(4*DIGIT_BIT))*/
51 /* will this lib work? */
53 #error CHAR_BIT must be a multiple of eight.
55 #if FP_MAX_SIZE % CHAR_BIT
56 #error FP_MAX_SIZE must be a multiple of CHAR_BIT
59 /* autodetect x86-64 and make sure we are using 64-bit digits with x86-64 asm */
60 #if defined(__x86_64__)
61 #if defined(TFM_X86) || defined(TFM_SSE2) || defined(TFM_ARM)
62 #error x86-64 detected, x86-32/SSE2/ARM optimizations are not valid!
64 #if !defined(TFM_X86_64) && !defined(TFM_NO_ASM)
68 #if defined(TFM_X86_64)
69 #if !defined(FP_64BIT)
74 /* try to detect x86-32 */
75 #if defined(__i386__) && !defined(TFM_SSE2)
76 #if defined(TFM_X86_64) || defined(TFM_ARM)
77 #error x86-32 detected, x86-64/ARM optimizations are not valid!
79 #if !defined(TFM_X86) && !defined(TFM_NO_ASM)
84 /* make sure we're 32-bit for x86-32/sse/arm */
85 #if (defined(TFM_X86) || defined(TFM_SSE2) || defined(TFM_ARM)) && defined(FP_64BIT)
86 #warning x86-32, SSE2 and ARM optimizations require 32-bit digits (undefining)
96 #error TFM_ASM already defined!
102 #error TFM_ASM already defined!
108 #error TFM_ASM already defined!
113 /* we want no asm? */
122 /* some default configurations.
124 #if defined(FP_64BIT)
125 /* for GCC only on supported platforms */
127 typedef unsigned long ulong64;
129 typedef ulong64 fp_digit;
130 typedef unsigned long fp_word __attribute__ ((mode(TI)));
132 /* this is to make porting into LibTomCrypt easier :-) */
134 #if defined(_MSC_VER) || defined(__BORLANDC__)
135 typedef unsigned __int64 ulong64;
136 typedef signed __int64 long64;
138 typedef unsigned long long ulong64;
139 typedef signed long long long64;
142 typedef unsigned long fp_digit;
143 typedef ulong64 fp_word;
146 /* # of digits this is */
147 #define DIGIT_BIT (int)((CHAR_BIT) * sizeof(fp_digit))
148 #define FP_MASK (fp_digit)(-1)
149 #define FP_SIZE (FP_MAX_SIZE/DIGIT_BIT)
161 #define FP_LT -1 /* less than */
162 #define FP_EQ 0 /* equal to */
163 #define FP_GT 1 /* greater than */
166 #define FP_YES 1 /* yes response */
167 #define FP_NO 0 /* no response */
171 fp_digit dp[FP_SIZE];
178 /* returns a TFM ident string useful for debugging... */
179 const char *fp_ident(void);
181 /* initialize [or zero] an fp int */
182 #define fp_init(a) (void)memset((a), 0, sizeof(fp_int))
183 #define fp_zero(a) fp_init(a)
185 /* zero/even/odd ? */
186 #define fp_iszero(a) (((a)->used == 0) ? FP_YES : FP_NO)
187 #define fp_iseven(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 0)) ? FP_YES : FP_NO)
188 #define fp_isodd(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 1)) ? FP_YES : FP_NO)
190 /* set to a small digit */
191 void fp_set(fp_int *a, fp_digit b);
193 /* copy from a to b */
194 #define fp_copy(a, b) (void)(((a) != (b)) && memcpy((b), (a), sizeof(fp_int)))
195 #define fp_init_copy(a, b) fp_copy(b, a)
197 /* negate and absolute */
198 #define fp_neg(a, b) { fp_copy(a, b); (b)->sign ^= 1; }
199 #define fp_abs(a, b) { fp_copy(a, b); (b)->sign = 0; }
202 #define fp_clamp(a) { while ((a)->used && (a)->dp[(a)->used-1] == 0) --((a)->used); (a)->sign = (a)->used ? (a)->sign : FP_ZPOS; }
204 /* right shift x digits */
205 void fp_rshd(fp_int *a, int x);
207 /* left shift x digits */
208 void fp_lshd(fp_int *a, int x);
210 /* signed comparison */
211 int fp_cmp(fp_int *a, fp_int *b);
213 /* unsigned comparison */
214 int fp_cmp_mag(fp_int *a, fp_int *b);
216 /* power of 2 operations */
217 void fp_div_2d(fp_int *a, int b, fp_int *c, fp_int *d);
218 void fp_mod_2d(fp_int *a, int b, fp_int *c);
219 void fp_mul_2d(fp_int *a, int b, fp_int *c);
220 void fp_2expt (fp_int *a, int b);
221 void fp_mul_2(fp_int *a, fp_int *c);
222 void fp_div_2(fp_int *a, fp_int *c);
224 /* Counts the number of lsbs which are zero before the first zero bit */
225 int fp_cnt_lsb(fp_int *a);
228 void fp_add(fp_int *a, fp_int *b, fp_int *c);
231 void fp_sub(fp_int *a, fp_int *b, fp_int *c);
234 void fp_mul(fp_int *a, fp_int *b, fp_int *c);
237 void fp_sqr(fp_int *a, fp_int *b);
239 /* a/b => cb + d == a */
240 int fp_div(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
242 /* c = a mod b, 0 <= c < b */
243 int fp_mod(fp_int *a, fp_int *b, fp_int *c);
245 /* compare against a single digit */
246 int fp_cmp_d(fp_int *a, fp_digit b);
249 void fp_add_d(fp_int *a, fp_digit b, fp_int *c);
252 void fp_sub_d(fp_int *a, fp_digit b, fp_int *c);
255 void fp_mul_d(fp_int *a, fp_digit b, fp_int *c);
257 /* a/b => cb + d == a */
258 int fp_div_d(fp_int *a, fp_digit b, fp_int *c, fp_digit *d);
260 /* c = a mod b, 0 <= c < b */
261 int fp_mod_d(fp_int *a, fp_digit b, fp_digit *c);
263 /* ---> number theory <--- */
264 /* d = a + b (mod c) */
265 int fp_addmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
267 /* d = a - b (mod c) */
268 int fp_submod(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
270 /* d = a * b (mod c) */
271 int fp_mulmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
273 /* c = a * a (mod b) */
274 int fp_sqrmod(fp_int *a, fp_int *b, fp_int *c);
276 /* c = 1/a (mod b) */
277 int fp_invmod(fp_int *a, fp_int *b, fp_int *c);
280 void fp_gcd(fp_int *a, fp_int *b, fp_int *c);
283 void fp_lcm(fp_int *a, fp_int *b, fp_int *c);
285 /* setups the montgomery reduction */
286 int fp_montgomery_setup(fp_int *a, fp_digit *mp);
288 /* computes a = B**n mod b without division or multiplication useful for
289 * normalizing numbers in a Montgomery system.
291 void fp_montgomery_calc_normalization(fp_int *a, fp_int *b);
293 /* computes x/R == x (mod N) via Montgomery Reduction */
294 void fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp);
296 /* d = a**b (mod c) */
297 int fp_exptmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
299 /* primality stuff */
301 /* perform a Miller-Rabin test of a to the base b and store result in "result" */
302 void fp_prime_miller_rabin (fp_int * a, fp_int * b, int *result);
304 /* 256 trial divisions + 8 Miller-Rabins, returns FP_YES if probable prime */
305 int fp_isprime(fp_int *a);
307 /* Primality generation flags */
308 #define TFM_PRIME_BBS 0x0001 /* BBS style prime */
309 #define TFM_PRIME_SAFE 0x0002 /* Safe prime (p-1)/2 == prime */
310 #define TFM_PRIME_2MSB_OFF 0x0004 /* force 2nd MSB to 0 */
311 #define TFM_PRIME_2MSB_ON 0x0008 /* force 2nd MSB to 1 */
313 /* callback for fp_prime_random, should fill dst with random bytes and return how many read [upto len] */
314 typedef int tfm_prime_callback(unsigned char *dst, int len, void *dat);
316 #define fp_prime_random(a, t, size, bbs, cb, dat) fp_prime_random_ex(a, t, ((size) * 8) + 1, (bbs==1)?TFM_PRIME_BBS:0, cb, dat)
318 int fp_prime_random_ex(fp_int *a, int t, int size, int flags, tfm_prime_callback cb, void *dat);
320 /* radix conersions */
321 int fp_count_bits(fp_int *a);
323 int fp_unsigned_bin_size(fp_int *a);
324 void fp_read_unsigned_bin(fp_int *a, unsigned char *b, int c);
325 void fp_to_unsigned_bin(fp_int *a, unsigned char *b);
327 int fp_signed_bin_size(fp_int *a);
328 void fp_read_signed_bin(fp_int *a, unsigned char *b, int c);
329 void fp_to_signed_bin(fp_int *a, unsigned char *b);
331 int fp_read_radix(fp_int *a, char *str, int radix);
332 int fp_toradix(fp_int *a, char *str, int radix);
333 int fp_toradix_n(fp_int * a, char *str, int radix, int maxlen);
334 int fp_radix_size(fp_int *a, int radix, int *size);
336 /* VARIOUS LOW LEVEL STUFFS */
337 void s_fp_add(fp_int *a, fp_int *b, fp_int *c);
338 void s_fp_sub(fp_int *a, fp_int *b, fp_int *c);
339 void bn_reverse(unsigned char *s, int len);
340 void fp_mul_comba(fp_int *A, fp_int *B, fp_int *C);
342 void fp_mul_comba32(fp_int *A, fp_int *B, fp_int *C);
345 void fp_mul_comba16(fp_int *A, fp_int *B, fp_int *C);
347 void fp_mul_comba8(fp_int *A, fp_int *B, fp_int *C);
348 void fp_mul_comba4(fp_int *A, fp_int *B, fp_int *C);
350 void fp_sqr_comba(fp_int *A, fp_int *B);
351 void fp_sqr_comba4(fp_int *A, fp_int *B);
352 void fp_sqr_comba8(fp_int *A, fp_int *B);
354 void fp_sqr_comba16(fp_int *A, fp_int *B);
357 void fp_sqr_comba32(fp_int *A, fp_int *B);
358 void fp_sqr_comba64(fp_int *A, fp_int *B);
360 extern const char *fp_s_rmap;