1 /* $NetBSD: bn_mp_prime_next_prime.c,v 1.2 2017/01/28 21:31:47 christos 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 */
mp_prime_next_prime(mp_int * a,int t,int bbs_style)25 int mp_prime_next_prime(mp_int *a, int t, int bbs_style)
26 {
27 int err, res = MP_NO, 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[x]);
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