2 #ifdef BN_MP_PRIME_NEXT_PRIME_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality.
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.
12 * The library is free for all purposes without any express
15 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
18 /* finds the next prime after the number "a" using "t" trials
21 * bbs_style = 1 means the prime must be congruent to 3 mod 4
23 int mp_prime_next_prime(mp_int
*a
, int t
, int bbs_style
)
26 mp_digit res_tab
[PRIME_SIZE
], step
, kstep
;
29 /* ensure t is valid */
30 if (t
<= 0 || t
> PRIME_SIZE
) {
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
) {
43 /* ok we found a prime smaller or
44 * equal [so the next is larger]
46 * however, the prime must be
47 * congruent to 3 mod 4
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
]);
59 mp_set(a
, ltm_prime_tab
[x
+ 1]);
64 /* at this point a maybe 1 */
65 if (mp_cmp_d(a
, 1) == MP_EQ
) {
69 /* fall through to the sieve */
72 /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
79 /* at this point we will use a combination of a sieve and Miller-Rabin */
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
; };
87 if (mp_iseven(a
) == 1) {
89 if ((err
= mp_sub_d(a
, 1, a
)) != MP_OKAY
) {
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
) {
102 /* init temp used for Miller-Rabin Testing */
103 if ((err
= mp_init(&b
)) != MP_OKAY
) {
108 /* skip to the next non-trivially divisible candidate */
111 /* y == 1 if any residue was zero [e.g. cannot be prime] */
114 /* increase step to next candidate */
117 /* compute the new residue without using division */
118 for (x
= 1; x
< PRIME_SIZE
; x
++) {
119 /* add the step to each residue */
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
];
127 /* set flag if zero */
128 if (res_tab
[x
] == 0) {
132 } while (y
== 1 && step
< ((((mp_digit
)1)<<DIGIT_BIT
) - kstep
));
135 if ((err
= mp_add_d(a
, step
, a
)) != MP_OKAY
) {
139 /* if didn't pass sieve and step == MAX then skip test */
140 if (y
== 1 && step
>= ((((mp_digit
)1)<<DIGIT_BIT
) - kstep
)) {
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
) {
168 /* $Source: /cvs/libtom/libtommath/bn_mp_prime_next_prime.c,v $ */
169 /* $Revision: 1.4 $ */
170 /* $Date: 2006/12/28 01:25:13 $ */