Skip to content

A prime sieve #197

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

A prime sieve #197

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
178 changes: 177 additions & 1 deletion demo/test.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -811,7 +811,7 @@ static int test_mp_prime_rand(void)

/* test for size */
for (ix = 10; ix < 128; ix++) {
printf("Testing (not safe-prime): %9d bits \n", ix);
printf("\rTesting (not safe-prime): %9d bits ", ix);
fflush(stdout);
DO(mp_prime_rand(&a, 8, ix, (rand_int() & 1) ? 0 : MP_PRIME_2MSB_ON));
EXPECT(mp_count_bits(&a) == ix);
Expand Down Expand Up @@ -1529,6 +1529,177 @@ static int test_mp_decr(void)
return EXIT_FAILURE;
}

static int test_mp_is_small_prime(void)
{
mp_sieve sieve;
mp_err e;
int i, test_size;

const mp_sieve_prime to_test[] = {
52, 137, 153, 179, 6, 153, 53, 132, 150, 65,
27414, 36339, 36155, 11067, 52060, 5741,
29755, 2698, 52572, 13053, 9375, 47241,39626,
207423, 128857, 37419, 141696, 189465,
41503, 127370, 91673, 8473, 479142414, 465566339,
961126169, 1057886067, 1222702060, 1017450741,
1019879755, 72282698, 2048787577, 2058368053
};
const bool tested[] = {
false, true, false, true, false, false, true, false, false,
false, false, false, false, false, false, true, false, false,
false, false, false, false, false, false, true, false, false,
false, false, false, true, false, false, false, true, false,
false, false, false, false, true, false
};
bool result;

mp_sieve_init(&sieve);

test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime));

for (i = 0; i < test_size; i++) {
if ((e = mp_is_small_prime(to_test[i], &result, &sieve)) != MP_OKAY) {
fprintf(stderr,"mp_is_small_prime failed: \"%s\"\n",mp_error_to_string(e));
goto LTM_ERR;
}
if (result != tested[i]) {
fprintf(stderr,"mp_is_small_prime failed for %lu. Said %lu but is %lu \n",
(unsigned long)to_test[i], (unsigned long)result, (unsigned long)tested[i]);
goto LTM_ERR;
}
}
mp_sieve_clear(&sieve);
return EXIT_SUCCESS;
LTM_ERR:
mp_sieve_clear(&sieve);
return EXIT_FAILURE;
}

static int test_mp_next_small_prime(void)
{
mp_sieve sieve;
mp_sieve_prime ret = 0lu, p;
mp_int primesum, t;
mp_err e;
int i, test_size;

/* Jumping wildly to and fro */
const mp_sieve_prime to_test[] = {
52, 137, 153, 179, 6, 153, 53, 132, 150, 65,
27414, 36339, 36155, 11067, 52060, 5741,
29755, 2698, 52572, 13053, 9375, 47241,
39626, 207423, 128857, 37419, 141696, 189465,
41503, 127370, 91673, 8473, 479142414, 465566339,
961126169, 1057886067, 1222702060, 1017450741,
1019879755, 72282698, 2048787577, 2058368053
};
const mp_sieve_prime tested[] = {
53, 137, 157, 179, 7, 157, 53, 137, 151, 67,
27427, 36341, 36161, 11069, 52067, 5741,
29759, 2699, 52579, 13063, 9377, 47251,
39631, 207433, 128857, 37423, 141697, 189467,
41507, 127373, 91673, 8501, 479142427, 465566393,
961126169, 1057886083, 1222702081, 1017450823,
1019879761, 72282701, 2048787577, 2058368113
};
const char *primesum_32 = "202259606268580";

mp_sieve_init(&sieve);

test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime));

for (i = 0; i < test_size; i++) {
if ((e = mp_next_small_prime(to_test[i], &ret, &sieve)) != MP_OKAY) {
fprintf(stderr,"mp_next_small_prime failed with \"%s\" at index %d\n",
mp_error_to_string(e), i);
goto LBL_ERR;
}
if (ret != tested[i]) {
fprintf(stderr,"mp_next_small_prime failed for %lu. Said %lu but is %lu \n",
(unsigned long)to_test[i], (unsigned long)ret, (unsigned long)tested[i]);
goto LBL_ERR;
}
}

