Mercurial > dropbear

comparison bn_mp_prime_next_prime.c @ 1:22d5cf7d4b1a libtommath

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