HEIMDAL: move code from source4/heimdal* to third_party/heimdal*
[samba.git] / third_party / heimdal / lib / hcrypto / libtommath / bn_mp_prime_next_prime.c
1 #include "tommath_private.h"
2 #ifdef BN_MP_PRIME_NEXT_PRIME_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
4 /* SPDX-License-Identifier: Unlicense */
5
6 /* finds the next prime after the number "a" using "t" trials
7  * of Miller-Rabin.
8  *
9  * bbs_style = 1 means the prime must be congruent to 3 mod 4
10  */
11 mp_err mp_prime_next_prime(mp_int *a, int t, int bbs_style)
12 {
13    int      x, y;
14    mp_ord   cmp;
15    mp_err   err;
16    mp_bool  res = MP_NO;
17    mp_digit res_tab[PRIVATE_MP_PRIME_TAB_SIZE], step, kstep;
18    mp_int   b;
19
20    /* force positive */
21    a->sign = MP_ZPOS;
22
23    /* simple algo if a is less than the largest prime in the table */
24    if (mp_cmp_d(a, s_mp_prime_tab[PRIVATE_MP_PRIME_TAB_SIZE-1]) == MP_LT) {
25       /* find which prime it is bigger than "a" */
26       for (x = 0; x < PRIVATE_MP_PRIME_TAB_SIZE; x++) {
27          cmp = mp_cmp_d(a, s_mp_prime_tab[x]);
28          if (cmp == MP_EQ) {
29             continue;
30          }
31          if (cmp != MP_GT) {
32             if ((bbs_style == 1) && ((s_mp_prime_tab[x] & 3u) != 3u)) {
33                /* try again until we get a prime congruent to 3 mod 4 */
34                continue;
35             } else {
36                mp_set(a, s_mp_prime_tab[x]);
37                return MP_OKAY;
38             }
39          }
40       }
41       /* fall through to the sieve */
42    }
43
44    /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
45    if (bbs_style == 1) {
46       kstep   = 4;
47    } else {
48       kstep   = 2;
49    }
50
51    /* at this point we will use a combination of a sieve and Miller-Rabin */
52
53    if (bbs_style == 1) {
54       /* if a mod 4 != 3 subtract the correct value to make it so */
55       if ((a->dp[0] & 3u) != 3u) {
56          if ((err = mp_sub_d(a, (a->dp[0] & 3u) + 1u, a)) != MP_OKAY) {
57             return err;
58          }
59       }
60    } else {
61       if (MP_IS_EVEN(a)) {
62          /* force odd */
63          if ((err = mp_sub_d(a, 1uL, a)) != MP_OKAY) {
64             return err;
65          }
66       }
67    }
68
69    /* generate the restable */
70    for (x = 1; x < PRIVATE_MP_PRIME_TAB_SIZE; x++) {
71       if ((err = mp_mod_d(a, s_mp_prime_tab[x], res_tab + x)) != MP_OKAY) {
72          return err;
73       }
74    }
75
76    /* init temp used for Miller-Rabin Testing */
77    if ((err = mp_init(&b)) != MP_OKAY) {
78       return err;
79    }
80
81    for (;;) {
82       /* skip to the next non-trivially divisible candidate */
83       step = 0;
84       do {
85          /* y == 1 if any residue was zero [e.g. cannot be prime] */
86          y     =  0;
87
88          /* increase step to next candidate */
89          step += kstep;
90
91          /* compute the new residue without using division */
92          for (x = 1; x < PRIVATE_MP_PRIME_TAB_SIZE; x++) {
93             /* add the step to each residue */
94             res_tab[x] += kstep;
95
96             /* subtract the modulus [instead of using division] */
97             if (res_tab[x] >= s_mp_prime_tab[x]) {
98                res_tab[x]  -= s_mp_prime_tab[x];
99             }
100
101             /* set flag if zero */
102             if (res_tab[x] == 0u) {
103                y = 1;
104             }
105          }
106       } while ((y == 1) && (step < (((mp_digit)1 << MP_DIGIT_BIT) - kstep)));
107
108       /* add the step */
109       if ((err = mp_add_d(a, step, a)) != MP_OKAY) {
110          goto LBL_ERR;
111       }
112
113       /* if didn't pass sieve and step == MP_MAX then skip test */
114       if ((y == 1) && (step >= (((mp_digit)1 << MP_DIGIT_BIT) - kstep))) {
115          continue;
116       }
117
118       if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
119          goto LBL_ERR;
120       }
121       if (res == MP_YES) {
122          break;
123       }
124    }
125
126    err = MP_OKAY;
127 LBL_ERR:
128    mp_clear(&b);
129    return err;
130 }
131
132 #endif