/src/dropbear/libtommath/bn_mp_prime_is_prime.c

Source (jump to first uncovered line)
#include "tommath_private.h"
#ifdef BN_MP_PRIME_IS_PRIME_C
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
/* SPDX-License-Identifier: Unlicense */

/* portable integer log of two with small footprint */
static unsigned int s_floor_ilog2(int value)
{
   unsigned int r = 0;
   while ((value >>= 1) != 0) {
      r++;
   }
   return r;
}


mp_err mp_prime_is_prime(const mp_int *a, int t, mp_bool *result)
{
   mp_int  b;
   int     ix, p_max = 0, size_a, len;
   mp_bool res;
   mp_err  err;
   unsigned int fips_rand, mask;

   /* default to no */
   *result = MP_NO;

   /* Some shortcuts */
   /* N > 3 */
   if (a->used == 1) {
      if ((a->dp[0] == 0u) || (a->dp[0] == 1u)) {
         *result = MP_NO;
         return MP_OKAY;
      }
      if (a->dp[0] == 2u) {
         *result = MP_YES;
         return MP_OKAY;
      }
   }

   /* N must be odd */
   if (MP_IS_EVEN(a)) {
      return MP_OKAY;
   }
   /* N is not a perfect square: floor(sqrt(N))^2 != N */
   if ((err = mp_is_square(a, &res)) != MP_OKAY) {
      return err;
   }
   if (res != MP_NO) {
      return MP_OKAY;
   }

   /* is the input equal to one of the primes in the table? */
   for (ix = 0; ix < PRIVATE_MP_PRIME_TAB_SIZE; ix++) {
      if (mp_cmp_d(a, s_mp_prime_tab[ix]) == MP_EQ) {
         *result = MP_YES;
         return MP_OKAY;
      }
   }
#ifdef MP_8BIT
   /* The search in the loop above was exhaustive in this case */
   if ((a->used == 1) && (PRIVATE_MP_PRIME_TAB_SIZE >= 31)) {
      return MP_OKAY;
   }
#endif

   /* first perform trial division */
   if ((err = s_mp_prime_is_divisible(a, &res)) != MP_OKAY) {
      return err;
   }

   /* return if it was trivially divisible */
   if (res == MP_YES) {
      return MP_OKAY;
   }

   /*
       Run the Miller-Rabin test with base 2 for the BPSW test.
    */
   if ((err = mp_init_set(&b, 2uL)) != MP_OKAY) {
      return err;
   }

   if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
      goto LBL_B;
   }
   if (res == MP_NO) {
      goto LBL_B;
   }
   /*
      Rumours have it that Mathematica does a second M-R test with base 3.
      Other rumours have it that their strong L-S test is slightly different.
      It does not hurt, though, beside a bit of extra runtime.
   */
   b.dp[0]++;
   if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
      goto LBL_B;
   }
   if (res == MP_NO) {
      goto LBL_B;
   }

   /*
    * Both, the Frobenius-Underwood test and the the Lucas-Selfridge test are quite
    * slow so if speed is an issue, define LTM_USE_ONLY_MR to use M-R tests with
    * bases 2, 3 and t random bases.
    */
#ifndef LTM_USE_ONLY_MR
   if (t >= 0) {
      /*
       * Use a Frobenius-Underwood test instead of the Lucas-Selfridge test for
       * MP_8BIT (It is unknown if the Lucas-Selfridge test works with 16-bit
       * integers but the necesssary analysis is on the todo-list).
       */
#if defined (MP_8BIT) || defined (LTM_USE_FROBENIUS_TEST)
      err = mp_prime_frobenius_underwood(a, &res);
      if ((err != MP_OKAY) && (err != MP_ITER)) {
         goto LBL_B;
      }
      if (res == MP_NO) {
         goto LBL_B;
      }
#else
      if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
         goto LBL_B;
      }
      if (res == MP_NO) {
         goto LBL_B;
      }
#endif
   }