DOR(mp_init_multi(&primesum, &t, NULL));

for (p = 4293918720lu; ret < (mp_sieve_prime)MP_SIEVE_BIGGEST_PRIME;) {
DO(mp_next_small_prime(p, &ret, &sieve));
p = ret + 1u;
#ifdef MP_64BIT
mp_set_u64(&t, ret);
#else
mp_set_u32(&t, ret);
#endif
DO(mp_add(&primesum, &t, &primesum));
}
printf("Primesum computed: ");
DO(mp_fwrite(&primesum, 10, stdout));
puts("");
DO(mp_read_radix(&t, primesum_32, 10));
EXPECT(mp_cmp(&primesum, &t) == MP_EQ);

mp_sieve_clear(&sieve);
mp_clear_multi(&primesum, &t, NULL);
return EXIT_SUCCESS;
LBL_ERR:
mp_clear_multi(&primesum, &t, NULL);
mp_sieve_clear(&sieve);
return EXIT_FAILURE;
}

static int test_mp_prec_small_prime(void)
{
mp_sieve sieve;
mp_sieve_prime ret;
mp_err e;
int i, test_size;

const mp_sieve_prime to_test[] = {
52, 137, 153, 179, 6, 153, 53, 132, 150, 65,
27414, 36339, 36155, 11067, 52060, 5741,
29755, 2698, 52572, 13053, 9375, 47241,
39626, 207423, 128857, 37419, 141696, 189465,
41503, 127370, 91673, 8473, 479142414, 465566339,
961126169, 1057886067, 1222702060, 1017450741,
1019879755, 72282698, 2048787577, 2058368053
};
const mp_sieve_prime tested[] = {
47, 137, 151, 179, 5, 151, 53, 131, 149, 61,
27409, 36319, 36151, 11059, 52057, 5741,
29753, 2693, 52571, 13049, 9371, 47237,
39623, 207409, 128857, 37409, 141689, 189463,
41491, 127363, 91673, 8467, 479142413, 465566323,
961126169, 1057886029, 1222702051, 1017450739,
1019879717, 72282697, 2048787577, 2058368051
};

mp_sieve_init(&sieve);

test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime));

for (i = 0; i < test_size; i++) {
if ((e = mp_prec_small_prime(to_test[i], &ret, &sieve)) != MP_OKAY) {
fprintf(stderr,"mp_prec_small_prime failed with \"%s\" at index %d\n",
mp_error_to_string(e), i);
goto LTM_ERR;
}
if (ret != tested[i]) {
fprintf(stderr,"mp_prec_small_prime failed for %lu. Said %lu but is %lu \n",
(unsigned long)to_test[i], (unsigned long)ret, (unsigned long)tested[i]);
goto LTM_ERR;
}
}

mp_sieve_clear(&sieve);
return EXIT_SUCCESS;
LTM_ERR:
mp_sieve_clear(&sieve);
return EXIT_FAILURE;
}



