Mercurial > dropbear

comparison libtommath/bn_mp_prime_next_prime.c @ 399:a707e6148060

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge of '5fdf69ca60d1683cdd9f4c2595134bed26394834' and '6b61c50f4cf888bea302ac8fcf5dbb573b443251'
author Matt Johnston <matt@ucc.asn.au>
date Sat, 03 Feb 2007 08:20:34 +0000
parents 5ff8218bcee9
children a55b97f5a485
comparison
equal deleted inserted replaced
394:17d097fc111c 399:a707e6148060
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, [email protected], http://math.libtomcrypt.com
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 */
23 int mp_prime_next_prime(mp_int *a, int t, int bbs_style)
24 {
25 int err, res, 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;
161 LBL_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.3 $ */
170 /* $Date: 2006/03/31 14:18:44 $ */