tommath/bn_mp_prime_next_prime.c

171 lines
4.6 KiB
C
Raw Normal View History

2004-10-29 18:07:18 -04:00
#include <tommath.h>
#ifdef BN_MP_PRIME_NEXT_PRIME_C
2003-03-22 10:10:20 -05:00
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
2003-08-04 21:24:44 -04:00
* LibTomMath is a library that provides multiple-precision
2003-03-22 10:10:20 -05:00
* integer arithmetic as well as number theoretic functionality.
*
2003-08-04 21:24:44 -04:00
* The library was designed directly after the MPI library by
2003-03-22 10:10:20 -05:00
* 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.
*
2007-04-18 05:58:18 -04:00
* Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2003-03-22 10:10:20 -05:00
*/
/* finds the next prime after the number "a" using "t" trials
* of Miller-Rabin.
2003-07-12 10:31:43 -04:00
*
* bbs_style = 1 means the prime must be congruent to 3 mod 4
2003-03-22 10:10:20 -05:00
*/
2003-07-12 10:31:43 -04:00
int mp_prime_next_prime(mp_int *a, int t, int bbs_style)
2003-03-22 10:10:20 -05:00
{
int err, res = MP_NO, x, y;
2003-07-12 10:31:43 -04:00
mp_digit res_tab[PRIME_SIZE], step, kstep;
mp_int b;
2003-05-17 08:33:54 -04:00
2003-07-12 10:31:43 -04:00
/* ensure t is valid */
if (t <= 0 || t > PRIME_SIZE) {
return MP_VAL;
}
/* force positive */
2003-08-04 21:24:44 -04:00
a->sign = MP_ZPOS;
2003-07-12 10:31:43 -04:00
/* simple algo if a is less than the largest prime in the table */
2004-12-22 21:40:37 -05:00
if (mp_cmp_d(a, ltm_prime_tab[PRIME_SIZE-1]) == MP_LT) {
2003-07-12 10:31:43 -04:00
/* find which prime it is bigger than */
for (x = PRIME_SIZE - 2; x >= 0; x--) {
2004-12-22 21:40:37 -05:00
if (mp_cmp_d(a, ltm_prime_tab[x]) != MP_LT) {
2003-07-12 10:31:43 -04:00
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
*/
2004-12-22 21:40:37 -05:00
if ((ltm_prime_tab[x + 1] & 3) != 3) {
2003-07-12 10:31:43 -04:00
/* scan upwards for a prime congruent to 3 mod 4 */
for (y = x + 1; y < PRIME_SIZE; y++) {
2004-12-22 21:40:37 -05:00
if ((ltm_prime_tab[y] & 3) == 3) {
mp_set(a, ltm_prime_tab[y]);
2003-07-12 10:31:43 -04:00
return MP_OKAY;
}
}
}
} else {
2004-12-22 21:40:37 -05:00
mp_set(a, ltm_prime_tab[x + 1]);
2003-07-12 10:31:43 -04:00
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; };
2003-03-22 10:10:20 -05:00
}
} else {
2003-07-12 10:31:43 -04:00
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++) {
2004-12-22 21:40:37 -05:00
if ((err = mp_mod_d(a, ltm_prime_tab[x], res_tab + x)) != MP_OKAY) {
2003-03-22 10:10:20 -05:00
return err;
}
2003-05-17 08:33:54 -04:00
}
2003-07-12 10:31:43 -04:00
/* init temp used for Miller-Rabin Testing */
if ((err = mp_init(&b)) != MP_OKAY) {
return err;
}
2003-03-22 10:10:20 -05:00
for (;;) {
2003-07-12 10:31:43 -04:00
/* 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;
2003-07-15 20:26:58 -04:00
/* increase step to next candidate */
2003-07-12 10:31:43 -04:00
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] */
2004-12-22 21:40:37 -05:00
if (res_tab[x] >= ltm_prime_tab[x]) {
res_tab[x] -= ltm_prime_tab[x];
2003-07-12 10:31:43 -04:00
}
/* 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) {
2004-12-22 21:40:37 -05:00
goto LBL_ERR;
2003-07-12 10:31:43 -04:00
}
2003-08-04 21:24:44 -04:00
/* if didn't pass sieve and step == MAX then skip test */
if (y == 1 && step >= ((((mp_digit)1)<<DIGIT_BIT) - kstep)) {
2003-07-12 10:31:43 -04:00
continue;
}
2003-03-22 10:10:20 -05:00
/* is this prime? */
2003-07-12 10:31:43 -04:00
for (x = 0; x < t; x++) {
mp_set(&b, ltm_prime_tab[x]);
2003-07-12 10:31:43 -04:00
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_ERR;
2003-07-12 10:31:43 -04:00
}
2003-12-24 13:59:22 -05:00
if (res == MP_NO) {
2003-07-12 10:31:43 -04:00
break;
}
2003-03-22 10:10:20 -05:00
}
2003-05-17 08:33:54 -04:00
2003-12-24 13:59:22 -05:00
if (res == MP_YES) {
2003-03-22 10:10:20 -05:00
break;
}
}
2003-05-17 08:33:54 -04:00
2003-07-12 10:31:43 -04:00
err = MP_OKAY;
2004-12-22 21:40:37 -05:00
LBL_ERR:
2003-07-12 10:31:43 -04:00
mp_clear(&b);
return err;
2003-03-22 10:10:20 -05:00
}
2004-10-29 18:07:18 -04:00
#endif
2005-08-01 12:37:28 -04:00
/* $Source$ */
/* $Revision$ */
/* $Date$ */