/*
Cannot test mp_exp(_d) without mp_root_n and vice versa.
So one of the two has to be tested from scratch.
Expand Down Expand Up @@ -2240,6 +2411,11 @@ static int unit_tests(int argc, char **argv)
T1(mp_prime_next_prime, MP_PRIME_NEXT_PRIME),
T1(mp_prime_rand, MP_PRIME_RAND),
T1(mp_rand, MP_RAND),

T1(mp_is_small_prime, MP_IS_SMALL_PRIME),
T1(mp_next_small_prime, MP_NEXT_SMALL_PRIME),
T1(mp_prec_small_prime, MP_PREC_SMALL_PRIME),

T1(mp_read_radix, MP_READ_RADIX),
T1(mp_read_write_ubin, MP_TO_UBIN),
T1(mp_read_write_sbin, MP_TO_SBIN),
Expand Down
170 changes: 170 additions & 0 deletions doc/bn.tex
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2347,6 +2347,173 @@ \section{Random Primes}
\label{fig:primeopts}
\end{figure}


\section{Prime Sieve}
The prime sieve is implemented as a simple segmented Sieve of Eratosthenes.

This library needs some small sequential amounts for divisibility tests starting at two and
some random small primes for the Miller--Rabin test. That means we need a small base sieve
for quick results with a cold start and small segments to get random small primes fast.

The base sieve holds odd values only and even with a size of \texttt{4096} bytes it is
small enough to get build quickly, in about 50 $\mu$sec on the author's old machine.

The choice of the size for the individual segments varies more in the results. See table
\ref{table:sievebenchmarks} for some data. Size must be of the form $-1 + 2^n$.
Default size is \texttt{0xFFF}. It might be of interest that the biggest prime gap
below $2^{32}$ is $335$.

\begin{table}[h]
\begin{center}
\begin{tabular}{r c l}
\textbf{Segment Size (bits)} & \textbf{Random Access in $\mu$sec} & \textbf{Full Primesum} \\
\texttt{ 0x1F} & (60) & - \\
\texttt{ 0x3F} & (75) & - \\
\texttt{ 0x7F} & 63 & - \\
\texttt{ 0xFF} & 70 & 26 min \\
\texttt{ 0x1FF} & 73 & 13 min \\
\texttt{ 0x3FF} & 75 & 7 min 17 sec \\
\texttt{ 0x7FF} & 85 & 4 min 10 sec \\
\texttt{ 0xFFF} & 100 & 2 min 52 sec \\
\texttt{ 0x1FFF} & 120 & 1 min 45 sec \\
\texttt{ 0x3FFF} & 145 & 1 min 19 sec \\
\texttt{ 0x7FFF} & 200 & 1 min 04 sec \\
\texttt{ 0xFFFF} & 350 & 0 min 57 sec \\
\texttt{ 0xFFFFFFFF} & 300 & 0 min 35 sec
\end{tabular}
\caption{Average access times (warm) with the default base sieve and \texttt{MP-64BIT}} \label{table:sievebenchmarks}
\end{center}
\end{table}

The default sizes are in \texttt{tommath\_private.h}: \texttt{MP\_SIEVE\_PRIME\_MAX\_SQRT} is
the size of the base sieve and \texttt{MP\_SIEVE\_RANGE\_A\_B} is the size of the segment.

\subsection{Initialization and Clearing}
Initializing. It cannot fail because it only sets some default values. Memory is allocated later according to needs.
\index{mp\_sieve\_init}
\begin{alltt}
void mp_sieve_init(mp_sieve *sieve);
\end{alltt}
The function \texttt{mp\_sieve\_init} is equivalent to
\begin{alltt}
mp_sieve sieve = {NULL, NULL, 0};
\end{alltt}

Free the memory used by the sieve with
\index{mp\_sieve\_clear}
\begin{alltt}
void mp_sieve_clear(mp_sieve *sieve);
\end{alltt}

\subsection{Primality Test of Small Numbers}
Individual small numbers can be tested for primality with
\index{mp\_is\_small\_prime}
\begin{alltt}
int mp_is_small_prime(mp_sieve_prime n, bool *result,
mp_sieve *sieve);
\end{alltt}
It sets the variable \texttt{result} to \texttt{false} if the given number in \texttt{n} is a composite and
\texttt{true} if it is a prime.

It will return \texttt{MP\_SIEVE\_MAX\_REACHED} to flag the content of \texttt{result} as the last valid one.

The implementation of this function does all the heavy lifting--the building of the base sieve and the segment
if one is necessary. The base sieve stays, so this function can be used to ``warm up'' the sieve and, if
\texttt{n} is slightly larger than the upper limit of the base sieve, ``warm up'' the first segment, too.

\subsection{Find Adjacent Primes}
To find the prime bigger than a number \texttt{n} use
\index{mp\_next\_small\_prime}
\begin{alltt}
int mp_next_small_prime(mp_sieve_prime n, mp_sieve_prime *result,
mp_sieve *sieve);
\end{alltt}
and to find the one smaller than \texttt{n}
\index{mp\_prec\_small\_prime}
\begin{alltt}
int mp_prec_small_prime(mp_sieve_prime n, mp_sieve_prime *result,
mp_sieve *sieve);
\end{alltt}
Both functions return \texttt{n} if \texttt{n} is prime. Use \texttt{n + 1} to get
the prime after \texttt{n} if \texttt{n} is prime and \texttt{n - 1} to get the
the prime preceding \texttt{n} if \texttt{n} is prime.

\subsection{Useful Constants}
\begin{description}
\item[\texttt{MP\_SIEVE\_BIGGEST\_PRIME}] \texttt{read-only} The biggest prime the sieve can offer.
It is $4\,294\,967\,291$ for all \texttt{MP\_xBIT}.

\item[\texttt{mp\_sieve\_prime}] \texttt{read-only} The basic type for the numbers in the sieve. It
is \texttt{uint32\_t} for \texttt{MP\_16BIT}, \texttt{MP\_32BIT}; and
\texttt{uint64\_t} for \texttt{MP\_64BIT}..

\item[\texttt{MP\_SIEVE\_PR\_UINT}] Choses the correct macro from \texttt{inttypes.h} to print a\\
\texttt{mp\_sieve\_prime}. The header \texttt{inttypes.h} must be included before\\
\texttt{tommath.h} for it to work.
\end{description}
\subsection{Examples}\label{sec:spnexamples}
\subsubsection{Initialization and Clearing}
Using a sieve follows the same procedure as using a bigint:
\begin{alltt}
/* Declaration */
mp_sieve sieve;

