Addition of a small segmented sieve

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<n<2^32 */
+   for (i = 0; i < 100; i++) {
+      n = (uint32_t)rand_uint();
+      if ((n < 2) || (n > 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),

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}

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

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*\(/;

libtommath_VS2008.vcproj

@@ -748,6 +748,26 @@
 			RelativePath="mp_signed_rsh.c"
 			>
 		</File>
+		<File
+			RelativePath="mp_small_prime_sieve_clear.c"
+			>
+		</File>
+		<File
+			RelativePath="mp_small_prime_sieve_init.c"
+			>
+		</File>
+		<File
+			RelativePath="mp_small_prime_sieve_is_small_prime.c"
+			>
+		</File>
+		<File
+			RelativePath="mp_small_prime_sieve_next_prime.c"
+			>
+		</File>
+		<File
+			RelativePath="mp_small_prime_sieve_prec_prime.c"
+			>
+		</File>
 		<File
 			RelativePath="mp_sqrmod.c"
 			>
@@ -916,6 +936,10 @@
 			RelativePath="s_mp_rand_source.c"
 			>
 		</File>
+		<File
+			RelativePath="s_mp_small_prime_sieve.c"
+			>
+		</File>
 		<File
 			RelativePath="s_mp_sqr.c"
 			>

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
 

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)