Complete rewrite.

1 files changed, 311 insertions, 77 deletions

diff --git a/demos/primes.c b/demos/primes.c
index af9dce3da..1a50402ab 100644
--- a/demos/primes.c
+++ b/demos/primes.c
@@ -1,6 +1,8 @@
 /* List and count primes.
+   Written by tege while on holiday in Rodupp, August 2001.
+   Between 10 and 500 times faster than previous program.
 
-Copyright 2000 Free Software Foundation, Inc.
+Copyright 2001 Free Software Foundation, Inc.
 
 This program is free software; you can redistribute it and/or modify it under
 the terms of the GNU General Public License as published by the Free Software
@@ -16,29 +18,92 @@ this program; if not, write to the Free Software Foundation, Inc., 59 Temple
 Place - Suite 330, Boston, MA 02111-1307, USA.  */
 
 #include <stdlib.h>
-#include <string.h>
 #include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+/* IDEAS:
+ * Do not fill primes[] with real primes when the range [fr,to] is small,
+   when fr,to are relatively large.  Fill primes[] with odd numbers instead.
+   [Probably a bad idea, since the primes[] array would become very large.]
+ * Separate small primes and large primes when sieving.  Either the Montgomery
+   way (i.e., having a large array a multiple of L1 cache size), or just
+   separate loops for primes <= S and primes > S.  The latter primes do not
+   require an inner loop, since they will touch the sieving array at most once.
+ * Pre-fill sieving array with an appropriately aligned ...00100100... pattern,
+   then omit 3 from primes array.  (May require similar special handling of 3
+   as we now have for 2.)
+ * A large SIEVE_LIMIT currently implies very large memory usage, mainly due
+   to the sieving array in make_primelist, but also because of the primes[]
+   array.  We might want to stage the program, using sieve_region/find_primes
+   to build primes[].  Make report() a function pointer, as part of achieving
+   this.
+ * Store primes[] as two arrays, one array with primes represented as delta
+   values using just 8 bits (if gaps are too big, store bogus primes!)
+   and one array with "rem" values.  The latter needs 32-bit values.
+ * A new entry point, mpz_probab_prime_likely_p, would be useful.
+ * Improve command line syntax and versaitlity.  "primes -f FROM -t TO",
+   allow either to be omitted for open interval.  (But disallow
+   "primes -c -f FROM" since that would be infinity.)  Allow printing a
+   limited *number* of primes using syntax like "primes -f FROM -n NUMBER".
+ */
+
 #include "gmp.h"
 
-/* Sieve table size */
-#define ST_SIZE 50000
-/* Largest prime to sieve with */
-#define MAX_S_PRIME 32000
+struct primes
+{
+  unsigned int prime;
+  signed int rem;
+};
+
+struct primes *primes;
+unsigned long n_primes;
+
+void find_primes (unsigned char *, mpz_t, unsigned long, mpz_t);
+void sieve_region (unsigned char *, mpz_t, unsigned long);
+void make_primelist (unsigned long);
+
+int flag_print = 1;
+int flag_count = 0;
+int flag_maxgap = 0;
+unsigned long maxgap = 0;
+unsigned long total_primes = 0;
 
-main (int argc, char **argv)
+void
+report (mpz_t prime)
+{
+  total_primes += 1;
+  if (flag_print)
+    {
+      mpz_out_str (stdout, 10, prime);
+      printf ("\n");
+    }
+  if (flag_maxgap)
+    {
+      static unsigned long prev_prime_low = 0;
+      unsigned long gap;
+      if (prev_prime_low != 0)
+	{
+	  gap = mpz_get_ui (prime) - prev_prime_low;
+	  if (maxgap < gap)
+	    maxgap = gap;
+	}
+      prev_prime_low = mpz_get_ui (prime);
+    }
+}
+
+int
+main (int argc, char *argv[])
 {
   char *progname = argv[0];
-  mpz_t r0, r1;			/* range */
-  mpz_t cur;
-  mpz_t max_s_prime_squared;
-  unsigned char *st;
-  unsigned long i, ii;
-  unsigned long nprimes = 0;
-  unsigned long last;
-  int flag_print = 1;
-  int flag_count = 0;
-
-  st = malloc (ST_SIZE);
+  mpz_t fr, to;
+  mpz_t fr2, to2;
+  unsigned long sieve_lim;
+  unsigned long est_n_primes;
+  unsigned char *s;
+  mpz_t tmp;
+  mpz_t siev_sqr_lim;
 
   while (argc != 1)
     {
@@ -54,99 +119,268 @@ main (int argc, char **argv)
 	  argv++;
 	  argc--;
 	}
+      else if (strcmp (argv[1], "-g") == 0)
+	{
+	  flag_maxgap = 1;
+	  argv++;
+	  argc--;
+	}
       else
 	break;
     }
 
