gen-sieve.c: New file to generate a small pre-sieved array

1 files changed, 123 insertions, 0 deletions

diff --git a/gen-sieve.c b/gen-sieve.c
new file mode 100644
index 000000000..8fac3413f
--- /dev/null
+++ b/gen-sieve.c
@@ -0,0 +1,123 @@
+/* Generate primesieve data.
+
+   Contributed to the GNU project by Marco Bodrato.
+
+Copyright 2021 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of either:
+
+  * the GNU Lesser General Public License as published by the Free
+    Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+or
+
+  * the GNU General Public License as published by the Free Software
+    Foundation; either version 2 of the License, or (at your option) any
+    later version.
+
+or both in parallel, as here.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received copies of the GNU General Public License and the
+GNU Lesser General Public License along with the GNU MP Library.  If not,
+see https://www.gnu.org/licenses/.  */
+
+#include <stdio.h>
+#include "bootstrap.c"
+
+static int
+bit_to_n (int bit) { return (bit*3+4)|1; }
+
+int
+generate (int limb_bits, int limit)
+{
+  mpz_t  limb;
+  int i, lb, pc, c, totpc, maxprime;
+
+  mpz_init (limb);
+
+  printf ("/* This file generated by gen-sieve.c - DO NOT EDIT. */\n");
+  printf ("\n");
+  printf ("#if GMP_LIMB_BITS != %d\n", limb_bits);
+  printf ("Error, error, this data is for %d bits\n", limb_bits);
+  printf ("#endif\n");
+  printf ("\n");
+  printf ("#define PRIMESIEVE_INIT_TABLE ");
+
+  maxprime = 3;
+  lb = pc = c = totpc = 0;
+  for (i = 0; i < limit; i++)
+    {
+      if (! isprime (bit_to_n (i)))
+	mpz_setbit (limb, lb);
+      else
+	maxprime = bit_to_n (i), ++pc;
+      ++lb;
+      if (lb == limb_bits)
+	{
+	  ++c;
+	  printf ("\\\n\tCNST_LIMB (0x");
+	  mpz_out_str (stdout, -16, limb);
+	  printf ("),\t/* %d - %d (%d primes) */\t", bit_to_n (i + 1 - limb_bits),
+		  bit_to_n (i + 1) - 1, pc);
+	  totpc += pc;
+	  lb = pc = 0;
+	  mpz_set_ui (limb, 0);
+	}
+    }
+
+  if ((mpz_sgn (limb) | lb | pc) != 0)
+    {
+      printf ("\ngen-sieve: Internal error, during generate (%d, %d).\n", limb_bits,  limit);
+      abort();
+    }
+
+  mpz_clear (limb);
+
+  printf ("\n");
+  printf ("#define PRIMESIEVE_NUMBEROF_TABLE %d\n", c);
+
+  printf ("/* #define PRIMESIEVE_PRIMES_IN_TABLE %d */\n", totpc);
+  printf ("/* #define PRIMESIEVE_HIGHEST_PRIME %d */\n", maxprime);
+  printf ("/* #define PRIMESIEVE_FIRST_UNCHECKED %d */\n", bit_to_n (limit));
+
+  return c;
+}
+
+/* 5*2 = 10
+   7*2 = 14
+   5*7*2 = 70 (2*35, 3*24, 4*18, 5*14...)
+   5*11*2 = 110 (2*55, 3*37, 4*28, 5*22...)
+   5*13*2 = 130 (2*65, 3*44, 4*33, 5*26...)
+   7*11*2 = 154 (2*77, 3*52, 4*39, 5*31...)
+   7*13*2 = 182 (2*91, 3*61, 4*46, 5*37...)
+*/
+
+int
+main (int argc, char *argv[])
+{
+  int  limb_bits, limit;
+
+  if (argc != 2)
+    {
+      fprintf (stderr, "Usage: gen-sieve <limbbits>\n");
+      exit (1);
+    }
+
+  limb_bits = atoi (argv[1]);
+
+  limit = 64 * 28; /* bits in the presieved sieve */
+  if (limit % limb_bits != 0)
+    limit += limb_bits - limit % limb_bits;
+  generate (limb_bits, limit);
+
+  return 0;
+}