From 4f809ca13e7d484837de91cc53331be2ce14a66b Mon Sep 17 00:00:00 2001 From: Christoph Zurnieden Date: Fri, 18 Sep 2026 20:22:03 +0200 Subject: [PATCH] Addition of a small segmented sieve --- demo/test.c | 167 +++++++++++++++++++++ doc/bn.tex | 42 ++++++ doc/tommath.3 | 33 +++++ helper.pl | 2 +- libtommath_VS2008.vcproj | 24 +++ makefile | 22 +-- makefile.mingw | 22 +-- makefile.msvc | 22 +-- makefile.shared | 22 +-- makefile.unix | 22 +-- mp_small_prime_sieve_clear.c | 24 +++ mp_small_prime_sieve_init.c | 25 ++++ mp_small_prime_sieve_is_small_prime.c | 109 ++++++++++++++ mp_small_prime_sieve_next_prime.c | 58 ++++++++ mp_small_prime_sieve_prec_prime.c | 48 ++++++ s_mp_small_prime_sieve.c | 205 ++++++++++++++++++++++++++ sources.cmake | 6 + tommath.def | 5 + tommath.h | 42 ++++++ tommath_class.h | 41 ++++++ tommath_private.h | 20 +++ 21 files changed, 910 insertions(+), 51 deletions(-) create mode 100644 mp_small_prime_sieve_clear.c create mode 100644 mp_small_prime_sieve_init.c create mode 100644 mp_small_prime_sieve_is_small_prime.c create mode 100644 mp_small_prime_sieve_next_prime.c create mode 100644 mp_small_prime_sieve_prec_prime.c create mode 100644 s_mp_small_prime_sieve.c diff --git a/demo/test.c b/demo/test.c index 2fa6e08d..bdc5747c 100644 --- a/demo/test.c +++ b/demo/test.c @@ -469,6 +469,170 @@ static int test_mp_signed_rsh(void) } +static int test_mp_small_prime_sieve_is_small_prime(void) +{ + mp_erat_sieve sieve; + mp_int N; + uint32_t n; + int i = 0; + bool isprime = false, isreallyprime = false; + + /* Cold start to test a cold start */ + DO(mp_small_prime_sieve_init(&sieve, false)); + DOR(mp_init(&N)); + /* Generate a handful of random numbers 1 4294967291u)) { + i--; + continue; + } + DO(mp_small_prime_sieve_is_small_prime(n, &isprime, &sieve)); + + mp_set_u32(&N, n); + DO(mp_prime_is_prime(&N, 1, &isreallyprime)); + EXPECT(isprime == isreallyprime); + } + mp_small_prime_sieve_clear(&sieve); + mp_clear(&N); + return EXIT_SUCCESS; +LBL_ERR: + mp_small_prime_sieve_clear(&sieve); + mp_clear(&N); + return EXIT_FAILURE; +} + + + + +static int test_mp_small_prime_sieve_prec_prime(void) +{ + mp_erat_sieve sieve; + mp_int N; +#ifdef MP_SMALL_SIEVE_FULL_ROUND_DECREASING + uint64_t primesum; + uint32_t prime; +#endif + uint32_t n = 0, precprime, nextprime ; + int i = 0; + bool isreallyprime = false; + + /* Warm start to test a warm start */ + DO(mp_small_prime_sieve_init(&sieve, true)); + DOR(mp_init(&N)); + + for (i = 0; i < 100; i++) { + /* rand_uint() is a PRNG so it won't loop forever */ + n = (uint32_t)rand_uint(); + if (n < 3) { + i--; + continue; + } + DO(mp_small_prime_sieve_prec_prime(n - 1, &precprime, &sieve)); + mp_set_u32(&N, precprime); + DO(mp_prime_is_prime(&N, 1, &isreallyprime)); + EXPECT(isreallyprime); + + /* precprime <= n, so nextprime(precprime) >= n */ + DO(mp_prime_next_prime(&N, 1, false)); + nextprime = mp_get_u32(&N); + EXPECT(nextprime >= n); + } + + /* + A full round takes quite some time, but is checking everything + A long time is about 30 seconds, but with valgrind about 30 minutes! + */ +#ifdef MP_SMALL_SIEVE_FULL_ROUND_DECREASING + prime = 4294967291u; + primesum = 0u; + while (prime >= 2) { + primesum = primesum + (uint64_t)prime; + DO(mp_small_prime_sieve_prec_prime(prime - 1, &prime, &sieve)); + } + EXPECT(primesum == 425649736193687430u); +#endif + mp_small_prime_sieve_clear(&sieve); + mp_clear(&N); + return EXIT_SUCCESS; +LBL_ERR: + mp_small_prime_sieve_clear(&sieve); + mp_clear(&N); + return EXIT_FAILURE; +} + + + + + +static int test_mp_small_prime_sieve_next_prime(void) +{ + mp_erat_sieve sieve; + mp_int N; + mp_err err2; +#ifdef MP_SMALL_SIEVE_FULL_ROUND_INCREASING + uint64_t primesum; + uint32_t prime; +#endif + uint32_t n, nextprime = 0u, reallynextprime = 0u; + int i = 0; + bool isreallyprime = false; + + /* Cold start to test a cold start */ + DO(mp_small_prime_sieve_init(&sieve, false)); + DOR(mp_init(&N)); + + for (i = 0; i < 100; i++) { + /* rand_uint() is a PRNG so it won't loop forever */ + n = (uint32_t)rand_uint(); + if ((n < 2) || (n > 4294967279u)) { + i--; + continue; + } + if ((err2 = mp_small_prime_sieve_next_prime(n + 1, &nextprime, &sieve)) != MP_OKAY) { + if (err2 != MP_OVF) { + goto LBL_ERR; + } else { + i--; + continue; + } + } + mp_set_u32(&N, nextprime); + DO(mp_prime_is_prime(&N, 1, &isreallyprime)); + EXPECT(isreallyprime); + + mp_set_u32(&N, n); + DO(mp_prime_next_prime(&N, 1, false)); + reallynextprime = mp_get_u32(&N); + EXPECT(nextprime == reallynextprime); + } + + /* A full round takes quite some time, but is checking everything */ +#ifdef MP_SMALL_SIEVE_FULL_ROUND_INCREASING + prime = 0u; + primesum = 0u; + while (prime <= 4294967291u) { + if ((err2 = mp_small_prime_sieve_next_prime(prime + 1, &prime, &sieve)) != MP_OKAY) { + /* Will throw ERAT_OVL which is the place we stop */ + if (err2 != MP_OVF) { + goto LBL_ERR; + } else { + break; + } + } + primesum = primesum + (uint64_t)prime; + } + EXPECT(primesum == 425649736193687430u); +#endif + mp_small_prime_sieve_clear(&sieve); + mp_clear(&N); + return EXIT_SUCCESS; +LBL_ERR: + mp_small_prime_sieve_clear(&sieve); + mp_clear(&N); + return EXIT_FAILURE; +} + static int test_mp_xor(void) { int i; @@ -2594,6 +2758,9 @@ static int unit_tests(int argc, char **argv) T1(mp_set_double, MP_SET_DOUBLE), #endif T1(mp_signed_rsh, MP_SIGNED_RSH), + T2(mp_small_prime_sieve_is_small_prime, MP_SMALL_PRIME_SIEVE_IS_SMALL_PRIME, S_MP_SMALL_PRIME_SIEVE), + T2(mp_small_prime_sieve_prec_prime, MP_SMALL_PRIME_SIEVE_PREC_PRIME, S_MP_SMALL_PRIME_SIEVE), + T2(mp_small_prime_sieve_next_prime, MP_SMALL_PRIME_SIEVE_NEXT_PRIME, S_MP_SMALL_PRIME_SIEVE), T2(mp_sqrt, MP_SQRT, MP_ROOT_N), T1(mp_sqrtmod_prime, MP_SQRTMOD_PRIME), T1(mp_xor, MP_XOR), diff --git a/doc/bn.tex b/doc/bn.tex index 63e71633..c4627b01 100644 --- a/doc/bn.tex +++ b/doc/bn.tex @@ -2396,6 +2396,47 @@ \section{Random Primes} \label{fig:primeopts} \end{figure} +\subsection{Small Prime Sieve} +A small segmented sieve (unsigned 32-bit max.) is available via three functions. +The type {\texttt{ERAT\_UINT}}, defined (\emph{not} {\texttt{typedef}}) in {\texttt{tommath.h}} +is normally {\texttt{uint32\_t}} but should be checked before use to be sure. + + +\index{mp\_small\_prime\_sieve\_init} +\begin{alltt} +mp_err mp_small_prime_sieve_init(mp_erat_sieve *sieve, bool warmup); +\end{alltt} +Initialize sieve to default values, build base-sieve, too, if the parameter {\texttt{warmup}} is true +otherwise the first call to either {\texttt{mp\_small\_prime\_sieve\_is\_small\_prime}}, +{\texttt{mp\_small\_prime\_sieve\_next\_prime}}, or {\texttt{mp\_small\_prime\_sieve\_prec\_prime}} +do it at first call. + +\index{mp\_small\_prime\_sieve\_clear} +\begin{alltt} +mp_err mp_small_prime_sieve_clear(mp_erat_sieve *sieve); +\end{alltt} +Reset sieve and free memory . + +\index{mp\_small\_prime\_sieve\_is\_small\_prime} +\begin{alltt} +mp_err mp_small_prime_sieve_is_small_prime(ERAT_UINT n, bool *result, mp_erat_sieve *sieve); +\end{alltt} +Check primality of $n$ if $n < 2^32$. + + +\index{mp\_small\_prime\_sieve\_next\_prime} +\begin{alltt} +mp_err mp_small_prime_sieve_next_prime(ERAT_UINT n, ERAT_UINT *result, mp_erat_sieve *sieve); +\end{alltt} +Sets {\texttt{result}} to the prime {\texttt{p}} larger than or equal to {\texttt{n}} if $p > 2^32 - 5$. +Returns {\texttt{MP\_OVF}} if $p > 2^32 - 5$ as an information about that fact. + +\index{mp\_small\_prime\_sieve\_prec\_prime} +\begin{alltt} +mp_err mp_small_prime_sieve_prec_prime(ERAT_UINT n, ERAT_UINT *result, mp_erat_sieve *sieve); +\end{alltt} +Sets {\texttt{result}} to the prime {\texttt{p}} smaller than or equal to {\texttt{n}} if $p > 1$. + \chapter{Random Number Generation} \section{PRNG} \index{mp\_rand} @@ -2415,6 +2456,7 @@ \section{PRNG} void mp_rand_source(mp_err(*source)(void *out, size_t size)); \end{alltt} + \chapter{Input and Output} \section{ASCII Conversions} \subsection{To ASCII} diff --git a/doc/tommath.3 b/doc/tommath.3 index fabaaf5d..c4af914d 100644 --- a/doc/tommath.3 +++ b/doc/tommath.3 @@ -1137,6 +1137,39 @@ This method uses \fBrealloc\fP with all its problems, use with caution. Shifts right \fBa\fP by \fBb\fP bits with sign extension. .in -1i +.LP +.BI "mp_err mp_small_prime_sieve_init(mp_erat_sieve *" sieve ", bool " warmup ")" +.in 1i +Init sieve to default values, build base-sieve, too, if the parameter \fBwarmup\fP is true +otherwise the first call to either \fBmp_small_prime_sieve_is_small_prime\fP, +\fBmp_small_prime_sieve_next_prime\fP, or \fBmp_small_prime_sieve_prec_prime\fP do it at first call. +.in -1i + +.LP +.BI "mp_err mp_small_prime_sieve_clear(mp_erat_sieve *" sieve ")" +.in 1i +Reset sieve and clear memory. +.in -1i + +.LP +.BI "mp_err mp_small_prime_sieve_is_small_prime(ERAT_UINT " n ", bool *" result ", mp_erat_sieve * " sieve ")" +.in 1i +Check primality of \fBn\fP if \fBn < 2^32\fP. +.in -1i + +.LP +.BI "mp_err mp_small_prime_sieve_next_prime(ERAT_UINT " n ", ERAT_UINT *" result ", mp_erat_sieve *" sieve ")" +.in 1i +Sets \fBresult\fP to the prime \fBp\fP larger than or equal to \fBp\fP if \fBp > 2^32\fP. +Returns \fBMP_OVF\fP if \fBp > 2^32\fP as an information about that fact. +.in -1i + +.LP +.BI "mp_err mp_small_prime_sieve_prec_prime(ERAT_UINT " n ", ERAT_UINT *" result ", mp_erat_sieve *" sieve ")" +.in 1i +Sets \fBresult\fP to the prime \fBp\fP smaller than or equal to \fBn\fP if \fBp > 1\fP. +.in -1i + .LP .BI "mp_err mp_sqr (const mp_int *" a ", mp_int *" b ")" .in 1i diff --git a/helper.pl b/helper.pl index 8a430f73..0df67c00 100755 --- a/helper.pl +++ b/helper.pl @@ -58,7 +58,7 @@ sub check_source { push @{$troubles->{unwanted_free}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bfree\s*\(/; # and we probably want to also avoid the following push @{$troubles->{unwanted_memcpy}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemcpy\s*\(/ && $file !~ /s_mp_copy_digs.c/; - push @{$troubles->{unwanted_memset}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemset\s*\(/ && $file !~ /s_mp_zero_buf.c/ && $file !~ /s_mp_zero_digs.c/; + push @{$troubles->{unwanted_memset}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemset\s*\(/ && $file !~ /s_mp_zero_buf.c/ && $file !~ /s_mp_zero_digs.c/ && $file !~ /s_mp_small_prime_sieve.c/ ; push @{$troubles->{unwanted_memmove}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemmove\s*\(/; push @{$troubles->{unwanted_memcmp}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemcmp\s*\(/; push @{$troubles->{unwanted_strcmp}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bstrcmp\s*\(/; diff --git a/libtommath_VS2008.vcproj b/libtommath_VS2008.vcproj index 080b7e0d..d3b00a56 100644 --- a/libtommath_VS2008.vcproj +++ b/libtommath_VS2008.vcproj @@ -748,6 +748,26 @@ RelativePath="mp_signed_rsh.c" > + + + + + + + + + + @@ -916,6 +936,10 @@ RelativePath="s_mp_rand_source.c" > + + diff --git a/makefile b/makefile index 6daa522b..d21b8301 100644 --- a/makefile +++ b/makefile @@ -41,16 +41,18 @@ mp_prime_rabin_miller_trials.o mp_prime_rand.o mp_prime_strong_lucas_selfridge.o mp_radix_size_overestimate.o mp_rand.o mp_rand_source.o mp_read_radix.o mp_reduce.o mp_reduce_2k.o \ mp_reduce_2k_l.o mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o mp_reduce_is_2k_l.o \ mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o \ -mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_sqrmod.o mp_sqrt.o \ -mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o \ -mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \ -s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \ -s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \ -s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \ -s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \ -s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \ -s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ -s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o +mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_small_prime_sieve_clear.o \ +mp_small_prime_sieve_init.o mp_small_prime_sieve_is_small_prime.o mp_small_prime_sieve_next_prime.o \ +mp_small_prime_sieve_prec_prime.o mp_sqrmod.o mp_sqrt.o mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o \ +mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o \ +s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o \ +s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o \ +s_mp_invmod_odd.o s_mp_log_2expt.o s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o \ +s_mp_mul_comba.o s_mp_mul_high.o s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o \ +s_mp_prime_is_divisible.o s_mp_prime_tab.o s_mp_radix_map.o s_mp_radix_size_overestimate.o \ +s_mp_rand_platform.o s_mp_rand_source.o s_mp_small_prime_sieve.o s_mp_sqr.o s_mp_sqr_comba.o \ +s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o s_mp_warray_put.o \ +s_mp_zero_buf.o s_mp_zero_digs.o #END_INS diff --git a/makefile.mingw b/makefile.mingw index 10f2bd59..c03e6b4d 100644 --- a/makefile.mingw +++ b/makefile.mingw @@ -43,16 +43,18 @@ mp_prime_rabin_miller_trials.o mp_prime_rand.o mp_prime_strong_lucas_selfridge.o mp_radix_size_overestimate.o mp_rand.o mp_rand_source.o mp_read_radix.o mp_reduce.o mp_reduce_2k.o \ mp_reduce_2k_l.o mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o mp_reduce_is_2k_l.o \ mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o \ -mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_sqrmod.o mp_sqrt.o \ -mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o \ -mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \ -s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \ -s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \ -s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \ -s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \ -s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \ -s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ -s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o +mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_small_prime_sieve_clear.o \ +mp_small_prime_sieve_init.o mp_small_prime_sieve_is_small_prime.o mp_small_prime_sieve_next_prime.o \ +mp_small_prime_sieve_prec_prime.o mp_sqrmod.o mp_sqrt.o mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o \ +mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o \ +s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o \ +s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o \ +s_mp_invmod_odd.o s_mp_log_2expt.o s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o \ +s_mp_mul_comba.o s_mp_mul_high.o s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o \ +s_mp_prime_is_divisible.o s_mp_prime_tab.o s_mp_radix_map.o s_mp_radix_size_overestimate.o \ +s_mp_rand_platform.o s_mp_rand_source.o s_mp_small_prime_sieve.o s_mp_sqr.o s_mp_sqr_comba.o \ +s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o s_mp_warray_put.o \ +s_mp_zero_buf.o s_mp_zero_digs.o HEADERS_PUB=tommath.h HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB) diff --git a/makefile.msvc b/makefile.msvc index f160f2cc..c3c9b725 100644 --- a/makefile.msvc +++ b/makefile.msvc @@ -39,16 +39,18 @@ mp_prime_rabin_miller_trials.obj mp_prime_rand.obj mp_prime_strong_lucas_selfrid mp_radix_size_overestimate.obj mp_rand.obj mp_rand_source.obj mp_read_radix.obj mp_reduce.obj mp_reduce_2k.obj \ mp_reduce_2k_l.obj mp_reduce_2k_setup.obj mp_reduce_2k_setup_l.obj mp_reduce_is_2k.obj mp_reduce_is_2k_l.obj \ mp_reduce_setup.obj mp_root_n.obj mp_rshd.obj mp_sbin_size.obj mp_set.obj mp_set_double.obj mp_set_i32.obj mp_set_i64.obj \ -mp_set_l.obj mp_set_u32.obj mp_set_u64.obj mp_set_ul.obj mp_shrink.obj mp_signed_rsh.obj mp_sqrmod.obj mp_sqrt.obj \ -mp_sqrtmod_prime.obj mp_sub.obj mp_sub_d.obj mp_submod.obj mp_to_radix.obj mp_to_sbin.obj mp_to_ubin.obj mp_ubin_size.obj \ -mp_unpack.obj mp_warray_free.obj mp_xor.obj mp_zero.obj s_mp_add.obj s_mp_copy_digs.obj s_mp_div_3.obj \ -s_mp_div_recursive.obj s_mp_div_school.obj s_mp_div_small.obj s_mp_exptmod.obj s_mp_exptmod_fast.obj s_mp_fp_log.obj \ -s_mp_fp_log_d.obj s_mp_get_bit.obj s_mp_invmod.obj s_mp_invmod_odd.obj s_mp_log_2expt.obj \ -s_mp_montgomery_reduce_comba.obj s_mp_mul.obj s_mp_mul_balance.obj s_mp_mul_comba.obj s_mp_mul_high.obj \ -s_mp_mul_high_comba.obj s_mp_mul_karatsuba.obj s_mp_mul_toom.obj s_mp_prime_is_divisible.obj s_mp_prime_tab.obj \ -s_mp_radix_map.obj s_mp_radix_size_overestimate.obj s_mp_rand_platform.obj s_mp_rand_source.obj s_mp_sqr.obj \ -s_mp_sqr_comba.obj s_mp_sqr_karatsuba.obj s_mp_sqr_toom.obj s_mp_sub.obj s_mp_warray.obj s_mp_warray_get.obj \ -s_mp_warray_put.obj s_mp_zero_buf.obj s_mp_zero_digs.obj +mp_set_l.obj mp_set_u32.obj mp_set_u64.obj mp_set_ul.obj mp_shrink.obj mp_signed_rsh.obj mp_small_prime_sieve_clear.obj \ +mp_small_prime_sieve_init.obj mp_small_prime_sieve_is_small_prime.obj mp_small_prime_sieve_next_prime.obj \ +mp_small_prime_sieve_prec_prime.obj mp_sqrmod.obj mp_sqrt.obj mp_sqrtmod_prime.obj mp_sub.obj mp_sub_d.obj mp_submod.obj \ +mp_to_radix.obj mp_to_sbin.obj mp_to_ubin.obj mp_ubin_size.obj mp_unpack.obj mp_warray_free.obj mp_xor.obj mp_zero.obj \ +s_mp_add.obj s_mp_copy_digs.obj s_mp_div_3.obj s_mp_div_recursive.obj s_mp_div_school.obj s_mp_div_small.obj \ +s_mp_exptmod.obj s_mp_exptmod_fast.obj s_mp_fp_log.obj s_mp_fp_log_d.obj s_mp_get_bit.obj s_mp_invmod.obj \ +s_mp_invmod_odd.obj s_mp_log_2expt.obj s_mp_montgomery_reduce_comba.obj s_mp_mul.obj s_mp_mul_balance.obj \ +s_mp_mul_comba.obj s_mp_mul_high.obj s_mp_mul_high_comba.obj s_mp_mul_karatsuba.obj s_mp_mul_toom.obj \ +s_mp_prime_is_divisible.obj s_mp_prime_tab.obj s_mp_radix_map.obj s_mp_radix_size_overestimate.obj \ +s_mp_rand_platform.obj s_mp_rand_source.obj s_mp_small_prime_sieve.obj s_mp_sqr.obj s_mp_sqr_comba.obj \ +s_mp_sqr_karatsuba.obj s_mp_sqr_toom.obj s_mp_sub.obj s_mp_warray.obj s_mp_warray_get.obj s_mp_warray_put.obj \ +s_mp_zero_buf.obj s_mp_zero_digs.obj HEADERS_PUB=tommath.h HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB) diff --git a/makefile.shared b/makefile.shared index 20bd8e03..b1a5527c 100644 --- a/makefile.shared +++ b/makefile.shared @@ -38,16 +38,18 @@ mp_prime_rabin_miller_trials.o mp_prime_rand.o mp_prime_strong_lucas_selfridge.o mp_radix_size_overestimate.o mp_rand.o mp_rand_source.o mp_read_radix.o mp_reduce.o mp_reduce_2k.o \ mp_reduce_2k_l.o mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o mp_reduce_is_2k_l.o \ mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o \ -mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_sqrmod.o mp_sqrt.o \ -mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o \ -mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \ -s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \ -s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \ -s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \ -s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \ -s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \ -s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ -s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o +mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_small_prime_sieve_clear.o \ +mp_small_prime_sieve_init.o mp_small_prime_sieve_is_small_prime.o mp_small_prime_sieve_next_prime.o \ +mp_small_prime_sieve_prec_prime.o mp_sqrmod.o mp_sqrt.o mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o \ +mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o \ +s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o \ +s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o \ +s_mp_invmod_odd.o s_mp_log_2expt.o s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o \ +s_mp_mul_comba.o s_mp_mul_high.o s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o \ +s_mp_prime_is_divisible.o s_mp_prime_tab.o s_mp_radix_map.o s_mp_radix_size_overestimate.o \ +s_mp_rand_platform.o s_mp_rand_source.o s_mp_small_prime_sieve.o s_mp_sqr.o s_mp_sqr_comba.o \ +s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o s_mp_warray_put.o \ +s_mp_zero_buf.o s_mp_zero_digs.o #END_INS diff --git a/makefile.unix b/makefile.unix index 06dadddb..bff49fcb 100644 --- a/makefile.unix +++ b/makefile.unix @@ -44,16 +44,18 @@ mp_prime_rabin_miller_trials.o mp_prime_rand.o mp_prime_strong_lucas_selfridge.o mp_radix_size_overestimate.o mp_rand.o mp_rand_source.o mp_read_radix.o mp_reduce.o mp_reduce_2k.o \ mp_reduce_2k_l.o mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o mp_reduce_is_2k_l.o \ mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o \ -mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_sqrmod.o mp_sqrt.o \ -mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o \ -mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \ -s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \ -s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \ -s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \ -s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \ -s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \ -s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ -s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o +mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_small_prime_sieve_clear.o \ +mp_small_prime_sieve_init.o mp_small_prime_sieve_is_small_prime.o mp_small_prime_sieve_next_prime.o \ +mp_small_prime_sieve_prec_prime.o mp_sqrmod.o mp_sqrt.o mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o \ +mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o \ +s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o \ +s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o \ +s_mp_invmod_odd.o s_mp_log_2expt.o s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o \ +s_mp_mul_comba.o s_mp_mul_high.o s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o \ +s_mp_prime_is_divisible.o s_mp_prime_tab.o s_mp_radix_map.o s_mp_radix_size_overestimate.o \ +s_mp_rand_platform.o s_mp_rand_source.o s_mp_small_prime_sieve.o s_mp_sqr.o s_mp_sqr_comba.o \ +s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o s_mp_warray_put.o \ +s_mp_zero_buf.o s_mp_zero_digs.o HEADERS_PUB=tommath.h diff --git a/mp_small_prime_sieve_clear.c b/mp_small_prime_sieve_clear.c new file mode 100644 index 00000000..89765f38 --- /dev/null +++ b/mp_small_prime_sieve_clear.c @@ -0,0 +1,24 @@ +#include "tommath_private.h" +#ifdef MP_SMALL_PRIME_SIEVE_CLEAR_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +static void s_mp_erat_clear_one(mp_erat_single_sieve *sieve) +{ + if (sieve->content != NULL) { + MP_FREE(sieve->content, sieve->alloc); + sieve->alloc = 0; + } + sieve->size = 0; +} + +void mp_small_prime_sieve_clear(mp_erat_sieve *sieve) +{ + s_mp_erat_clear_one(&(sieve->base)); + s_mp_erat_clear_one(&(sieve->segment)); +} + + +#endif diff --git a/mp_small_prime_sieve_init.c b/mp_small_prime_sieve_init.c new file mode 100644 index 00000000..c0e5789d --- /dev/null +++ b/mp_small_prime_sieve_init.c @@ -0,0 +1,25 @@ +#include "tommath_private.h" +#ifdef MP_SMALL_PRIME_SIEVE_INIT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +mp_err mp_small_prime_sieve_init(mp_erat_sieve *sieve, bool warmup) +{ + mp_err err = MP_OKAY; + sieve->base.content = NULL; + sieve->base.alloc = 0; + sieve->base.size = 0; + sieve->segment.content = NULL; + sieve->segment.alloc = 0; + sieve->segment.size = 0; + sieve->single_segment_a = 0; + if (warmup) { + err = s_mp_erat_eratosthenes_init(ERAT_UINT_MAX_SQRT, &(sieve->base)); + } + return err; +} + + +#endif diff --git a/mp_small_prime_sieve_is_small_prime.c b/mp_small_prime_sieve_is_small_prime.c new file mode 100644 index 00000000..1615b73d --- /dev/null +++ b/mp_small_prime_sieve_is_small_prime.c @@ -0,0 +1,109 @@ +#include "tommath_private.h" +#ifdef MP_SMALL_PRIME_SIEVE_IS_SMALL_PRIME_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* TODO: change name to "mp_small_prime_sieve_is_prime" ("small" is already in it) + leave the "sieve" part out completely in all of the public sieve functions? + Yes, programmers and naming things is a mess. + + "It was, it is, and will always be, my dear firstborn son: Son_0" +*/ +/* + * Sets "result" to true if n is prime or false respectively. + * Also sets "result" to zero in case of error. + * Worst case runtime is: building a base sieve O(n log log n) and a segment + * O(n log log n) and search the segment linearly O(n). + */ +mp_err mp_small_prime_sieve_is_small_prime(ERAT_UINT n, bool *result, mp_erat_sieve *sieve) +{ + mp_err err = MP_OKAY; + ERAT_UINT a = 0, b = 0; + + if ((n < 2) || (n > ERAT_BIGGEST_PRIME)) { + *result = false; + return err; + } + + if (n == 2) { + *result = true; + return err; + } + + if ((n & 1) == 0) { + *result = false; + return err; + } + /* n is not 2 and odd, leaves 3,4,5,6,7 both composites are even, rest is prime. + Found no difference in runtime, maybe even slightly slower. + But all the large tests are sequential, so it maybe useful for random very small input? + */ + if (n <= 7) { + *result = 1; + return err; + } + + if (sieve->base.content == NULL) { + if ((err = s_mp_erat_eratosthenes_init(ERAT_UINT_MAX_SQRT, &(sieve->base))) != MP_OKAY) { + *result = false; + return err; + } + } + + /* No need to generate a segment if n is in the base sieve */ + if (n < (ERAT_UINT_MAX_SQRT)) { + /* might have been a small sieve, so check size of sieve first */ + if (n < sieve->base.size) { + *result = (s_mp_mp_erat_sieve_get_bit(&(sieve->base), (n - 1) / 2) == 1)?true:false; + return err; + } + } + /* no further shortcuts to apply, build and search a segment */ + + /* we have a segment and may be able to use it */ + if (sieve->segment.content != NULL) { + a = sieve->single_segment_a; + /* last segment may not fit into range_a_b */ + if (a > (ERAT_BIGGEST_PRIME - ERAT_UINT_MAX_SQRT)) { + b = ERAT_BIGGEST_PRIME; + } else { + b = a + ERAT_UINT_MAX_SQRT; + } + /* check if n is inside the bounds of the segment */ + if (n >= a && n <= b) { + *result = (s_mp_mp_erat_sieve_get_bit(&(sieve->segment), n - a))?true:false; + return err; + } + } + + /* + * A bit of heuristics ( "heuristics" is a more pretentious word for the + * commonly known expression "wild guess") + * Based on the vague idea of the assumption that most sieves get used for + * sequential series of primes or a single test here and there, but not + * for massive amounts of random requests. + */ + if (n > a) { + if (n > (ERAT_UINT_MAX - ERAT_UINT_MAX_SQRT)) { + a = ERAT_UINT_MAX - ERAT_UINT_MAX_SQRT; + } else { + a = n; + } + } else { + a = n - ERAT_UINT_MAX_SQRT; + } + if ((err = s_mp_erat_init_single_segment_with_start(a, + &(sieve->base), &(sieve->segment), &(sieve->single_segment_a))) != MP_OKAY) { + *result = 0; + return err; + } + /* finally, check for primality */ + *result = (s_mp_mp_erat_sieve_get_bit(&(sieve->segment), n - a))?true:false; + return err; +} + + + +#endif diff --git a/mp_small_prime_sieve_next_prime.c b/mp_small_prime_sieve_next_prime.c new file mode 100644 index 00000000..fb93983d --- /dev/null +++ b/mp_small_prime_sieve_next_prime.c @@ -0,0 +1,58 @@ +#include "tommath_private.h" +#ifdef MP_SMALL_PRIME_SIEVE_NEXT_PRIME_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + + +mp_err mp_small_prime_sieve_next_prime(ERAT_UINT n, ERAT_UINT *result, mp_erat_sieve *sieve) +{ + bool ret = false; + mp_err err = MP_OKAY; + + /* Shortcut */ + if (n < 2) { + *result = 2; + return err; + } + /* Even primes get a special treatment */ + if (n == 2) { + *result = 3; + return err; + } + /* Shortcut */ + if ((n % 2) == 0) { + n++; + } + + /* Check for upper limit */ + if (n == ERAT_BIGGEST_PRIME) { + *result = ERAT_BIGGEST_PRIME; + return err; + } + + /* Shortcut */ + if (n > ERAT_BIGGEST_PRIME) { + *result = 0; + return MP_OVF; + } + + /* Now we are down to the range 3 < n < ERAT_BIGGEST_PRIME with n odd */ + + /* As long as the flag is false check the odd numbers >= n */ + for (; ret == false; n+=2) { + /* just call is_small_prime(), it does all of the heavy work */ + if ((err = mp_small_prime_sieve_is_small_prime(n, &ret, sieve)) != MP_OKAY) { + *result = 0; + return err; + } + } + /* The index overshots be two, subtract two unconditionally */ + *result = n - 2; + return err; +} + + + +#endif diff --git a/mp_small_prime_sieve_prec_prime.c b/mp_small_prime_sieve_prec_prime.c new file mode 100644 index 00000000..bd82d53b --- /dev/null +++ b/mp_small_prime_sieve_prec_prime.c @@ -0,0 +1,48 @@ +#include "tommath_private.h" +#ifdef MP_SMALL_PRIME_SIEVE_PREC_PRIME_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +mp_err mp_small_prime_sieve_prec_prime(ERAT_UINT n, ERAT_UINT *result, mp_erat_sieve *sieve) +{ + bool ret = false; + mp_err err = MP_OKAY; + + /* Shortcut */ + if (n < 2) { + *result = 0; + return err; + } + + /* Even primes get a special treatment */ + if (n == 2) { + *result = 2; + return err; + } + + /* The rest of the primes are odd */ + if ((n % 2) == 0) { + n--; + } + + /* Now we are down to the range 2 < n < ERAT_UINT_MAX */ + + /* As long as the flag is false check the odd numbers <= n */ + for (; ret == false; n -= 2) { + if ((err = mp_small_prime_sieve_is_small_prime(n, &ret, sieve)) != MP_OKAY) { + *result = 0; + return err; + } + } + + /* The index overshots be two, add two unconditionally */ + *result = n + 2; + + return err; +} + + + +#endif diff --git a/s_mp_small_prime_sieve.c b/s_mp_small_prime_sieve.c new file mode 100644 index 00000000..78a2821c --- /dev/null +++ b/s_mp_small_prime_sieve.c @@ -0,0 +1,205 @@ +#include "tommath_private.h" +#ifdef S_MP_SMALL_PRIME_SIEVE_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + +#ifdef MP_USE_MEMOPS +#include +#endif + +ERAT_UINT s_mp_erat_isqrt(ERAT_UINT n) +{ + ERAT_UINT s, rem, root; + + if (n < 1) { + return 0; + } + /* highest power of four <= n */ + s = ((ERAT_UINT) 1) << (BITS_IN_ERAT_UINT - 2); + + rem = n; + root = 0; + while (s > n) { + s >>= 2; + } + while (s != 0) { + if (rem >= (s | root)) { + rem -= (s | root); + root >>= 1; + root |= s; + } else { + root >>= 1; + } + s >>= 2; + } + return root; +} +void s_mp_mp_erat_sieve_setall(mp_erat_single_sieve *bst) +{ +#ifdef MP_USE_MEMOPS + memset(bst->content, 0xFF, bst->alloc); +#else + size_t i, bs_size; + bs_size = bst->alloc / sizeof(ERAT_UINT); + for (i = 0; i < bs_size; i++) { + bst->content[i] = ERAT_UINT_MAX; + } +#endif +} + + +/* TODO: inline? macros? */ + +/* Set bit at position n to zero */ +void s_mp_mp_mp_erat_sieve_clear_bit(mp_erat_single_sieve *bst, ERAT_UINT n) +{ + ((*((bst)->content+(n/BITS_IN_ERAT_UINT)) &= ~(1lu<<(n % BITS_IN_ERAT_UINT)))); +} + +/* Read bit at position n */ +ERAT_UINT s_mp_mp_erat_sieve_get_bit(mp_erat_single_sieve *bst, ERAT_UINT n) +{ + return (((*((bst)->content+(n/BITS_IN_ERAT_UINT)) & (1lu<<(n % BITS_IN_ERAT_UINT))) != 0)); +} + +/* Seek set bit after position n */ +ERAT_UINT s_mp_mp_erat_sieve_nextset(mp_erat_single_sieve *bst, ERAT_UINT n) +{ + while ((n < (ERAT_UINT)(((bst)->size)) && (!s_mp_mp_erat_sieve_get_bit(bst, n)))) { + n++; + } + return n; +} + +/* + * Initiate a sieve that stores the odd numbers only: + * allocate memory, set actual size and allocated size and fill it + */ +mp_err s_mp_erat_eratosthenes_init(ERAT_UINT n, mp_erat_single_sieve *bst) +{ + bst->size = n; + /* The number BITS we need */ + n = (n + 1) / 2; + /* The number of BYTES we need is an eighth of the above. */ + n = (n + 1) / 8; + + bst->content = (ERAT_UINT *)MP_MALLOC(n + sizeof(ERAT_UINT)); + if (bst->content == NULL) { + return MP_MEM; + } + /* Take note */ + bst->alloc = n + sizeof(ERAT_UINT); + /* Fill sieve */ + s_mp_erat_eratosthenes(bst); + return MP_OKAY; +} + +/* + * Build a segment sieve with the largest reasonable size. "a" is the start of + * the sieve, Size is MIN(range_a_b,ERAT_UINT_MAX-a) + */ +mp_err s_mp_erat_init_single_segment_with_start( + ERAT_UINT a, + mp_erat_single_sieve *base_sieve, + mp_erat_single_sieve *single_segment, + ERAT_UINT *single_segment_a) +{ + ERAT_UINT b, n; + mp_err e = MP_OKAY; + + if (a > (ERAT_BIGGEST_PRIME - ERAT_UINT_MAX_SQRT)) { + b = ERAT_BIGGEST_PRIME; + } else { + b = a + ERAT_UINT_MAX_SQRT; + } + + n = ((b - a)+1)/CHAR_BIT; + + /* + We are reusing that slice of heap for all segments. + Support of multithreading would need a bit of work. + */ + if (single_segment->content == NULL) { + single_segment->content = (ERAT_UINT *) MP_MALLOC(n + sizeof(ERAT_UINT)); + if (single_segment->content == NULL) { + return MP_MEM; + } + single_segment->alloc = n + sizeof(ERAT_UINT); + single_segment->size = n * CHAR_BIT; + } + + s_mp_erat_eratosthenes_segment(a, b, base_sieve, single_segment); + *single_segment_a = a; + + return e; +} + + +/* + * Simple Eratosthenes' sieve, starting at zero + * Keeping odd numbers only as the single optimization + */ +void s_mp_erat_eratosthenes(mp_erat_single_sieve *bst) +{ + ERAT_UINT n, k, r, j; + + n = (ERAT_UINT)bst->size; + r = s_mp_erat_isqrt(n); + s_mp_mp_erat_sieve_setall(bst); + + for (k = 1; k < ((r - 1) / 2); k += 1) { + if (s_mp_mp_erat_sieve_get_bit(bst, k)) { + for (j = k * 2 * (k + 1); j < (n - 1) / 2; j += 2 * k + 1) { + s_mp_mp_mp_erat_sieve_clear_bit(bst, j); + } + } + } +} + +ERAT_UINT s_mp_erat_nextprime(ERAT_UINT p, mp_erat_single_sieve *bst) +{ + ERAT_UINT ret; + if (p <= 1) { + return 2; + } + ret = s_mp_mp_erat_sieve_nextset(bst, ((p - 1) / 2) + 1); + return 2 * ret + 1; +} + +/* + * TODO: Not memory optimized, it stores the even numbers, too. + * Should be ok if the segments are small but needs to + * be done at some time in the near future. + * + * Fill sieve "segment" of the range [a,b] from the basic sieve "base" + * Default size of a segment (32-bit): 8192 bytes (plus some angst-allowance) + */ +void s_mp_erat_eratosthenes_segment(ERAT_UINT a, ERAT_UINT b, mp_erat_single_sieve *base, mp_erat_single_sieve *segment) +{ + ERAT_UINT r, j; + ERAT_UINT p; + + r = s_mp_erat_isqrt(b); + + s_mp_mp_erat_sieve_setall(segment); + p = 0; + + while (p <= r) { + p = s_mp_erat_nextprime(p, base); + j = p * p; + if (j < a) { + j = ((a + p - 1) / p) * p; + } + for (; j <= b; j += p) { + /* j+=p can overflow, so check size of j in relation to a */ + if (j >= a) { + s_mp_mp_mp_erat_sieve_clear_bit(segment, j - a); + } else { + break; + } + } + } +} + + +#endif diff --git a/sources.cmake b/sources.cmake index ba3012ef..dc65e424 100644 --- a/sources.cmake +++ b/sources.cmake @@ -111,6 +111,11 @@ mp_set_u64.c mp_set_ul.c mp_shrink.c mp_signed_rsh.c +mp_small_prime_sieve_clear.c +mp_small_prime_sieve_init.c +mp_small_prime_sieve_is_small_prime.c +mp_small_prime_sieve_next_prime.c +mp_small_prime_sieve_prec_prime.c mp_sqrmod.c mp_sqrt.c mp_sqrtmod_prime.c @@ -153,6 +158,7 @@ s_mp_radix_map.c s_mp_radix_size_overestimate.c s_mp_rand_platform.c s_mp_rand_source.c +s_mp_small_prime_sieve.c s_mp_sqr.c s_mp_sqr_comba.c s_mp_sqr_karatsuba.c diff --git a/tommath.def b/tommath.def index ed5aa8b0..a336b43f 100644 --- a/tommath.def +++ b/tommath.def @@ -114,6 +114,11 @@ EXPORTS mp_set_ul mp_shrink mp_signed_rsh + mp_small_prime_sieve_clear + mp_small_prime_sieve_init + mp_small_prime_sieve_is_small_prime + mp_small_prime_sieve_next_prime + mp_small_prime_sieve_prec_prime mp_sqrmod mp_sqrt mp_sqrtmod_prime diff --git a/tommath.h b/tommath.h index 1820d243..d240797e 100644 --- a/tommath.h +++ b/tommath.h @@ -566,6 +566,48 @@ mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style) MP_WUR; */ mp_err mp_prime_rand(mp_int *a, int t, int size, int flags) MP_WUR; + +/* A small segmented sieve */ + +/* Results in a sieve-size of 4,100 bytes (base) + 8,196 bytes (segment) */ +#ifndef ERAT_BIGGEST_PRIME +#define ERAT_BIGGEST_PRIME 0xFFFFFFFBlu +#define ERAT_BIGGEST_BASE_PRIME 0xFFF1lu +#define ERAT_UINT uint32_t +#define ERAT_UINT_MAX 0xFFFFFFFFlu +#define ERAT_UINT_MAX_SQRT 0xFFFFlu + +#define BITS_IN_ERAT_UINT (sizeof(ERAT_UINT)*CHAR_BIT) +#endif + +typedef struct { + ERAT_UINT *content; /* bitset holding the sieve */ + size_t size; /* number of entries (which is a slightly misleading description) */ + size_t alloc; /* size in bytes */ +} mp_erat_single_sieve; + +typedef struct { + mp_erat_single_sieve base; /* base sieve (0 -> ERAT_UINT_MAX_SQRT) */ + mp_erat_single_sieve segment; /* segment (range_a_b) */ + ERAT_UINT single_segment_a; /* startpoint of segment */ +} mp_erat_sieve; + + +/* Check primality of n if n < 2^32 */ +mp_err mp_small_prime_sieve_is_small_prime(ERAT_UINT n, bool *result, mp_erat_sieve *sieve) MP_WUR; +/* Returns next prime p if the next prime is p < 2^32 */ +mp_err mp_small_prime_sieve_next_prime(ERAT_UINT n, ERAT_UINT *result, mp_erat_sieve *sieve) MP_WUR; +/* Returns preceding prime p if the preceding prime is p > 1 */ +mp_err mp_small_prime_sieve_prec_prime(ERAT_UINT n, ERAT_UINT *result, mp_erat_sieve *sieve) MP_WUR; +/* Reset sieve and clear memory */ +void mp_small_prime_sieve_clear(mp_erat_sieve *sieve); +/* Init sieve to default values, build base-sieve, too, if the parameter warmup is true + otherwise the first call to either mp_small_prime_sieve_is_small_prime, + mp_small_prime_sieve_next_prime, or mp_small_prime_sieve_prec_prime do it at first call */ +mp_err mp_small_prime_sieve_init(mp_erat_sieve *sieve, bool warmup) MP_WUR; + + + /* ---> radix conversion <--- */ int mp_count_bits(const mp_int *a) MP_WUR; diff --git a/tommath_class.h b/tommath_class.h index 40ddfd4b..6e441713 100644 --- a/tommath_class.h +++ b/tommath_class.h @@ -120,6 +120,11 @@ # define MP_SET_UL_C # define MP_SHRINK_C # define MP_SIGNED_RSH_C +# define MP_SMALL_PRIME_SIEVE_CLEAR_C +# define MP_SMALL_PRIME_SIEVE_INIT_C +# define MP_SMALL_PRIME_SIEVE_IS_SMALL_PRIME_C +# define MP_SMALL_PRIME_SIEVE_NEXT_PRIME_C +# define MP_SMALL_PRIME_SIEVE_PREC_PRIME_C # define MP_SQRMOD_C # define MP_SQRT_C # define MP_SQRTMOD_PRIME_C @@ -162,6 +167,7 @@ # define S_MP_RADIX_SIZE_OVERESTIMATE_C # define S_MP_RAND_PLATFORM_C # define S_MP_RAND_SOURCE_C +# define S_MP_SMALL_PRIME_SIEVE_C # define S_MP_SQR_C # define S_MP_SQR_COMBA_C # define S_MP_SQR_KARATSUBA_C @@ -882,6 +888,28 @@ # define MP_SUB_D_C #endif +#if defined(MP_SMALL_PRIME_SIEVE_CLEAR_C) +# define S_MP_ERAT_CLEAR_ONE_C +#endif + +#if defined(MP_SMALL_PRIME_SIEVE_INIT_C) +# define S_MP_ERAT_ERATOSTHENES_INIT_C +#endif + +#if defined(MP_SMALL_PRIME_SIEVE_IS_SMALL_PRIME_C) +# define S_MP_ERAT_ERATOSTHENES_INIT_C +# define S_MP_ERAT_INIT_SINGLE_SEGMENT_WITH_START_C +# define S_MP_MP_ERAT_SIEVE_GET_BIT_C +#endif + +#if defined(MP_SMALL_PRIME_SIEVE_NEXT_PRIME_C) +# define MP_SMALL_PRIME_SIEVE_IS_SMALL_PRIME_C +#endif + +#if defined(MP_SMALL_PRIME_SIEVE_PREC_PRIME_C) +# define MP_SMALL_PRIME_SIEVE_IS_SMALL_PRIME_C +#endif + #if defined(MP_SQRMOD_C) # define MP_MOD_C # define MP_MUL_C @@ -1253,6 +1281,19 @@ # define S_MP_RAND_PLATFORM_C #endif +#if defined(S_MP_SMALL_PRIME_SIEVE_C) +# define S_MP_ERAT_ERATOSTHENES_C +# define S_MP_ERAT_ERATOSTHENES_INIT_C +# define S_MP_ERAT_ERATOSTHENES_SEGMENT_C +# define S_MP_ERAT_INIT_SINGLE_SEGMENT_WITH_START_C +# define S_MP_ERAT_ISQRT_C +# define S_MP_ERAT_NEXTPRIME_C +# define S_MP_MP_ERAT_SIEVE_GET_BIT_C +# define S_MP_MP_ERAT_SIEVE_NEXTSET_C +# define S_MP_MP_ERAT_SIEVE_SETALL_C +# define S_MP_MP_MP_ERAT_SIEVE_CLEAR_BIT_C +#endif + #if defined(S_MP_SQR_C) # define MP_CLAMP_C # define MP_CLEAR_C diff --git a/tommath_private.h b/tommath_private.h index be620dbc..60d387ff 100644 --- a/tommath_private.h +++ b/tommath_private.h @@ -234,6 +234,26 @@ MP_PRIVATE mp_err s_mp_radix_size_overestimate(const mp_int *a, const int radix, MP_PRIVATE mp_err s_mp_fp_log(const mp_int *a, mp_int *c) MP_WUR; MP_PRIVATE mp_err s_mp_fp_log_d(const mp_int *a, mp_word *c) MP_WUR; + + +MP_PRIVATE ERAT_UINT s_mp_erat_nextprime(ERAT_UINT p, mp_erat_single_sieve *bst); +MP_PRIVATE void s_mp_erat_eratosthenes(mp_erat_single_sieve *bst); +MP_PRIVATE mp_err s_mp_erat_eratosthenes_init(ERAT_UINT n, mp_erat_single_sieve *bst) MP_WUR; +MP_PRIVATE void s_mp_erat_eratosthenes_segment(ERAT_UINT a, ERAT_UINT b, mp_erat_single_sieve *base, + mp_erat_single_sieve *segment); +MP_PRIVATE void s_mp_erat_eratosthenes_segment_clear(mp_erat_single_sieve *segment, ERAT_UINT *single_segment_a); +MP_PRIVATE mp_err s_mp_erat_init_single_segment_with_start(ERAT_UINT a, mp_erat_single_sieve *base_sieve, + mp_erat_single_sieve *single_segment, ERAT_UINT *single_segment_a) MP_WUR; +MP_PRIVATE ERAT_UINT s_mp_erat_isqrt(ERAT_UINT n); +MP_PRIVATE ERAT_UINT s_mp_mp_erat_sieve_get_bit(mp_erat_single_sieve *bst, ERAT_UINT n); +MP_PRIVATE ERAT_UINT s_mp_mp_erat_sieve_nextset(mp_erat_single_sieve *bst, ERAT_UINT n); +MP_PRIVATE void s_mp_mp_erat_sieve_setall(mp_erat_single_sieve *bst); +MP_PRIVATE void s_mp_mp_mp_erat_sieve_clear_bit(mp_erat_single_sieve *bst, ERAT_UINT n); + + + + + #ifdef MP_SMALL_STACK_SIZE #if defined(__GNUC__)