Removed support or 64-bit large sieve. Was just too big.

Author: czurnieden

demo/test.c

@@ -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);
@@ -1623,31 +1623,13 @@ static int test_mp_next_small_prime(void)
       1019879761, 72282701, 2048787577, 2058368113
    };
 
-/* TODO: No, remove 64-bit version, way too big and a bit unwieldy, too! */
-#if ( (defined MP_64BIT) && (defined MP_SIEVE_USE_LARGE_SIEVE) )
-   /* primesum up to 2^64
-      $ time /home/czurnieden/GITHUB/primesum/primesum 2^64
-       3879578600671960388666457126750869198
-
-       real 37m6,448s
-       user 107m17,056s
-       sys  0m12,152s
-      I think we should use something smaller.
-   */
-   /* const char *primesum_64 = "3879578600671960388666457126750869198"; */
-
-   const char *primesum_64 = "208139436659661";
-#else
 #ifdef BENCHMARK_SIEVE
    const char *primesum_32 = "425649736193687430";
 #else
    const char *primesum_32 = "202259606268580";
 #endif
-#endif
-
    mp_sieve_init(&sieve);
 
-
    test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime));
 
    for (i = 0; i < test_size; i++) {
@@ -1666,7 +1648,7 @@ static int test_mp_next_small_prime(void)
       }
 #ifdef BENCHMARK_SIEVE
       stopgt = gettime();
-      printf("Single random access for prime: %*u: ", 10, ret);
+      printf("Single random access for prime: %*"MP_SIEVE_PR_UINT": ", 10, ret);
       print_timer("",startgt, stopgt);
 #endif
    }
@@ -1675,37 +1657,32 @@ static int test_mp_next_small_prime(void)
 #ifdef BENCHMARK_SIEVE
    startgt = gettime();
 #endif
