diff options
Diffstat (limited to 'libtommath/bn_mp_prime_next_prime.c')
| -rw-r--r-- | libtommath/bn_mp_prime_next_prime.c | 166 |
1 files changed, 166 insertions, 0 deletions
diff --git a/libtommath/bn_mp_prime_next_prime.c b/libtommath/bn_mp_prime_next_prime.c new file mode 100644 index 0000000..24f93c4 --- /dev/null +++ b/libtommath/bn_mp_prime_next_prime.c @@ -0,0 +1,166 @@ +#include <tommath.h> +#ifdef BN_MP_PRIME_NEXT_PRIME_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis + * + * LibTomMath is a library that provides multiple-precision + * integer arithmetic as well as number theoretic functionality. + * + * The library was designed directly after the MPI library by + * Michael Fromberger but has been written from scratch with + * additional optimizations in place. + * + * The library is free for all purposes without any express + * guarantee it works. + * + * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org + */ + +/* finds the next prime after the number "a" using "t" trials + * of Miller-Rabin. + * + * bbs_style = 1 means the prime must be congruent to 3 mod 4 + */ +int mp_prime_next_prime(mp_int *a, int t, int bbs_style) +{ + int err, res, x, y; + mp_digit res_tab[PRIME_SIZE], step, kstep; + mp_int b; + + /* ensure t is valid */ + if (t <= 0 || t > PRIME_SIZE) { + return MP_VAL; + } + + /* force positive */ + a->sign = MP_ZPOS; + + /* simple algo if a is less than the largest prime in the table */ + if (mp_cmp_d(a, ltm_prime_tab[PRIME_SIZE-1]) == MP_LT) { + /* find which prime it is bigger than */ + for (x = PRIME_SIZE - 2; x >= 0; x--) { + if (mp_cmp_d(a, ltm_prime_tab[x]) != MP_LT) { + if (bbs_style == 1) { + /* ok we found a prime smaller or + * equal [so the next is larger] + * + * however, the prime must be + * congruent to 3 mod 4 + */ + if ((ltm_prime_tab[x + 1] & 3) != 3) { + /* scan upwards for a prime congruent to 3 mod 4 */ + for (y = x + 1; y < PRIME_SIZE; y++) { + if ((ltm_prime_tab[y] & 3) == 3) { + mp_set(a, ltm_prime_tab[y]); + return MP_OKAY; + } + } + } + } else { + mp_set(a, ltm_prime_tab[x + 1]); + return MP_OKAY; + } + } + } + /* at this point a maybe 1 */ + if (mp_cmp_d(a, 1) == MP_EQ) { + mp_set(a, 2); + return MP_OKAY; + } + /* fall through to the sieve */ + } + + /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */ + if (bbs_style == 1) { + kstep = 4; + } else { + kstep = 2; + } + + /* at this point we will use a combination of a sieve and Miller-Rabin */ + + if (bbs_style == 1) { + /* if a mod 4 != 3 subtract the correct value to make it so */ + if ((a->dp[0] & 3) != 3) { + if ((err = mp_sub_d(a, (a->dp[0] & 3) + 1, a)) != MP_OKAY) { return err; }; + } + } else { + if (mp_iseven(a) == 1) { + /* force odd */ + if ((err = mp_sub_d(a, 1, a)) != MP_OKAY) { + return err; + } + } + } + + /* generate the restable */ + for (x = 1; x < PRIME_SIZE; x++) { + if ((err = mp_mod_d(a, ltm_prime_tab[x], res_tab + x)) != MP_OKAY) { + return err; + } + } + + /* init temp used for Miller-Rabin Testing */ + if ((err = mp_init(&b)) != MP_OKAY) { + return err; + } + + for (;;) { + /* skip to the next non-trivially divisible candidate */ + step = 0; + do { + /* y == 1 if any residue was zero [e.g. cannot be prime] */ + y = 0; + + /* increase step to next candidate */ + step += kstep; + + /* compute the new residue without using division */ + for (x = 1; x < PRIME_SIZE; x++) { + /* add the step to each residue */ + res_tab[x] += kstep; + + /* subtract the modulus [instead of using division] */ + if (res_tab[x] >= ltm_prime_tab[x]) { + res_tab[x] -= ltm_prime_tab[x]; + } + + /* set flag if zero */ + if (res_tab[x] == 0) { + y = 1; + } + } + } while (y == 1 && step < ((((mp_digit)1)<<DIGIT_BIT) - kstep)); + + /* add the step */ + if ((err = mp_add_d(a, step, a)) != MP_OKAY) { + goto LBL_ERR; + } + + /* if didn't pass sieve and step == MAX then skip test */ + if (y == 1 && step >= ((((mp_digit)1)<<DIGIT_BIT) - kstep)) { + continue; + } + + /* is this prime? */ + for (x = 0; x < t; x++) { + mp_set(&b, ltm_prime_tab[t]); + if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) { + goto LBL_ERR; + } + if (res == MP_NO) { + break; + } + } + + if (res == MP_YES) { + break; + } + } + + err = MP_OKAY; +LBL_ERR: + mp_clear(&b); + return err; +} + +#endif |