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