-#if ( (defined MP_64BIT) && (defined MP_SIEVE_USE_LARGE_SIEVE) )
-   for (p = 4293918720lu; ret < (mp_sieve_prime)4294997297lu;) {
-      DO(mp_next_small_prime(p, &ret, &sieve));
-      p = ret + 1;
-      mp_set_u64(&t, ret);
-      DO(mp_add(&primesum, &t, &primesum));
-   }
-#else
+
 #ifdef BENCHMARK_SIEVE
    for (p = 0ul; ret < (mp_sieve_prime)MP_SIEVE_BIGGEST_PRIME;) {
 #else
    for (p = 4293918720lu; ret < (mp_sieve_prime)MP_SIEVE_BIGGEST_PRIME;) {
 #endif
       DO(mp_next_small_prime(p, &ret, &sieve));
-      p = ret + 1;
+      p = ret + 1u;
+#ifdef MP_64BIT
+      mp_set_u64(&t, ret);
+#else
       mp_set_u32(&t, ret);
+#endif
       DO(mp_add(&primesum, &t, &primesum));
    }
-#endif
+
 #ifdef BENCHMARK_SIEVE
    stopgt = gettime();
-   printf("Sum of all primes between 0 and %u = ", MP_SIEVE_BIGGEST_PRIME);
+   printf("Sum of all primes between 0 and %"MP_SIEVE_PR_UINT" = ", (mp_sieve_prime)MP_SIEVE_BIGGEST_PRIME);
    DO(mp_fwrite(&primesum, 10, stdout));
-   print_timer(", computed in ",startgt, stopgt);
+   print_timer(", computed in \n",startgt, stopgt);
 #endif
-
-#if ( (defined MP_64BIT) && (defined MP_SIEVE_USE_LARGE_SIEVE) )
-   DO(mp_read_radix(&t, primesum_64, 10));
-#else
+   printf("Primesum computed: ");
+   DO(mp_fwrite(&primesum, 10, stdout));
+   puts("");
    DO(mp_read_radix(&t, primesum_32, 10));
-#endif
    EXPECT(mp_cmp(&primesum, &t) == MP_EQ);
 
    mp_sieve_clear(&sieve);

doc/bn.tex

@@ -2349,9 +2349,44 @@ \section{Random Primes}
 
 
 \section{Prime Sieve}
-The prime sieve is implemented as a simple segmented Sieve of Eratosthenes. It is only moderately optimized for
-space and runtime but should be small enough and also fast enough for almost all use-cases; quite quick for
-sequential access but relatively slow for random access.
+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 (\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.
@@ -2402,20 +2437,12 @@ \subsection{Find Adjacent Primes}
 
 \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 \texttt{MP\_16BIT}, \texttt{MP\_32BIT} and \texttt{MP\_64BIT}; and
- $18\,446\,744\,073\,709\,551\,557$ for \texttt{MP\_64BIT} if the macro\\
- \texttt{LTM\_SIEVE\_USE\_LARGE\_SIEVE} is defined.
+\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{MP\_64BIT}; and
- \texttt{uint64\_t} for \texttt{MP\_64BIT} if the macro \texttt{LTM\_SIEVE\_USE\_LARGE\_SIEVE} is defined.
-
-\item[\texttt{LTM\_SIEVE\_USE\_LARGE\_SIEVE}] \texttt{read-only} A compile time flag to make a large sieve.
- No advantage has been seen in using 64-bit integers if available except the ability to get a sieve up
-to $2^64$. But in this case the base sieve gets 0.25 Gibibytes large and the segments 0.5 Gibibytes
-(although you can change \texttt{LTM\_SIEVE\_RANGE\_A\_B} in \texttt{bn\_mp\_sieve.c} to get smaller segments)
-and it needs a long time to fill.
+ 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\\
@@ -2453,148 +2480,36 @@ \subsubsection{Primality Test}
 int main(int argc, char **argv)
 {
    mp_sieve_prime number;
-   mp_sieve *base = NULL;
-   mp_sieve *segment = NULL;
-   mp_sieve_prime single_segment_a = 0;
+   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}
-\subsubsection{Find Adjacent Primes}
-To sum up all primes up to and including \texttt{MP\_SIEVE\_BIGGEST\_PRIME} you might do something like:
-\begin{alltt}
-#include <stdlib.h>
-#include <stdio.h>
-#include <errno.h>
-#include <tommath.h>
-int main(int argc, char **argv)
-{
-   mp_sieve_prime number;
-   mp_sieve sieve;
-   mp_sieve_prime k, ret;
-   mp_int total, t;
-   int e;
-
-   if (argc != 2) {
-      fprintf(stderr,"Usage %s integer\textbackslash{}n", argv[0]);
-      exit(EXIT_FAILURE);
-   }
-
-   if ((e = mp_init_multi(&total, &t, NULL)) != MP_OKAY) {
-      fprintf(stderr,"mp_init_multi(segment): \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
-              mp_error_to_string(e));
-      goto LTM_ERR_1;
-   }
-   errno = 0;
-#if ( (defined MP_64BIT) && (defined LTM_SIEVE_USE_LARGE_SIEVE) )
-   number = (mp_sieve_prime) strtoull(argv[1],NULL, 10);
-#else
-   number = (mp_sieve_prime) strtoul(argv[1],NULL, 10);
-#endif
-   if (errno == ERANGE) {
-      fprintf(stderr,"strtoul(l) failed: input out of range\textbackslash{}n");
-      return EXIT_FAILURE
-   }
-
-   mp_sieve_init(&sieve);
-
-   for (k = 0, ret = 0; ret < number; k = ret) {
-      if ((e = mp_next_small_prime(k + 1, &ret, &sieve)) != MP_OKAY) {
-         if (e == LTM_SIEVE_MAX_REACHED) {
-#ifdef MP_64BIT
-            if ((e = mp_add_d(&total, (mp_digit) k, &total)) != MP_OKAY) {
-               fprintf(stderr,"mp_add_d (1) failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
-                       mp_error_to_string(e));
-               goto LTM_ERR;
-            }
-#else
-            if ((e = mp_set_long(&t, k)) != MP_OKAY) {
-               fprintf(stderr,"mp_set_long (1) failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
-                       mp_error_to_string(e));
-               goto LTM_ERR;
-            }
-            if ((e = mp_add(&total, &t, &total)) != MP_OKAY) {
-               fprintf(stderr,"mp_add (1) failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
-                       mp_error_to_string(e));
-               goto LTM_ERR;
-            }
-#endif
-            break;
-         }
-         fprintf(stderr,"mp_next_small_prime failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
-                 mp_error_to_string(e));
-         goto LTM_ERR;
-      }
-      /* The check if the prime is below the given limit
-       * cannot be done in the for-loop conditions because if
-       * it could we wouldn't need the sieve in the first place.
-       */
-      if (ret <= number) {
-#ifdef MP_64BIT
-         if ((e = mp_add_d(&total, (mp_digit) k, &total)) != MP_OKAY) {
-            fprintf(stderr,"mp_add_d failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
-                    mp_error_to_string(e));
-            goto LTM_ERR;
-         }
-#else
-         if ((e = mp_set_long(&t, k)) != MP_OKAY) {
-            fprintf(stderr,"mp_set_long failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
-                    mp_error_to_string(e));
-            goto LTM_ERR;
-         }
-         if ((e = mp_add(&total, &t, &total)) != MP_OKAY) {
-            fprintf(stderr,"mp_add failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n",
-            mp_error_to_string(e));
-            goto LTM_ERR;
-         }
-#endif
-      }
-   }
-   printf("total: ");
-   mp_fwrite(&total,10,stdout);
-   putchar('\textbackslash{}n');
-
-   mp_clear_multi(&total, &t, NULL);
    mp_sieve_clear(&sieve);
    exit(EXIT_SUCCESS);
 LTM_ERR:
-   mp_clear_multi(&total, &t, NULL);
    mp_sieve_clear(&sieve);
    exit(EXIT_FAILURE);
 }
 \end{alltt}
-It took about a minute on the authors machine from 2015 to get the expected $425\,649\,736\,193\,687\,430$ for the sum of all primes up to $2^{32}$, about the same runtime as Pari/GP version 2.9.4 (with a GMP-5.1.3 kernel).
 
 \chapter{Random Number Generation}
 \section{PRNG}