-  if (flag_count)
+  if (flag_count || flag_maxgap)
     flag_print--;		/* clear unless an explicit -p  */
 
-  mpz_init (r0);
-  mpz_init (r1);
-  mpz_init (cur);
-  mpz_init_set_ui (max_s_prime_squared, MAX_S_PRIME);
-  mpz_mul_ui (max_s_prime_squared, max_s_prime_squared, MAX_S_PRIME);
+  mpz_init (fr);
+  mpz_init (to);
+  mpz_init (fr2);
+  mpz_init (to2);
 
-  if (argc == 2)
-    {
-      mpz_set_ui (r0, 2);
-      mpz_set_str (r1, argv[1], 0);
-    }
-  else if (argc == 3)
+  if (argc == 3)
     {
-      mpz_set_str (r0, argv[1], 0);
+      mpz_set_str (fr, argv[1], 0);
       if (argv[2][0] == '+')
 	{
-	  mpz_set_str (r1, argv[2] + 1, 0);
-	  mpz_add (r1, r1, r0);
+	  mpz_set_str (to, argv[2] + 1, 0);
+	  mpz_add (to, to, fr);
 	}
       else
-	mpz_set_str (r1, argv[2], 0);
+	mpz_set_str (to, argv[2], 0);
+    }
+  else if (argc == 2)
+    {
+      mpz_set_ui (fr, 0);
+      mpz_set_str (to, argv[1], 0);
     }
   else
     {
-      fprintf (stderr, "usage: %s [-c] [-g] [from [+]]to\n", progname);
+      fprintf (stderr, "usage: %s [-c] [-p] [-g] [from [+]]to\n", progname);
       exit (1);
     }
 
-  if (mpz_cmp_ui (r0, 2) < 0)
-    mpz_set_ui (r0, 2);
+  mpz_set (fr2, fr);
+  if (mpz_cmp_ui (fr2, 3) < 0)
+    {
+      mpz_set_ui (fr2, 2);
+      report (fr2);
+      mpz_set_ui (fr2, 3);
+    }
+  mpz_setbit (fr2, 0);				/* make odd */
+  mpz_sub_ui (to2, to, 1);
+  mpz_setbit (to2, 0);				/* make odd */
+
+  mpz_init (tmp);
+  mpz_init (siev_sqr_lim);
 