#endif

   /* run at least one Miller-Rabin test with a random base */
   if (t == 0) {
      t = 1;
   }

   /*
      Only recommended if the input range is known to be < 3317044064679887385961981

      It uses the bases necessary for a deterministic M-R test if the input is
      smaller than  3317044064679887385961981
      The caller has to check the size.
      TODO: can be made a bit finer grained but comparing is not free.
   */
   if (t < 0) {
      /*
          Sorenson, Jonathan; Webster, Jonathan (2015).
           "Strong Pseudoprimes to Twelve Prime Bases".
       */
      /* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
      if ((err =   mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
         goto LBL_B;
      }

      if (mp_cmp(a, &b) == MP_LT) {
         p_max = 12;
      } else {
         /* 0x2be6951adc5b22410a5fd = 3317044064679887385961981 */
         if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
            goto LBL_B;
         }

         if (mp_cmp(a, &b) == MP_LT) {
            p_max = 13;
         } else {
            err = MP_VAL;
            goto LBL_B;
         }
      }

      /* we did bases 2 and 3  already, skip them */
      for (ix = 2; ix < p_max; ix++) {
         mp_set(&b, s_mp_prime_tab[ix]);
         if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
            goto LBL_B;
         }
         if (res == MP_NO) {
            goto LBL_B;
         }
      }
   }
   /*
       Do "t" M-R tests with random bases between 3 and "a".
       See Fips 186.4 p. 126ff
   */
   else if (t > 0) {
      /*
       * The mp_digit's have a defined bit-size but the size of the
       * array a.dp is a simple 'int' and this library can not assume full
       * compliance to the current C-standard (ISO/IEC 9899:2011) because
       * it gets used for small embeded processors, too. Some of those MCUs
       * have compilers that one cannot call standard compliant by any means.
       * Hence the ugly type-fiddling in the following code.
       */
      size_a = mp_count_bits(a);
      mask = (1u << s_floor_ilog2(size_a)) - 1u;
      /*
         Assuming the General Rieman hypothesis (never thought to write that in a
         comment) the upper bound can be lowered to  2*(log a)^2.
         E. Bach, "Explicit bounds for primality testing and related problems,"
         Math. Comp. 55 (1990), 355-380.

            size_a = (size_a/10) * 7;
            len = 2 * (size_a * size_a);

         E.g.: a number of size 2^2048 would be reduced to the upper limit

            floor(2048/10)*7 = 1428
            2 * 1428^2       = 4078368

         (would have been ~4030331.9962 with floats and natural log instead)
         That number is smaller than 2^28, the default bit-size of mp_digit.
      */

      /*
        How many tests, you might ask? Dana Jacobsen of Math::Prime::Util fame
        does exactly 1. In words: one. Look at the end of _GMP_is_prime() in
        Math-Prime-Util-GMP-0.50/primality.c if you do not believe it.

        The function mp_rand() goes to some length to use a cryptographically
        good PRNG. That also means that the chance to always get the same base
        in the loop is non-zero, although very low.
        If the BPSW test and/or the addtional Frobenious test have been
        performed instead of just the Miller-Rabin test with the bases 2 and 3,
        a single extra test should suffice, so such a very unlikely event
        will not do much harm.

        To preemptivly answer the dangling question: no, a witness does not
        need to be prime.
      */
      for (ix = 0; ix < t; ix++) {
         /* mp_rand() guarantees the first digit to be non-zero */
         if ((err = mp_rand(&b, 1)) != MP_OKAY) {
            goto LBL_B;
         }
         /*
          * Reduce digit before casting because mp_digit might be bigger than
          * an unsigned int and "mask" on the other side is most probably not.
          */
         fips_rand = (unsigned int)(b.dp[0] & (mp_digit) mask);
#ifdef MP_8BIT
         /*
          * One 8-bit digit is too small, so concatenate two if the size of
          * unsigned int allows for it.
          */
         if ((MP_SIZEOF_BITS(unsigned int)/2) >= MP_SIZEOF_BITS(mp_digit)) {
            if ((err = mp_rand(&b, 1)) != MP_OKAY) {
               goto LBL_B;
            }
            fips_rand <<= MP_SIZEOF_BITS(mp_digit);
            fips_rand |= (unsigned int) b.dp[0];
            fips_rand &= mask;
         }
#endif
         if (fips_rand > (unsigned int)(INT_MAX - MP_DIGIT_BIT)) {
            len = INT_MAX / MP_DIGIT_BIT;
         } else {
            len = (((int)fips_rand + MP_DIGIT_BIT) / MP_DIGIT_BIT);
         }
         /*  Unlikely. */
         if (len < 0) {
            ix--;
            continue;
         }
         /*
          * As mentioned above, one 8-bit digit is too small and
          * although it can only happen in the unlikely case that
          * an "unsigned int" is smaller than 16 bit a simple test
          * is cheap and the correction even cheaper.
          */
#ifdef MP_8BIT
         /* All "a" < 2^8 have been caught before */
         if (len == 1) {
            len++;
         }
#endif
         if ((err = mp_rand(&b, len)) != MP_OKAY) {
            goto LBL_B;
         }
         /*
          * That number might got too big and the witness has to be
          * smaller than "a"
          */
         len = mp_count_bits(&b);
         if (len >= size_a) {
            len = (len - size_a) + 1;
            if ((err = mp_div_2d(&b, len, &b, NULL)) != MP_OKAY) {
               goto LBL_B;
            }
         }
         /* Although the chance for b <= 3 is miniscule, try again. */
         if (mp_cmp_d(&b, 3uL) != MP_GT) {
            ix--;
            continue;
         }
         if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
            goto LBL_B;
         }
         if (res == MP_NO) {
            goto LBL_B;
         }
      }
   }

   /* passed the test */
   *result = MP_YES;