/*
Initialization.
Cannot fail, only sets a handful of default values.
Memory allocation is done in the library itself
according to needs.
*/
mp_sieve_init(&sieve);

/* use the sieve */

/* Clean up */
mp_sieve_clear(&sieve);
\end{alltt}
\subsubsection{Primality Test}
A small program to test the input of a small number for primality.
\begin{alltt}
#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
/*inttypes.h is needed for printing and must be included before tommath.h*/
#include <inttypes.h>
#include "tommath.h"
int main(int argc, char **argv)
{
mp_sieve_prime number;
mp_sieve sieve;
int e;
/* variable holding the result of mp_is_small_prime */
mp_sieve_prime result;

if (argc != 2) {
fprintf(stderr,"Usage %s number\textbackslash{}n", argv[0]);
exit(EXIT_FAILURE);
}
number = (mp_sieve_prime) strtoul(argv[1],NULL, 10);
if (errno == ERANGE) {
fprintf(stderr,"strtoul(l) failed: input out of range\textbackslash{}n");
goto LTM_ERR;
}
mp_sieve_init(&sieve);
if ((e = mp_is_small_prime(number, &result, &sieve)) != MP_OKAY) {
fprintf(stderr,"mp_is_small_prime failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
mp_error_to_string(e));
goto LTM_ERR;
}
printf("The number %" LTM_SIEVE_PR_UINT " is %s prime\textbackslash{}n",
number,(result)?"":"not");

mp_sieve_clear(&sieve);
exit(EXIT_SUCCESS);
LTM_ERR:
mp_sieve_clear(&sieve);
exit(EXIT_FAILURE);
}
\end{alltt}

\chapter{Random Number Generation}
\section{PRNG}
\index{mp\_rand}
Expand Down Expand Up @@ -2417,6 +2584,9 @@ \subsection{To ASCII}
mp_err mp_fwrite(const mp_int *a, int radix, FILE *stream);
\end{alltt}




\subsection{From ASCII}
\index{mp\_read\_radix}
\begin{alltt}
Expand Down
2 changes: 1 addition & 1 deletion helper.pl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_sieve_setall.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*\(/;
Expand Down
Loading