-  if ((mpz_get_ui (r0) & 1) == 0)
+  mpz_sqrt (tmp, to2);
+#define SIEVE_LIMIT 10000000
+  if (mpz_cmp_ui (tmp, SIEVE_LIMIT) < 0)
     {
-      if (mpz_cmp_ui (r0, 2) == 0)
+      sieve_lim = mpz_get_ui (tmp);
+    }
+  else
+    {
+      sieve_lim = SIEVE_LIMIT;
+      mpz_sub (tmp, to2, fr2);
+      if (mpz_cmp_ui (tmp, sieve_lim) < 0)
+	sieve_lim = mpz_get_ui (tmp);	/* limit sieving for small ranges */
+    }
+  mpz_set_ui (siev_sqr_lim, sieve_lim + 1);
+  mpz_mul_ui (siev_sqr_lim, siev_sqr_lim, sieve_lim + 1);
+
+  est_n_primes = (size_t) (sieve_lim / log((double) sieve_lim) * 1.13) + 10;
+  primes = malloc (est_n_primes * sizeof primes[0]);
+  make_primelist (sieve_lim);
+  assert (est_n_primes >= n_primes);
+
+#if DEBUG
+  printf ("sieve_lim = %lu\n", sieve_lim);
+  printf ("n_primes = %lu (3..%u)\n",
+	  n_primes, primes[n_primes - 1].prime);
+#endif
+
+#define S (1 << 15)		/* FIXME: Figure out L1 cache size */
+  s = malloc (S/2);
+  while (mpz_cmp (fr2, to2) <= 0)
+    {
+      unsigned long rsize;
+      rsize = S;
+      mpz_add_ui (tmp, fr2, rsize);
+      if (mpz_cmp (tmp, to2) > 0)
 	{
-	  if (flag_print)
-	    puts ("2");
-	  nprimes++;
+	  mpz_sub (tmp, to2, fr2);
+	  rsize = mpz_get_ui (tmp) + 2;
 	}
-      mpz_add_ui (r0, r0, 1);
+#if DEBUG
+      printf ("Sieving region ["); mpz_out_str (stdout, 10, fr2);
+      printf (","); mpz_add_ui (tmp, fr2, rsize - 2);
+      mpz_out_str (stdout, 10, tmp); printf ("]\n");
+#endif
+      sieve_region (s, fr2, rsize);
+      find_primes (s, fr2, rsize / 2, siev_sqr_lim);
+
+      mpz_add_ui (fr2, fr2, S);
     }
+  free (s);
 