LBL_B:
   mp_clear(&b);
   return err;
}

#endif

Coverage Report

Created: 2025-08-29 06:39

Line	Count	Source (jump to first uncovered line)
1		#include "tommath_private.h"
2		#ifdef BN_MP_PRIME_IS_PRIME_C
3		/* LibTomMath, multiple-precision integer library -- Tom St Denis */
4		/* SPDX-License-Identifier: Unlicense */
5
6		/* portable integer log of two with small footprint */
7		static unsigned int s_floor_ilog2(int value)
8	52	{
9	52	unsigned int r = 0;
10	416	while ((value >>= 1) != 0) {
11	364	r++;
12	364	}
13	52	return r;
14	52	}
15
16
17		mp_err mp_prime_is_prime(const mp_int a, int t, mp_bool result)
18	218	{
19	218	mp_int b;
20	218	int ix, p_max = 0, size_a, len;
21	218	mp_bool res;
22	218	mp_err err;
23	218	unsigned int fips_rand, mask;
24
25		/* default to no */
26	218	*result = MP_NO;
27
28		/* Some shortcuts */
29		/* N > 3 */
30	218	if (a->used == 1) {
31	0	if ((a->dp[0] == 0u) \|\| (a->dp[0] == 1u)) {
32	0	*result = MP_NO;
33	0	return MP_OKAY;
34	0	}
35	0	if (a->dp[0] == 2u) {
36	0	*result = MP_YES;
37	0	return MP_OKAY;
38	0	}
39	0	}
40
41		/* N must be odd */
42	218	if (MP_IS_EVEN(a)) {
43	82	return MP_OKAY;
44	82	}
45		/* N is not a perfect square: floor(sqrt(N))^2 != N */
46	136	if ((err = mp_is_square(a, &res)) != MP_OKAY) {
47	0	return err;
48	0	}
49	136	if (res != MP_NO) {
50	1	return MP_OKAY;
51	1	}
52
53		/* is the input equal to one of the primes in the table? */
54	34.6k	for (ix = 0; ix < PRIVATE_MP_PRIME_TAB_SIZE; ix++) {
55	34.5k	if (mp_cmp_d(a, s_mp_prime_tab[ix]) == MP_EQ) {
56	0	*result = MP_YES;
57	0	return MP_OKAY;
58	0	}
59	34.5k	}
60		#ifdef MP_8BIT
61		/* The search in the loop above was exhaustive in this case */
62		if ((a->used == 1) && (PRIVATE_MP_PRIME_TAB_SIZE >= 31)) {
63		return MP_OKAY;
64		}
65		#endif
66
67		/* first perform trial division */
68	135	if ((err = s_mp_prime_is_divisible(a, &res)) != MP_OKAY) {
69	0	return err;
70	0	}
71
72		/* return if it was trivially divisible */
73	135	if (res == MP_YES) {
74	46	return MP_OKAY;
75	46	}
76
77		/*
78		Run the Miller-Rabin test with base 2 for the BPSW test.
79		*/
80	89	if ((err = mp_init_set(&b, 2uL)) != MP_OKAY) {
81	0	return err;
82	0	}
83
84	89	if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
85	0	goto LBL_B;
86	0	}
87	89	if (res == MP_NO) {
88	37	goto LBL_B;
89	37	}
90		/*
91		Rumours have it that Mathematica does a second M-R test with base 3.
92		Other rumours have it that their strong L-S test is slightly different.
93		It does not hurt, though, beside a bit of extra runtime.
94		*/
95	52	b.dp[0]++;
96	52	if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
97	0	goto LBL_B;
98	0	}
99	52	if (res == MP_NO) {
100	0	goto LBL_B;
101	0	}
102
103		/*
104		* Both, the Frobenius-Underwood test and the the Lucas-Selfridge test are quite
105		* slow so if speed is an issue, define LTM_USE_ONLY_MR to use M-R tests with
106		* bases 2, 3 and t random bases.
107		*/
108		#ifndef LTM_USE_ONLY_MR
109		if (t >= 0) {
110		/*
111		* Use a Frobenius-Underwood test instead of the Lucas-Selfridge test for
112		* MP_8BIT (It is unknown if the Lucas-Selfridge test works with 16-bit
113		* integers but the necesssary analysis is on the todo-list).
114		*/
115		#if defined (MP_8BIT) \|\| defined (LTM_USE_FROBENIUS_TEST)
116		err = mp_prime_frobenius_underwood(a, &res);
117		if ((err != MP_OKAY) && (err != MP_ITER)) {
118		goto LBL_B;
119		}
120		if (res == MP_NO) {
121		goto LBL_B;
122		}
123		#else
124		if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
125		goto LBL_B;
126		}
127		if (res == MP_NO) {
128		goto LBL_B;
129		}
130		#endif
131		}
132		#endif
133
134		/* run at least one Miller-Rabin test with a random base */
135	52	if (t == 0) {
136	0	t = 1;
137	0	}
138
139		/*
140		Only recommended if the input range is known to be < 3317044064679887385961981
141
142		It uses the bases necessary for a deterministic M-R test if the input is
143		smaller than 3317044064679887385961981
144		The caller has to check the size.
145		TODO: can be made a bit finer grained but comparing is not free.
146		*/
147	52	if (t < 0) {
148		/*
149		Sorenson, Jonathan; Webster, Jonathan (2015).
150		"Strong Pseudoprimes to Twelve Prime Bases".
151		*/
152		/* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
153	0	if ((err = mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
154	0	goto LBL_B;
155	0	}
156
157	0	if (mp_cmp(a, &b) == MP_LT) {
158	0	p_max = 12;
159	0	} else {
160		/* 0x2be6951adc5b22410a5fd = 3317044064679887385961981 */
161	0	if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
162	0	goto LBL_B;
163	0	}
164
165	0	if (mp_cmp(a, &b) == MP_LT) {
166	0	p_max = 13;
167	0	} else {
168	0	err = MP_VAL;
169	0	goto LBL_B;
170	0	}
171	0	}
172
173		/* we did bases 2 and 3 already, skip them */
174	0	for (ix = 2; ix < p_max; ix++) {
175	0	mp_set(&b, s_mp_prime_tab[ix]);
176	0	if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
177	0	goto LBL_B;
178	0	}
179	0	if (res == MP_NO) {
180	0	goto LBL_B;
181	0	}
182	0	}
183	0	}
184		/*
185		Do "t" M-R tests with random bases between 3 and "a".
186		See Fips 186.4 p. 126ff
187		*/
188	52	else if (t > 0) {
189		/*
190		* The mp_digit's have a defined bit-size but the size of the
191		* array a.dp is a simple 'int' and this library can not assume full
192		* compliance to the current C-standard (ISO/IEC 9899:2011) because
193		* it gets used for small embeded processors, too. Some of those MCUs
194		* have compilers that one cannot call standard compliant by any means.
195		* Hence the ugly type-fiddling in the following code.
196		*/
197	52	size_a = mp_count_bits(a);
198	52	mask = (1u << s_floor_ilog2(size_a)) - 1u;
199		/*
200		Assuming the General Rieman hypothesis (never thought to write that in a
201		comment) the upper bound can be lowered to 2*(log a)^2.
202		E. Bach, "Explicit bounds for primality testing and related problems,"
203		Math. Comp. 55 (1990), 355-380.
204
205		size_a = (size_a/10) * 7;
206		len = 2 * (size_a * size_a);
207
208		E.g.: a number of size 2^2048 would be reduced to the upper limit
209
210		floor(2048/10)*7 = 1428
211		2 * 1428^2 = 4078368
212
213		(would have been ~4030331.9962 with floats and natural log instead)
214		That number is smaller than 2^28, the default bit-size of mp_digit.
215		*/
216
217		/*
218		How many tests, you might ask? Dana Jacobsen of Math::Prime::Util fame
219		does exactly 1. In words: one. Look at the end of _GMP_is_prime() in
220		Math-Prime-Util-GMP-0.50/primality.c if you do not believe it.
221
222		The function mp_rand() goes to some length to use a cryptographically
223		good PRNG. That also means that the chance to always get the same base
224		in the loop is non-zero, although very low.
225		If the BPSW test and/or the addtional Frobenious test have been
226		performed instead of just the Miller-Rabin test with the bases 2 and 3,
227		a single extra test should suffice, so such a very unlikely event
228		will not do much harm.
229
230		To preemptivly answer the dangling question: no, a witness does not
231		need to be prime.
232		*/
233	1.87k	for (ix = 0; ix < t; ix++) {
234		/* mp_rand() guarantees the first digit to be non-zero */
235	1.82k	if ((err = mp_rand(&b, 1)) != MP_OKAY) {
236	0	goto LBL_B;
237	0	}
238		/*
239		* Reduce digit before casting because mp_digit might be bigger than
240		* an unsigned int and "mask" on the other side is most probably not.
241		*/
242	1.82k	fips_rand = (unsigned int)(b.dp[0] & (mp_digit) mask);
243		#ifdef MP_8BIT
244		/*
245		* One 8-bit digit is too small, so concatenate two if the size of
246		* unsigned int allows for it.
247		*/
248		if ((MP_SIZEOF_BITS(unsigned int)/2) >= MP_SIZEOF_BITS(mp_digit)) {
249		if ((err = mp_rand(&b, 1)) != MP_OKAY) {
250		goto LBL_B;
251		}
252		fips_rand <<= MP_SIZEOF_BITS(mp_digit);
253		fips_rand \|= (unsigned int) b.dp[0];
254		fips_rand &= mask;
255		}
256		#endif
257	1.82k	if (fips_rand > (unsigned int)(INT_MAX - MP_DIGIT_BIT)) {
258	0	len = INT_MAX / MP_DIGIT_BIT;
259	1.82k	} else {
260	1.82k	len = (((int)fips_rand + MP_DIGIT_BIT) / MP_DIGIT_BIT);
261	1.82k	}
262		/* Unlikely. */
263	1.82k	if (len < 0) {
264	0	ix--;
265	0	continue;
266	0	}
267		/*
268		* As mentioned above, one 8-bit digit is too small and
269		* although it can only happen in the unlikely case that
270		* an "unsigned int" is smaller than 16 bit a simple test
271		* is cheap and the correction even cheaper.
272		*/
273		#ifdef MP_8BIT
274		/* All "a" < 2^8 have been caught before */
275		if (len == 1) {
276		len++;
277		}
278		#endif
279	1.82k	if ((err = mp_rand(&b, len)) != MP_OKAY) {
280	0	goto LBL_B;
281	0	}
282		/*
283		* That number might got too big and the witness has to be
284		* smaller than "a"
285		*/
286	1.82k	len = mp_count_bits(&b);
287	1.82k	if (len >= size_a) {
288	203	len = (len - size_a) + 1;
289	203	if ((err = mp_div_2d(&b, len, &b, NULL)) != MP_OKAY) {
290	0	goto LBL_B;
291	0	}
292	203	}
293		/* Although the chance for b <= 3 is miniscule, try again. */
294	1.82k	if (mp_cmp_d(&b, 3uL) != MP_GT) {
295	0	ix--;
296	0	continue;
297	0	}
298	1.82k	if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
299	0	goto LBL_B;
300	0	}
301	1.82k	if (res == MP_NO) {
302	0	goto LBL_B;
303	0	}
304	1.82k	}
305	52	}
306
307		/* passed the test */
308	52	*result = MP_YES;
309	89	LBL_B:
310	89	mp_clear(&b);
311	89	return err;
312	52	}
313
314		#endif