1 /* $NetBSD: bn_mp_prime_next_prime.c,v 1.1.1.1 2011/04/13 18:14:54 elric Exp $ */ 2 3 #include <tommath.h> 4 #ifdef BN_MP_PRIME_NEXT_PRIME_C 5 /* LibTomMath, multiple-precision integer library -- Tom St Denis 6 * 7 * LibTomMath is a library that provides multiple-precision 8 * integer arithmetic as well as number theoretic functionality. 9 * 10 * The library was designed directly after the MPI library by 11 * Michael Fromberger but has been written from scratch with 12 * additional optimizations in place. 13 * 14 * The library is free for all purposes without any express 15 * guarantee it works. 16 * 17 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org 18 */ 19 20 /* finds the next prime after the number "a" using "t" trials 21 * of Miller-Rabin. 22 * 23 * bbs_style = 1 means the prime must be congruent to 3 mod 4 24 */ 25 int mp_prime_next_prime(mp_int *a, int t, int bbs_style) 26 { 27 int err, res, x, y; 28 mp_digit res_tab[PRIME_SIZE], step, kstep; 29 mp_int b; 30 31 /* ensure t is valid */ 32 if (t <= 0 || t > PRIME_SIZE) { 33 return MP_VAL; 34 } 35 36 /* force positive */ 37 a->sign = MP_ZPOS; 38 39 /* simple algo if a is less than the largest prime in the table */ 40 if (mp_cmp_d(a, ltm_prime_tab[PRIME_SIZE-1]) == MP_LT) { 41 /* find which prime it is bigger than */ 42 for (x = PRIME_SIZE - 2; x >= 0; x--) { 43 if (mp_cmp_d(a, ltm_prime_tab[x]) != MP_LT) { 44 if (bbs_style == 1) { 45 /* ok we found a prime smaller or 46 * equal [so the next is larger] 47 * 48 * however, the prime must be 49 * congruent to 3 mod 4 50 */ 51 if ((ltm_prime_tab[x + 1] & 3) != 3) { 52 /* scan upwards for a prime congruent to 3 mod 4 */ 53 for (y = x + 1; y < PRIME_SIZE; y++) { 54 if ((ltm_prime_tab[y] & 3) == 3) { 55 mp_set(a, ltm_prime_tab[y]); 56 return MP_OKAY; 57 } 58 } 59 } 60 } else { 61 mp_set(a, ltm_prime_tab[x + 1]); 62 return MP_OKAY; 63 } 64 } 65 } 66 /* at this point a maybe 1 */ 67 if (mp_cmp_d(a, 1) == MP_EQ) { 68 mp_set(a, 2); 69 return MP_OKAY; 70 } 71 /* fall through to the sieve */ 72 } 73 74 /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */ 75 if (bbs_style == 1) { 76 kstep = 4; 77 } else { 78 kstep = 2; 79 } 80 81 /* at this point we will use a combination of a sieve and Miller-Rabin */ 82 83 if (bbs_style == 1) { 84 /* if a mod 4 != 3 subtract the correct value to make it so */ 85 if ((a->dp[0] & 3) != 3) { 86 if ((err = mp_sub_d(a, (a->dp[0] & 3) + 1, a)) != MP_OKAY) { return err; }; 87 } 88 } else { 89 if (mp_iseven(a) == 1) { 90 /* force odd */ 91 if ((err = mp_sub_d(a, 1, a)) != MP_OKAY) { 92 return err; 93 } 94 } 95 } 96 97 /* generate the restable */ 98 for (x = 1; x < PRIME_SIZE; x++) { 99 if ((err = mp_mod_d(a, ltm_prime_tab[x], res_tab + x)) != MP_OKAY) { 100 return err; 101 } 102 } 103 104 /* init temp used for Miller-Rabin Testing */ 105 if ((err = mp_init(&b)) != MP_OKAY) { 106 return err; 107 } 108 109 for (;;) { 110 /* skip to the next non-trivially divisible candidate */ 111 step = 0; 112 do { 113 /* y == 1 if any residue was zero [e.g. cannot be prime] */ 114 y = 0; 115 116 /* increase step to next candidate */ 117 step += kstep; 118 119 /* compute the new residue without using division */ 120 for (x = 1; x < PRIME_SIZE; x++) { 121 /* add the step to each residue */ 122 res_tab[x] += kstep; 123 124 /* subtract the modulus [instead of using division] */ 125 if (res_tab[x] >= ltm_prime_tab[x]) { 126 res_tab[x] -= ltm_prime_tab[x]; 127 } 128 129 /* set flag if zero */ 130 if (res_tab[x] == 0) { 131 y = 1; 132 } 133 } 134 } while (y == 1 && step < ((((mp_digit)1)<<DIGIT_BIT) - kstep)); 135 136 /* add the step */ 137 if ((err = mp_add_d(a, step, a)) != MP_OKAY) { 138 goto LBL_ERR; 139 } 140 141 /* if didn't pass sieve and step == MAX then skip test */ 142 if (y == 1 && step >= ((((mp_digit)1)<<DIGIT_BIT) - kstep)) { 143 continue; 144 } 145 146 /* is this prime? */ 147 for (x = 0; x < t; x++) { 148 mp_set(&b, ltm_prime_tab[t]); 149 if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) { 150 goto LBL_ERR; 151 } 152 if (res == MP_NO) { 153 break; 154 } 155 } 156 157 if (res == MP_YES) { 158 break; 159 } 160 } 161 162 err = MP_OKAY; 163 LBL_ERR: 164 mp_clear(&b); 165 return err; 166 } 167 168 #endif 169 170 /* Source: /cvs/libtom/libtommath/bn_mp_prime_next_prime.c,v */ 171 /* Revision: 1.4 */ 172 /* Date: 2006/12/28 01:25:13 */ 173