1/*	$NetBSD$	*/
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 */
25int 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;
163LBL_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