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