From 7853f587e90f5e9897872b8adb496855eb2fc0c4 Mon Sep 17 00:00:00 2001
From: "Michael I. Bushnell" <mib@gnu.org>
Date: Tue, 22 Mar 1994 22:58:41 +0000
Subject: Formerly primes.c.~2~

---
 proc/primes.c | 77 ++++++++++++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 71 insertions(+), 6 deletions(-)

(limited to 'proc')

diff --git a/proc/primes.c b/proc/primes.c
index cf6d4ba7..896f84fb 100644
--- a/proc/primes.c
+++ b/proc/primes.c
@@ -15,6 +15,9 @@
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
 
+#include <stdlib.h>
+#include <string.h>
+
 /* Array of prime numbers */
 int *primes;
 
@@ -26,16 +29,78 @@ int nprimes;
 
 /* Initialize primes  */
 void
-initprimes (void)
+initprimes ()
 {
-  primessize = 1;
-  nprimes = 1;
-  primes = malloc (sizeof (int) * 1);
-  *primes = 2;
+  primessize = 2;
+  nprimes = 2;
+  primes = malloc (sizeof (int) * 2);
+  primes[0] = 2;
+  primes[1] = 3;
 }
 
+/* Make the array of primes larger than it is right now.  */
+void
+growprimes ()
+{
+  int *iscomp;
+  int nints;
+  int lastprime = primes[nprimes - 1];
+  int i, j;
+  
+  nints = lastprime * lastprime;
+  iscomp = alloca (sizeof (int) * nints);
+  bzero (iscomp, sizeof (int) * nints);
+
+  for (i = 0; i < nprimes; i++)
+    for (j = primes[i] * 2; j < nints; j += primes[i])
+      iscomp[j] = 1;
+
+  for (i = lastprime; i < nints; i++)
+    {
+      if (nprimes == primessize)
+	{
+	  primes = realloc (primes, primessize * sizeof (int) * 2);
+	  primessize *= 2;
+	}
+      if (!iscomp[i])
+	primes[nprimes++] = i;
+    }
+}  
+
 /* Return the next prime greater than or equal to n. */
 int 
 nextprime (int n)
 {
-  if (n >= primes[nprimes])
+  int p;
+  int low, high;
+  
+  if (n < primes[0])
+    return primes[0];
+
+  while (n > primes[nprimes - 1])
+    growprimes ();
+
+  /* Binary search */
+  low = 0;
+  high = nprimes - 1;
+  p = high / 2;
+
+  /* This works because nprimes is always at least 2. */
+  while (primes[p - 1] >= n || primes[p] < n)
+    {
+      if (n > primes[p])
+	low = p;
+      else
+	high = p;
+      p = ((high - low) / 2) + low;
+    }
+  
+  return primes[p];
+}
+
+  
+
+
+  
+    
+  
-- 
cgit v1.2.3