Mercurial > dropbear

comparison libtommath/bn_mp_prime_next_prime.c @ 284:eed26cff980b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
propagate from branch 'au.asn.ucc.matt.ltm.dropbear' (head 6c790cad5a7fa866ad062cb3a0c279f7ba788583) to branch 'au.asn.ucc.matt.dropbear' (head fff0894a0399405a9410ea1c6d118f342cf2aa64)
author Matt Johnston <matt@ucc.asn.au>
date Wed, 08 Mar 2006 13:23:49 +0000
parents
children 5ff8218bcee9
comparison
equal deleted inserted replaced
283:bd240aa12ba7 284:eed26cff980b
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.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 */
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