-  mpz_set (cur, r0);
+  if (flag_count)
+    printf ("Pi(interval) = %lu\n", total_primes);
+
+  if (flag_maxgap)
+    printf ("max gap: %lu\n", maxgap);
+
+  return 0;
+}
+
+/* Find primes in region [fr,fr+rsize).  Requires that fr is odd and that
+   rsize is even.  The sieving array s should be aligned for "long int" and
+   have rsize/2 entries, rounded up to the nearest multiple of "long int".  */
+void
+sieve_region (unsigned char *s, mpz_t fr, unsigned long rsize)
+{
+  unsigned long ssize = rsize / 2;
+  unsigned long start, start2, prime;
+  unsigned long i;
+  mpz_t tmp;
+
+  mpz_init (tmp);
+
+#if 0
+  /* initialize sieving array */
+  for (ii = 0; ii < (ssize + sizeof (long) - 1) / sizeof (long); ii++)
+    ((long *) s) [ii] = ~0L;
+#else
+  {
+    signed long k;
+    long *se = (long *) (s + ((ssize + sizeof (long) - 1) & -sizeof (long)));
+    for (k = -((ssize + sizeof (long) - 1) / sizeof (long)); k < 0; k++)
+      se[k] = ~0L;
+  }
+#endif
 
-  while (mpz_cmp (cur, r1) <= 0)
+  for (i = 0; i < n_primes; i++)
     {
-      memset (st, 1, ST_SIZE);
-      for (i = 3; i < MAX_S_PRIME; i += 2)
+      prime = primes[i].prime;
+
+      if (primes[i].rem >= 0)
 	{
-	  unsigned long start;
-	  start = mpz_tdiv_ui (cur, i);
-	  if (start != 0)
-	    start = i - start;
-	  
-	  if (mpz_cmp_ui (cur, i - start) == 0)
-	    start += i;
-	  for (ii = start; ii < ST_SIZE; ii += i)
-	    st[ii] = 0;
+	  start2 = primes[i].rem;
 	}
-      last = 0;
-      for (ii = 0; ii < ST_SIZE; ii += 2)
+      else
+	{
+	  mpz_set_ui (tmp, prime);
+	  mpz_mul_ui (tmp, tmp, prime);
+	  if (mpz_cmp (fr, tmp) <= 0)
+	    {
+	      mpz_sub (tmp, tmp, fr);
+	      if (mpz_cmp_ui (tmp, 2 * ssize) > 0)
+		break;		/* avoid overflow at next line, also speedup */
+	      start = mpz_get_ui (tmp);
+	    }
+	  else
+	    {
+	      start = (prime - mpz_tdiv_ui (fr, prime)) % prime;
+	      if (start % 2 != 0)
+		start += prime;		/* adjust if even divisable */
+	    }
+	  start2 = start / 2;
+	}
+
+#if 0
+      for (ii = start2; ii < ssize; ii += prime)
+	s[ii] = 0;
+      primes[i].rem = ii - ssize;
+#else
+      {
+	signed long k;
+	unsigned char *se = s + ssize; /* point just beyond sieving range */
+	for (k = start2 - ssize; k < 0; k += prime)
+	  se[k] = 0;
+	primes[i].rem = k;
+      }
+#endif
+    }
+  mpz_clear (tmp);
+}
+
+/* Find primes in region [fr,fr+rsize), using the previously sieved s[].  */
+void
+find_primes (unsigned char *s, mpz_t  fr, unsigned long ssize,
+	     mpz_t siev_sqr_lim)
+{
+  unsigned long j, ij;
+  mpz_t tmp;
+
+  mpz_init (tmp);
+  for (j = 0; j < (ssize + sizeof (long) - 1) / sizeof (long); j++)
+    {
+      if (((long *) s) [j] != 0)
 	{
-	  if (st[ii] != 0)
+	  for (ij = 0; ij < sizeof (long); ij++)
 	    {
-	      mpz_add_ui (cur, cur, ii - last);
-	      last = ii;
-	      if (mpz_cmp (cur, r1) > 0)
-		goto done;
-	      /* Perform probabilistic prime test, but only if we're in a range
-		 were the sieving hasn't already proved this number prime.  */
-	      if (mpz_cmp (cur, max_s_prime_squared) < 0
-		  || mpz_probab_prime_p (cur, 3))
+	      if (s[j * sizeof (long) + ij] != 0)
 		{
-		  nprimes++;
-		  if (flag_print)
-		    {
-		      mpz_out_str (stdout, 10, cur); puts ("");
-		    }
+		  if (j * sizeof (long) + ij >= ssize)
+		    goto out;
+		  mpz_add_ui (tmp, fr, (j * sizeof (long) + ij) * 2);
+		  if (mpz_cmp (tmp, siev_sqr_lim) < 0 ||
+		      mpz_probab_prime_p (tmp, 3))
+		    report (tmp);
 		}
 	    }
 	}
-      mpz_add_ui (cur, cur, ST_SIZE - last);
     }
- done:
-  if (flag_count)
-    printf ("Pi(interval) = %lu\n", nprimes);
+ out:
+  mpz_clear (tmp);
+}
+
+/* Generate a lits of primes and store in the global array primes[].  */
+void
+make_primelist (unsigned long maxprime)
+{
+#if 1
+  unsigned char *s;
+  unsigned long ssize = maxprime / 2;
+  unsigned long i, ii, j;
 
-  exit (0);
+  s = malloc (ssize);
+  memset (s, ~0, ssize);
+  for (i = 3; ; i += 2)
+    {
+      unsigned long isqr = i * i;
+      if (isqr >= maxprime)
+	break;
+      if (s[i * i / 2 - 1] == 0)
+	continue;				/* only sieve with primes */
+      for (ii = i * i / 2 - 1; ii < ssize; ii += i)
+	s[ii] = 0;
+    }
+  n_primes = 0;
+  for (j = 0; j < ssize; j++)
+    {
+      if (s[j] != 0)
+	{
+	  primes[n_primes].prime = j * 2 + 3;
+	  primes[n_primes].rem = -1;
+	  n_primes++;
+	}
+    }
+  /* FIXME: This should not be needed if fencepost errors were fixed... */
+  if (primes[n_primes - 1].prime > maxprime)
+    n_primes--;
+  free (s);
+#else
+  unsigned long i;
+  n_primes = 0;
+  for (i = 3; i <= maxprime; i += 2)
+    {
+      if (i < 7 || (i % 3 != 0 && i % 5 != 0 && i % 7 != 0))
+	{
+	  primes[n_primes].prime = i;
+	  primes[n_primes].rem = -1;
+	  n_primes++;
+	}
+    }
+#endif
 }