Move bitwise_primesieve from mpz/oddfac to a generally available function.

5 files changed, 292 insertions, 199 deletions

diff --git a/ChangeLog b/ChangeLog
index 5f15fa9e4..dde8939a3 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2012-04-19 Marco Bodrato <bodrato@mail.dm.unipi.it>
+
+	* primesieve.c: New file, with functions from mpz/oddfac_1.c .
+	* mpz/oddfac_1.c (bitwise_primesieve): Re-moved.
+	* Makefile.am (libgmp_la_SOURCES): Add primesieve.c .
+	* gmp-impl.h (gmp_primesieve): Declare.
+
 2012-04-17  Torbjorn Granlund  <tege@gmplib.org>
 
 	* mpn/x86_64/coreisbr/aorsmul_1.asm: Fix some DOS64 issues.
diff --git a/Makefile.am b/Makefile.am
index c2e84e806..4dfb3ed56 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -251,7 +251,7 @@ CXX_OBJECTS =								\
 libgmp_la_SOURCES = gmp-impl.h longlong.h				\
   assert.c compat.c errno.c extract-dbl.c invalid.c memory.c		\
   mp_bpl.c mp_clz_tab.c mp_dv_tab.c mp_minv_tab.c mp_get_fns.c mp_set_fns.c \
-  version.c nextprime.c
+  version.c nextprime.c primesieve.c
 EXTRA_libgmp_la_SOURCES = tal-debug.c tal-notreent.c tal-reent.c
 libgmp_la_DEPENDENCIES = @TAL_OBJECT@		\
   $(MPF_OBJECTS) $(MPZ_OBJECTS) $(MPQ_OBJECTS)	\
diff --git a/gmp-impl.h b/gmp-impl.h
index 528827821..602f2413a 100644
--- a/gmp-impl.h
+++ b/gmp-impl.h
@@ -1922,6 +1922,9 @@ __GMP_DECLSPEC void gmp_init_primesieve (gmp_primesieve_t *);
 #define gmp_nextprime __gmp_nextprime
 __GMP_DECLSPEC unsigned long int gmp_nextprime (gmp_primesieve_t *);
 
+#define gmp_primesieve __gmp_primesieve
+__GMP_DECLSPEC mp_limb_t gmp_primesieve (mp_ptr, mp_limb_t);
+
 
 #ifndef MUL_TOOM22_THRESHOLD
 #define MUL_TOOM22_THRESHOLD             30
diff --git a/mpz/oddfac_1.c b/mpz/oddfac_1.c
index 2dcb12827..5858088fe 100644
--- a/mpz/oddfac_1.c
+++ b/mpz/oddfac_1.c
@@ -101,202 +101,6 @@ n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
 static mp_size_t
 primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; }
 
