Skip to content

Commit c3216ae

Browse files
czurniedenczurnieden
committed
A prime sieve
1 parent 96f9edf commit c3216ae

30 files changed

+1200
-84
lines changed

‎demo/test.c

+177-1
Original file line numberDiff line numberDiff line change
@@ -811,7 +811,7 @@ static int test_mp_prime_rand(void)
811811

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

1532+
static int test_mp_is_small_prime(void)
1533+
{
1534+
mp_sieve sieve;
1535+
mp_err e;
1536+
int i, test_size;
1537+
1538+
const mp_sieve_prime to_test[] = {
1539+
52, 137, 153, 179, 6, 153, 53, 132, 150, 65,
1540+
27414, 36339, 36155, 11067, 52060, 5741,
1541+
29755, 2698, 52572, 13053, 9375, 47241,39626,
1542+
207423, 128857, 37419, 141696, 189465,
1543+
41503, 127370, 91673, 8473, 479142414, 465566339,
1544+
961126169, 1057886067, 1222702060, 1017450741,
1545+
1019879755, 72282698, 2048787577, 2058368053
1546+
};
1547+
const bool tested[] = {
1548+
false, true, false, true, false, false, true, false, false,
1549+
false, false, false, false, false, false, true, false, false,
1550+
false, false, false, false, false, false, true, false, false,
1551+
false, false, false, true, false, false, false, true, false,
1552+
false, false, false, false, true, false
1553+
};
1554+
bool result;
1555+
1556+
mp_sieve_init(&sieve);
1557+
1558+
test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime));
1559+
1560+
for (i = 0; i < test_size; i++) {
1561+
if ((e = mp_is_small_prime(to_test[i], &result, &sieve)) != MP_OKAY) {
1562+
fprintf(stderr,"mp_is_small_prime failed: \"%s\"\n",mp_error_to_string(e));
1563+
goto LTM_ERR;
1564+
}
1565+
if (result != tested[i]) {
1566+
fprintf(stderr,"mp_is_small_prime failed for %lu. Said %lu but is %lu \n",
1567+
(unsigned long)to_test[i], (unsigned long)result, (unsigned long)tested[i]);
1568+
goto LTM_ERR;
1569+
}
1570+
}
1571+
mp_sieve_clear(&sieve);
1572+
return EXIT_SUCCESS;
1573+
LTM_ERR:
1574+
mp_sieve_clear(&sieve);
1575+
return EXIT_FAILURE;
1576+
}
1577+
1578+
static int test_mp_next_small_prime(void)
1579+
{
1580+
mp_sieve sieve;
1581+
mp_sieve_prime ret = 0lu, p;
1582+
mp_int primesum, t;
1583+
mp_err e;
1584+
int i, test_size;
1585+
1586+
/* Jumping wildly to and fro */
1587+
const mp_sieve_prime to_test[] = {
1588+
52, 137, 153, 179, 6, 153, 53, 132, 150, 65,
1589+
27414, 36339, 36155, 11067, 52060, 5741,
1590+
29755, 2698, 52572, 13053, 9375, 47241,
1591+
39626, 207423, 128857, 37419, 141696, 189465,
1592+
41503, 127370, 91673, 8473, 479142414, 465566339,
1593+
961126169, 1057886067, 1222702060, 1017450741,
1594+
1019879755, 72282698, 2048787577, 2058368053
1595+
};
1596+
const mp_sieve_prime tested[] = {
1597+
53, 137, 157, 179, 7, 157, 53, 137, 151, 67,
1598+
27427, 36341, 36161, 11069, 52067, 5741,
1599+
29759, 2699, 52579, 13063, 9377, 47251,
1600+
39631, 207433, 128857, 37423, 141697, 189467,
1601+
41507, 127373, 91673, 8501, 479142427, 465566393,
1602+
961126169, 1057886083, 1222702081, 1017450823,
1603+
1019879761, 72282701, 2048787577, 2058368113
1604+
};
1605+
const char *primesum_32 = "202259606268580";
1606+
1607+
mp_sieve_init(&sieve);
1608+
1609+
test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime));
1610+
1611+
for (i = 0; i < test_size; i++) {
1612+
if ((e = mp_next_small_prime(to_test[i], &ret, &sieve)) != MP_OKAY) {
1613+
fprintf(stderr,"mp_next_small_prime failed with \"%s\" at index %d\n",
1614+
mp_error_to_string(e), i);
1615+
goto LBL_ERR;
1616+
}
1617+
if (ret != tested[i]) {
1618+
fprintf(stderr,"mp_next_small_prime failed for %lu. Said %lu but is %lu \n",
1619+
(unsigned long)to_test[i], (unsigned long)ret, (unsigned long)tested[i]);
1620+
goto LBL_ERR;
1621+
}
1622+
}
1623+
1624+
DOR(mp_init_multi(&primesum, &t, NULL));
1625+
1626+
for (p = 4293918720lu; ret < (mp_sieve_prime)MP_SIEVE_BIGGEST_PRIME;) {
1627+
DO(mp_next_small_prime(p, &ret, &sieve));
1628+
p = ret + 1u;
1629+
#ifdef MP_64BIT
1630+
mp_set_u64(&t, ret);
1631+
#else
1632+
mp_set_u32(&t, ret);
1633+
#endif
1634+
DO(mp_add(&primesum, &t, &primesum));
1635+
}
1636+
printf("Primesum computed: ");
1637+
DO(mp_fwrite(&primesum, 10, stdout));
1638+
puts("");
1639+
DO(mp_read_radix(&t, primesum_32, 10));
1640+
EXPECT(mp_cmp(&primesum, &t) == MP_EQ);
1641+
1642+
mp_sieve_clear(&sieve);
1643+
mp_clear_multi(&primesum, &t, NULL);
1644+
return EXIT_SUCCESS;
1645+
LBL_ERR:
1646+
mp_clear_multi(&primesum, &t, NULL);
1647+
mp_sieve_clear(&sieve);
1648+
return EXIT_FAILURE;
1649+
}
1650+
1651+
static int test_mp_prec_small_prime(void)
1652+
{
1653+
mp_sieve sieve;
1654+
mp_sieve_prime ret;
1655+
mp_err e;
1656+
int i, test_size;
1657+
1658+
const mp_sieve_prime to_test[] = {
1659+
52, 137, 153, 179, 6, 153, 53, 132, 150, 65,
1660+
27414, 36339, 36155, 11067, 52060, 5741,
1661+
29755, 2698, 52572, 13053, 9375, 47241,
1662+
39626, 207423, 128857, 37419, 141696, 189465,
1663+
41503, 127370, 91673, 8473, 479142414, 465566339,
1664+
961126169, 1057886067, 1222702060, 1017450741,
1665+
1019879755, 72282698, 2048787577, 2058368053
1666+
};
1667+
const mp_sieve_prime tested[] = {
1668+
47, 137, 151, 179, 5, 151, 53, 131, 149, 61,
1669+
27409, 36319, 36151, 11059, 52057, 5741,
1670+
29753, 2693, 52571, 13049, 9371, 47237,
1671+
39623, 207409, 128857, 37409, 141689, 189463,
1672+
41491, 127363, 91673, 8467, 479142413, 465566323,
1673+
961126169, 1057886029, 1222702051, 1017450739,
1674+
1019879717, 72282697, 2048787577, 2058368051
1675+
};
1676+
1677+
mp_sieve_init(&sieve);
1678+
1679+
test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime));
1680+
1681+
for (i = 0; i < test_size; i++) {
1682+
if ((e = mp_prec_small_prime(to_test[i], &ret, &sieve)) != MP_OKAY) {
1683+
fprintf(stderr,"mp_prec_small_prime failed with \"%s\" at index %d\n",
1684+
mp_error_to_string(e), i);
1685+
goto LTM_ERR;
1686+
}
1687+
if (ret != tested[i]) {
1688+
fprintf(stderr,"mp_prec_small_prime failed for %lu. Said %lu but is %lu \n",
1689+
(unsigned long)to_test[i], (unsigned long)ret, (unsigned long)tested[i]);
1690+
goto LTM_ERR;
1691+
}
1692+
}
1693+
1694+
mp_sieve_clear(&sieve);
1695+
return EXIT_SUCCESS;
1696+
LTM_ERR:
1697+
mp_sieve_clear(&sieve);
1698+
return EXIT_FAILURE;
1699+
}
1700+
1701+
1702+
15321703
/*
15331704
Cannot test mp_exp(_d) without mp_root_n and vice versa.
15341705
So one of the two has to be tested from scratch.
@@ -2240,6 +2411,11 @@ static int unit_tests(int argc, char **argv)
22402411
T1(mp_prime_next_prime, MP_PRIME_NEXT_PRIME),
22412412
T1(mp_prime_rand, MP_PRIME_RAND),
22422413
T1(mp_rand, MP_RAND),
2414+
2415+
T1(mp_is_small_prime, MP_IS_SMALL_PRIME),
2416+
T1(mp_next_small_prime, MP_NEXT_SMALL_PRIME),
2417+
T1(mp_prec_small_prime, MP_PREC_SMALL_PRIME),
2418+
22432419
T1(mp_read_radix, MP_READ_RADIX),
22442420
T1(mp_read_write_ubin, MP_TO_UBIN),
22452421
T1(mp_read_write_sbin, MP_TO_SBIN),

‎doc/bn.tex

+170
Original file line numberDiff line numberDiff line change
@@ -2347,6 +2347,173 @@ \section{Random Primes}
23472347
\label{fig:primeopts}
23482348
\end{figure}
23492349

2350+
2351+
\section{Prime Sieve}
2352+
The prime sieve is implemented as a simple segmented Sieve of Eratosthenes.
2353+
2354+
This library needs some small sequential amounts for divisibility tests starting at two and
2355+
some random small primes for the Miller--Rabin test. That means we need a small base sieve
2356+
for quick results with a cold start and small segments to get random small primes fast.
2357+
2358+
The base sieve holds odd values only and even with a size of \texttt{4096} bytes it is
2359+
small enough to get build quickly, in about 50 $\mu$sec on the author's old machine.
2360+
2361+
The choice of the size for the individual segments varies more in the results. See table
2362+
\ref{table:sievebenchmarks} for some data. Size must be of the form $-1 + 2^n$.
2363+
Default size is \texttt{0xFFF}. It might be of interest that the biggest prime gap
2364+
below $2^{32}$ is $335$.
2365+
2366+
\begin{table}[h]
2367+
\begin{center}
2368+
\begin{tabular}{r c l}
2369+
\textbf{Segment Size (bits)} & \textbf{Random Access in $\mu$sec} & \textbf{Full Primesum} \\
2370+
\texttt{ 0x1F} & (60) & - \\
2371+
\texttt{ 0x3F} & (75) & - \\
2372+
\texttt{ 0x7F} & 63 & - \\
2373+
\texttt{ 0xFF} & 70 & 26 min \\
2374+
\texttt{ 0x1FF} & 73 & 13 min \\
2375+
\texttt{ 0x3FF} & 75 & 7 min 17 sec \\
2376+
\texttt{ 0x7FF} & 85 & 4 min 10 sec \\
2377+
\texttt{ 0xFFF} & 100 & 2 min 52 sec \\
2378+
\texttt{ 0x1FFF} & 120 & 1 min 45 sec \\
2379+
\texttt{ 0x3FFF} & 145 & 1 min 19 sec \\
2380+
\texttt{ 0x7FFF} & 200 & 1 min 04 sec \\
2381+
\texttt{ 0xFFFF} & 350 & 0 min 57 sec \\
2382+
\texttt{ 0xFFFFFFFF} & 300 & 0 min 35 sec
2383+
\end{tabular}
2384+
\caption{Average access times (warm) with the default base sieve and \texttt{MP-64BIT}} \label{table:sievebenchmarks}
2385+
\end{center}
2386+
\end{table}
2387+
2388+
The default sizes are in \texttt{tommath\_private.h}: \texttt{MP\_SIEVE\_PRIME\_MAX\_SQRT} is
2389+
the size of the base sieve and \texttt{MP\_SIEVE\_RANGE\_A\_B} is the size of the segment.
2390+
2391+
\subsection{Initialization and Clearing}
2392+
Initializing. It cannot fail because it only sets some default values. Memory is allocated later according to needs.
2393+
\index{mp\_sieve\_init}
2394+
\begin{alltt}
2395+
void mp_sieve_init(mp_sieve *sieve);
2396+
\end{alltt}
2397+
The function \texttt{mp\_sieve\_init} is equivalent to
2398+
\begin{alltt}
2399+
mp_sieve sieve = {NULL, NULL, 0};
2400+
\end{alltt}
2401+
2402+
Free the memory used by the sieve with
2403+
\index{mp\_sieve\_clear}
2404+
\begin{alltt}
2405+
void mp_sieve_clear(mp_sieve *sieve);
2406+
\end{alltt}
2407+
2408+
\subsection{Primality Test of Small Numbers}
2409+
Individual small numbers can be tested for primality with
2410+
\index{mp\_is\_small\_prime}
2411+
\begin{alltt}
2412+
int mp_is_small_prime(mp_sieve_prime n, bool *result,
2413+
mp_sieve *sieve);
2414+
\end{alltt}
2415+
It sets the variable \texttt{result} to \texttt{false} if the given number in \texttt{n} is a composite and
2416+
\texttt{true} if it is a prime.
2417+
2418+
It will return \texttt{MP\_SIEVE\_MAX\_REACHED} to flag the content of \texttt{result} as the last valid one.
2419+
2420+
The implementation of this function does all the heavy lifting--the building of the base sieve and the segment
2421+
if one is necessary. The base sieve stays, so this function can be used to ``warm up'' the sieve and, if
2422+
\texttt{n} is slightly larger than the upper limit of the base sieve, ``warm up'' the first segment, too.
2423+
2424+
\subsection{Find Adjacent Primes}
2425+
To find the prime bigger than a number \texttt{n} use
2426+
\index{mp\_next\_small\_prime}
2427+
\begin{alltt}
2428+
int mp_next_small_prime(mp_sieve_prime n, mp_sieve_prime *result,
2429+
mp_sieve *sieve);
2430+
\end{alltt}
2431+
and to find the one smaller than \texttt{n}
2432+
\index{mp\_prec\_small\_prime}
2433+
\begin{alltt}
2434+
int mp_prec_small_prime(mp_sieve_prime n, mp_sieve_prime *result,
2435+
mp_sieve *sieve);
2436+
\end{alltt}
2437+
Both functions return \texttt{n} if \texttt{n} is prime. Use \texttt{n + 1} to get
2438+
the prime after \texttt{n} if \texttt{n} is prime and \texttt{n - 1} to get the
2439+
the prime preceding \texttt{n} if \texttt{n} is prime.
2440+
2441+
\subsection{Useful Constants}
2442+
\begin{description}
2443+
\item[\texttt{MP\_SIEVE\_BIGGEST\_PRIME}] \texttt{read-only} The biggest prime the sieve can offer.
2444+
It is $4\,294\,967\,291$ for all \texttt{MP\_xBIT}.
2445+
2446+
\item[\texttt{mp\_sieve\_prime}] \texttt{read-only} The basic type for the numbers in the sieve. It
2447+
is \texttt{uint32\_t} for \texttt{MP\_16BIT}, \texttt{MP\_32BIT}; and
2448+
\texttt{uint64\_t} for \texttt{MP\_64BIT}..
2449+
2450+
\item[\texttt{MP\_SIEVE\_PR\_UINT}] Choses the correct macro from \texttt{inttypes.h} to print a\\
2451+
\texttt{mp\_sieve\_prime}. The header \texttt{inttypes.h} must be included before\\
2452+
\texttt{tommath.h} for it to work.
2453+
\end{description}
2454+
\subsection{Examples}\label{sec:spnexamples}
2455+
\subsubsection{Initialization and Clearing}
2456+
Using a sieve follows the same procedure as using a bigint:
2457+
\begin{alltt}
2458+
/* Declaration */
2459+
mp_sieve sieve;
2460+
2461+
/*
2462+
Initialization.
2463+
Cannot fail, only sets a handful of default values.
2464+
Memory allocation is done in the library itself
2465+
according to needs.
2466+
*/
2467+
mp_sieve_init(&sieve);
2468+
2469+
/* use the sieve */
2470+
2471+
/* Clean up */
2472+
mp_sieve_clear(&sieve);
2473+
\end{alltt}
2474+
\subsubsection{Primality Test}
2475+
A small program to test the input of a small number for primality.
2476+
\begin{alltt}
2477+
#include <stdlib.h>
2478+
#include <stdio.h>
2479+
#include <errno.h>
2480+
/*inttypes.h is needed for printing and must be included before tommath.h*/
2481+
#include <inttypes.h>
2482+
#include "tommath.h"
2483+
int main(int argc, char **argv)
2484+
{
2485+
mp_sieve_prime number;
2486+
mp_sieve sieve;
2487+
int e;
2488+
/* variable holding the result of mp_is_small_prime */
2489+
mp_sieve_prime result;
2490+
2491+
if (argc != 2) {
2492+
fprintf(stderr,"Usage %s number\textbackslash{}n", argv[0]);
2493+
exit(EXIT_FAILURE);
2494+
}
2495+
number = (mp_sieve_prime) strtoul(argv[1],NULL, 10);
2496+
if (errno == ERANGE) {
2497+
fprintf(stderr,"strtoul(l) failed: input out of range\textbackslash{}n");
2498+
goto LTM_ERR;
2499+
}
2500+
mp_sieve_init(&sieve);
2501+
if ((e = mp_is_small_prime(number, &result, &sieve)) != MP_OKAY) {
2502+
fprintf(stderr,"mp_is_small_prime failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
2503+
mp_error_to_string(e));
2504+
goto LTM_ERR;
2505+
}
2506+
printf("The number %" LTM_SIEVE_PR_UINT " is %s prime\textbackslash{}n",
2507+
number,(result)?"":"not");
2508+
2509+
mp_sieve_clear(&sieve);
2510+
exit(EXIT_SUCCESS);
2511+
LTM_ERR:
2512+
mp_sieve_clear(&sieve);
2513+
exit(EXIT_FAILURE);
2514+
}
2515+
\end{alltt}
2516+
23502517
\chapter{Random Number Generation}
23512518
\section{PRNG}
23522519
\index{mp\_rand}
@@ -2417,6 +2584,9 @@ \subsection{To ASCII}
24172584
mp_err mp_fwrite(const mp_int *a, int radix, FILE *stream);
24182585
\end{alltt}
24192586

2587+
2588+
2589+
24202590
\subsection{From ASCII}
24212591
\index{mp\_read\_radix}
24222592
\begin{alltt}

‎helper.pl

+1-1
Original file line numberDiff line numberDiff line change
@@ -58,7 +58,7 @@ sub check_source {
5858
push @{$troubles->{unwanted_free}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bfree\s*\(/;
5959
# and we probably want to also avoid the following
6060
push @{$troubles->{unwanted_memcpy}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemcpy\s*\(/ && $file !~ /s_mp_copy_digs.c/;
61-
push @{$troubles->{unwanted_memset}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemset\s*\(/ && $file !~ /s_mp_zero_buf.c/ && $file !~ /s_mp_zero_digs.c/;
61+
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/;
6262
push @{$troubles->{unwanted_memmove}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemmove\s*\(/;
6363
push @{$troubles->{unwanted_memcmp}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bmemcmp\s*\(/;
6464
push @{$troubles->{unwanted_strcmp}}, $lineno if $file =~ /^[^\/]+\.c$/ && $l =~ /\bstrcmp\s*\(/;

0 commit comments

Comments
 (0)