-#if GMP_LIMB_BITS > 61
-#define SIEVE_SEED CNST_LIMB(0x3294C9E069128480)
-#define SEED_LIMIT 202
-#else
-#if GMP_LIMB_BITS > 30
-#define SIEVE_SEED CNST_LIMB(0x69128480)
-#define SEED_LIMIT 114
-#else
-#if GMP_LIMB_BITS > 15
-#define SIEVE_SEED CNST_LIMB(0x8480)
-#define SEED_LIMIT 54
-#else
-#if GMP_LIMB_BITS > 7
-#define SIEVE_SEED CNST_LIMB(0x80)
-#define SEED_LIMIT 34
-#else
-#define SIEVE_SEED CNST_LIMB(0x0)
-#define SEED_LIMIT 24
-#endif /* 7 */
-#endif /* 15 */
-#endif /* 30 */
-#endif /* 61 */
-
-static void
-first_block_primesieve (mp_ptr bit_array, mp_limb_t n)
-{
-  mp_size_t bits, limbs;
-
-  ASSERT (n > 4);
-
-  bits  = n_to_bit(n);
-  limbs = bits / GMP_LIMB_BITS + 1;
-
-  /* FIXME: We can skip 5 too, filling with a 5-part pattern. */
-  MPN_ZERO (bit_array, limbs);
-  bit_array[0] = SIEVE_SEED;
-
-  if ((bits + 1) % GMP_LIMB_BITS != 0)
-    bit_array[limbs-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
-
-  if (n > SEED_LIMIT) {
-    mp_limb_t mask, index, i;
-
-    ASSERT (n > 49);
-
-    mask = 1;
-    index = 0;
-    i = 1;
-    do {
-      if ((bit_array[index] & mask) == 0)
-	{
-	  mp_size_t step, lindex;
-	  mp_limb_t lmask;
-	  unsigned  maskrot;
-
-	  step = id_to_n(i);
-/*	  lindex = n_to_bit(id_to_n(i)*id_to_n(i)); */
-	  lindex = i*(step+1)-1+(-(i&1)&(i+1));
-/*	  lindex = i*(step+1+(i&1))-1+(i&1); */
-	  if (lindex > bits)
-	    break;
-
-	  step <<= 1;
-	  maskrot = step % GMP_LIMB_BITS;
-
-	  lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
-	  do {
-	    bit_array[lindex / GMP_LIMB_BITS] |= lmask;
-	    lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
-	    lindex += step;
-	  } while (lindex <= bits);
-
-/*	  lindex = n_to_bit(id_to_n(i)*bit_to_n(i)); */
-	  lindex = i*(i*3+6)+(i&1);
-
-	  lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
-	  for ( ; lindex <= bits; lindex += step) {
-	    bit_array[lindex / GMP_LIMB_BITS] |= lmask;
-	    lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
-	  };
-	}
-      mask = mask << 1 | mask >> (GMP_LIMB_BITS-1);
-      index += mask & 1;
-      i++;
-    } while (1);
-  }
-}
-
-static void
-block_resieve (mp_ptr bit_array, mp_size_t limbs, mp_limb_t offset,
-		      mp_srcptr sieve, mp_limb_t sieve_bits)
-{
-  mp_size_t bits, step;
-
-  ASSERT (limbs > 0);
-
-  bits = limbs * GMP_LIMB_BITS - 1;
-
-  /* FIXME: We can skip 5 too, filling with a 5-part pattern. */
-  MPN_ZERO (bit_array, limbs);
-
-  LOOP_ON_SIEVE_BEGIN(step,0,sieve_bits,0,sieve);
-  {
-    mp_size_t lindex;
-    mp_limb_t lmask;
-    unsigned  maskrot;
-
-/*  lindex = n_to_bit(id_to_n(i)*id_to_n(i)); */
-    lindex = __i*(step+1)-1+(-(__i&1)&(__i+1));
-/*  lindex = __i*(step+1+(__i&1))-1+(__i&1); */
-    if (lindex > bits + offset)
-      break;
-
-    step <<= 1;
-    maskrot = step % GMP_LIMB_BITS;
-
-    if (lindex < offset)
-      lindex += step * ((offset - lindex - 1) / step + 1);
-
-    lindex -= offset;
-
-    lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
-    for ( ; lindex <= bits; lindex += step) {
-      bit_array[lindex / GMP_LIMB_BITS] |= lmask;
-      lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
-    };
-
-/*  lindex = n_to_bit(id_to_n(i)*bit_to_n(i)); */
-    lindex = __i*(__i*3+6)+(__i&1);
-    if (lindex > bits + offset)
-      continue;
-
-    if (lindex < offset)
-      lindex += step * ((offset - lindex - 1) / step + 1);
-
-    lindex -= offset;
-
-    lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
-    for ( ; lindex <= bits; lindex += step) {
-      bit_array[lindex / GMP_LIMB_BITS] |= lmask;
-      lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
-    };
-  }
-  LOOP_ON_SIEVE_END;
-}
-
-#define BLOCK_SIZE 2048
-
-/* Fills bit_array with the characteristic function of composite
-   numbers up to the parameter n. I.e. a bit set to "1" represent a
-   composite, a "0" represent a prime.
-
-   The primesieve_size(n) limbs pointed to by bit_array are
-   overwritten. The returned value counts prime integers in the
-   interval [4, n]. Note that n > 4.
-
-   Even numbers and multiples of 3 are excluded "a priori", only
-   numbers equivalent to +/- 1 mod 6 have their bit in the array.
-
-   Once sieved, if the bit b is ZERO it represent a prime, the
-   represented prime is bit_to_n(b), if the LSbit is bit 0, or
-   id_to_n(b), if you call "1" the first bit.
- */
-
-static mp_limb_t
-bitwise_primesieve (mp_ptr bit_array, mp_limb_t n)
-{
-  mp_size_t size;
-  mp_limb_t bits;
-
-  ASSERT (n > 4);
-
-  bits = n_to_bit(n);
-  size = bits / GMP_LIMB_BITS + 1;
-
-  if (size > BLOCK_SIZE * 2) {
-    mp_size_t off;
-    off = BLOCK_SIZE + (size % BLOCK_SIZE);
-    first_block_primesieve (bit_array, id_to_n (off * GMP_LIMB_BITS));
-    for ( ; off < size; off += BLOCK_SIZE)
-      block_resieve (bit_array + off, BLOCK_SIZE, off * GMP_LIMB_BITS, bit_array, off * GMP_LIMB_BITS - 1);
-  } else {
-    first_block_primesieve (bit_array, n);
-  }
-
-  if ((bits + 1) % GMP_LIMB_BITS != 0)
-    bit_array[size-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
-
-
-  return size * GMP_LIMB_BITS - mpn_popcount (bit_array, size);
-}
-
-#undef BLOCK_SIZE
-#undef SEED_LIMIT
-#undef SIEVE_SEED
-
 /*********************************************************/
 /* Section mswing: 2-multiswing factorial                 */
 /*********************************************************/
@@ -591,7 +395,7 @@ mpz_oddfac_1 (mpz_ptr x, mp_limb_t n, unsigned flag)
 	  /* Put the sieve on the second half, it will be overwritten by the last mswing. */
 	  sieve = PTR (mswing) + size / 2 + 1;
 
-	  size = (bitwise_primesieve (sieve, n - 1) + 1) / log_n_max (n) + 1;
+	  size = (gmp_primesieve (sieve, n - 1) + 1) / log_n_max (n) + 1;
 
 	  factors = TMP_ALLOC_LIMBS (size);
 	  do {
@@ -633,4 +437,3 @@ mpz_oddfac_1 (mpz_ptr x, mp_limb_t n, unsigned flag)
 #undef LOG2C
 #undef FACTORS_PER_LIMB
 #undef FACTOR_LIST_STORE
-
diff --git a/primesieve.c b/primesieve.c
new file mode 100644
index 000000000..ea894dd77
--- /dev/null
+++ b/primesieve.c
@@ -0,0 +1,280 @@
+/* primesieve (BIT_ARRAY, N) -- Fills the BIT_ARRAY with a mask for primes up to N.
+
+Contributed to the GNU project by Marco Bodrato.
+
+THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.
+IT IS ONLY SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.
+IN FACT, IT IS ALMOST GUARANTEED THAT IT WILL CHANGE OR
+DISAPPEAR IN A FUTURE GNU MP RELEASE.
+
+Copyright 2010, 2011, 2012 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 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.
+
+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 Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/**************************************************************/
+/* Section macros: common macros, for mswing/fac/bin (&sieve) */
+/**************************************************************/
+
+#define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)			\
+    __max_i = (end);						\
+								\
+    do {							\
+      ++__i;							\
+      if (((sieve)[__index] & __mask) == 0)			\
+	{							\
+	  (prime) = id_to_n(__i)
+
+#define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve)		\
+  do {								\
+    mp_limb_t __mask, __index, __max_i, __i;			\
+								\
+    __i = (start)-(off);					\
+    __index = __i / GMP_LIMB_BITS;				\
+    __mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS);		\
+    __i += (off);						\
+								\
+    LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)
+
+#define LOOP_ON_SIEVE_STOP					\
+	}							\
+      __mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1);	\
+      __index += __mask & 1;					\
+    }  while (__i <= __max_i)					\
+
+#define LOOP_ON_SIEVE_END					\
+    LOOP_ON_SIEVE_STOP;						\
+  } while (0)
+
+/*********************************************************/
+/* Section sieve: sieving functions and tools for primes */
+/*********************************************************/
+
+static mp_limb_t
+bit_to_n (mp_limb_t bit) { return (bit*3+4)|1; }
+
+/* id_to_n (x) = bit_to_n (x-1) = (id*3+1)|1*/
+static mp_limb_t
+id_to_n  (mp_limb_t id)  { return id*3+1+(id&1); }
+
+/* n_to_bit (n) = ((n-1)&(-CNST_LIMB(2)))/3U-1 */
+static mp_limb_t
+n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
+
+static mp_size_t
+primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; }
+
+#if GMP_LIMB_BITS > 61
+#define SIEVE_SEED CNST_LIMB(0x3294C9E069128480)
+#define SEED_LIMIT 202
+#else
+#if GMP_LIMB_BITS > 30
+#define SIEVE_SEED CNST_LIMB(0x69128480)
+#define SEED_LIMIT 114
+#else
+#if GMP_LIMB_BITS > 15
+#define SIEVE_SEED CNST_LIMB(0x8480)
+#define SEED_LIMIT 54
+#else
+#if GMP_LIMB_BITS > 7
+#define SIEVE_SEED CNST_LIMB(0x80)
+#define SEED_LIMIT 34
+#else
+#define SIEVE_SEED CNST_LIMB(0x0)
+#define SEED_LIMIT 24
+#endif /* 7 */
+#endif /* 15 */
+#endif /* 30 */
+#endif /* 61 */
+
+static void
+first_block_primesieve (mp_ptr bit_array, mp_limb_t n)
+{
+  mp_size_t bits, limbs;
+
+  ASSERT (n > 4);
+
+  bits  = n_to_bit(n);
+  limbs = bits / GMP_LIMB_BITS + 1;
+
+  /* FIXME: We can skip 5 too, filling with a 5-part pattern. */
+  MPN_ZERO (bit_array, limbs);
+  bit_array[0] = SIEVE_SEED;
+
+  if ((bits + 1) % GMP_LIMB_BITS != 0)
+    bit_array[limbs-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
+
+  if (n > SEED_LIMIT) {
+    mp_limb_t mask, index, i;
+
+    ASSERT (n > 49);
+
+    mask = 1;
+    index = 0;
+    i = 1;
+    do {
+      if ((bit_array[index] & mask) == 0)
+	{
+	  mp_size_t step, lindex;
+	  mp_limb_t lmask;
+	  unsigned  maskrot;
+
+	  step = id_to_n(i);
+/*	  lindex = n_to_bit(id_to_n(i)*id_to_n(i)); */
+	  lindex = i*(step+1)-1+(-(i&1)&(i+1));
+/*	  lindex = i*(step+1+(i&1))-1+(i&1); */
+	  if (lindex > bits)
+	    break;
+
+	  step <<= 1;
+	  maskrot = step % GMP_LIMB_BITS;
+
+	  lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
+	  do {
+	    bit_array[lindex / GMP_LIMB_BITS] |= lmask;
+	    lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
+	    lindex += step;
+	  } while (lindex <= bits);
+
+/*	  lindex = n_to_bit(id_to_n(i)*bit_to_n(i)); */
+	  lindex = i*(i*3+6)+(i&1);
+
+	  lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
+	  for ( ; lindex <= bits; lindex += step) {
+	    bit_array[lindex / GMP_LIMB_BITS] |= lmask;
+	    lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
+	  };
+	}
+      mask = mask << 1 | mask >> (GMP_LIMB_BITS-1);
+      index += mask & 1;
+      i++;
+    } while (1);
+  }
+}
+
+static void
+block_resieve (mp_ptr bit_array, mp_size_t limbs, mp_limb_t offset,
+		      mp_srcptr sieve, mp_limb_t sieve_bits)
+{
+  mp_size_t bits, step;
+
+  ASSERT (limbs > 0);
+
+  bits = limbs * GMP_LIMB_BITS - 1;
+
+  /* FIXME: We can skip 5 too, filling with a 5-part pattern. */
+  MPN_ZERO (bit_array, limbs);
+
+  LOOP_ON_SIEVE_BEGIN(step,0,sieve_bits,0,sieve);
+  {
+    mp_size_t lindex;
+    mp_limb_t lmask;
+    unsigned  maskrot;
+
+/*  lindex = n_to_bit(id_to_n(i)*id_to_n(i)); */
+    lindex = __i*(step+1)-1+(-(__i&1)&(__i+1));
+/*  lindex = __i*(step+1+(__i&1))-1+(__i&1); */
+    if (lindex > bits + offset)
+      break;
+
+    step <<= 1;
+    maskrot = step % GMP_LIMB_BITS;
+
+    if (lindex < offset)
+      lindex += step * ((offset - lindex - 1) / step + 1);
+
+    lindex -= offset;
+
+    lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
+    for ( ; lindex <= bits; lindex += step) {
+      bit_array[lindex / GMP_LIMB_BITS] |= lmask;
+      lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
+    };
+
+/*  lindex = n_to_bit(id_to_n(i)*bit_to_n(i)); */
+    lindex = __i*(__i*3+6)+(__i&1);
+    if (lindex > bits + offset)
+      continue;
+
+    if (lindex < offset)
+      lindex += step * ((offset - lindex - 1) / step + 1);
+
+    lindex -= offset;
+
+    lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
+    for ( ; lindex <= bits; lindex += step) {
+      bit_array[lindex / GMP_LIMB_BITS] |= lmask;
+      lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
+    };
+  }
+  LOOP_ON_SIEVE_END;
+}
+
+#define BLOCK_SIZE 2048
+
+/* Fills bit_array with the characteristic function of composite
+   numbers up to the parameter n. I.e. a bit set to "1" represent a
+   composite, a "0" represent a prime.
+
+   The primesieve_size(n) limbs pointed to by bit_array are
+   overwritten. The returned value counts prime integers in the
+   interval [4, n]. Note that n > 4.
+
+   Even numbers and multiples of 3 are excluded "a priori", only
+   numbers equivalent to +/- 1 mod 6 have their bit in the array.
+
+   Once sieved, if the bit b is ZERO it represent a prime, the
+   represented prime is bit_to_n(b), if the LSbit is bit 0, or
+   id_to_n(b), if you call "1" the first bit.
+ */
+
+mp_limb_t
+gmp_primesieve (mp_ptr bit_array, mp_limb_t n)
+{
+  mp_size_t size;
+  mp_limb_t bits;
+
+  ASSERT (n > 4);
+
+  bits = n_to_bit(n);
+  size = bits / GMP_LIMB_BITS + 1;
+
+  if (size > BLOCK_SIZE * 2) {
+    mp_size_t off;
+    off = BLOCK_SIZE + (size % BLOCK_SIZE);
+    first_block_primesieve (bit_array, id_to_n (off * GMP_LIMB_BITS));
+    for ( ; off < size; off += BLOCK_SIZE)
+      block_resieve (bit_array + off, BLOCK_SIZE, off * GMP_LIMB_BITS, bit_array, off * GMP_LIMB_BITS - 1);
+  } else {
+    first_block_primesieve (bit_array, n);
+  }
+
+  if ((bits + 1) % GMP_LIMB_BITS != 0)
+    bit_array[size-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
+
+
+  return size * GMP_LIMB_BITS - mpn_popcount (bit_array, size);
+}
+
+#undef BLOCK_SIZE
+#undef SEED_LIMIT
+#undef SIEVE_SEED
+#undef LOOP_ON_SIEVE_END
+#undef LOOP_ON_SIEVE_STOP
+#undef LOOP_ON_SIEVE_BEGIN
+#undef LOOP_ON_SIEVE_CONTINUE