primesieve.c

   1 /* primesieve (BIT_ARRAY, N) -- Fills the BIT_ARRAY with a mask for primes up to N.
   2
   3 Contributed to the GNU project by Marco Bodrato.
   4
   5 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.
   6 IT IS ONLY SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.
   7 IN FACT, IT IS ALMOST GUARANTEED THAT IT WILL CHANGE OR
   8 DISAPPEAR IN A FUTURE GNU MP RELEASE.
   9
  10 Copyright 2010, 2011, 2012 Free Software Foundation, Inc.
  11
  12 This file is part of the GNU MP Library.
  13
  14 The GNU MP Library is free software; you can redistribute it and/or modify
  15 it under the terms of the GNU Lesser General Public License as published by
  16 the Free Software Foundation; either version 3 of the License, or (at your
  17 option) any later version.
  18
  19 The GNU MP Library is distributed in the hope that it will be useful, but
  20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  21 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  22 License for more details.
  23
  24 You should have received a copy of the GNU Lesser General Public License
  25 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  26
  27 #include "gmp.h"
  28 #include "gmp-impl.h"
  29
  30 /**************************************************************/
  31 /* Section macros: common macros, for mswing/fac/bin (&sieve) */
  32 /**************************************************************/
  33
  34 #define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)                 \
  35     __max_i = (end);                                            \
  36                                                                 \
  37     do {                                                        \
  38       ++__i;                                                    \
  39       if (((sieve)[__index] & __mask) == 0)                     \
  40         {                                                       \
  41           (prime) = id_to_n(__i)
  42
  43 #define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve)          \
  44   do {                                                          \
  45     mp_limb_t __mask, __index, __max_i, __i;                    \
  46                                                                 \
  47     __i = (start)-(off);                                        \
  48     __index = __i / GMP_LIMB_BITS;                              \
  49     __mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS);             \
  50     __i += (off);                                               \
  51                                                                 \
  52     LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)
  53
  54 #define LOOP_ON_SIEVE_STOP                                      \
  55         }                                                       \
  56       __mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1);       \
  57       __index += __mask & 1;                                    \
  58     }  while (__i <= __max_i)                                   \
  59
  60 #define LOOP_ON_SIEVE_END                                       \
  61     LOOP_ON_SIEVE_STOP;                                         \
  62   } while (0)
  63
  64 /*********************************************************/
  65 /* Section sieve: sieving functions and tools for primes */
  66 /*********************************************************/
  67
  68 #if 0
  69 static mp_limb_t
  70 bit_to_n (mp_limb_t bit) { return (bit*3+4)|1; }
  71 #endif
  72
  73 /* id_to_n (x) = bit_to_n (x-1) = (id*3+1)|1*/
  74 static mp_limb_t
  75 id_to_n  (mp_limb_t id)  { return id*3+1+(id&1); }
  76
  77 /* n_to_bit (n) = ((n-1)&(-CNST_LIMB(2)))/3U-1 */
  78 static mp_limb_t
  79 n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
  80
  81 #if 0
  82 static mp_size_t
  83 primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; }
  84 #endif
  85
  86 #if GMP_LIMB_BITS > 61
  87 #define SIEVE_SEED CNST_LIMB(0x3294C9E069128480)
  88 #define SEED_LIMIT 202
  89 #else
  90 #if GMP_LIMB_BITS > 30
  91 #define SIEVE_SEED CNST_LIMB(0x69128480)
  92 #define SEED_LIMIT 114
  93 #else
  94 #if GMP_LIMB_BITS > 15
  95 #define SIEVE_SEED CNST_LIMB(0x8480)
  96 #define SEED_LIMIT 54
  97 #else
  98 #if GMP_LIMB_BITS > 7
  99 #define SIEVE_SEED CNST_LIMB(0x80)
 100 #define SEED_LIMIT 34
 101 #else
 102 #define SIEVE_SEED CNST_LIMB(0x0)
 103 #define SEED_LIMIT 24
 104 #endif /* 7 */
 105 #endif /* 15 */
 106 #endif /* 30 */
 107 #endif /* 61 */
 108
 109 static void
 110 first_block_primesieve (mp_ptr bit_array, mp_limb_t n)
 111 {
 112   mp_size_t bits, limbs;
 113
 114   ASSERT (n > 4);
 115
 116   bits  = n_to_bit(n);
 117   limbs = bits / GMP_LIMB_BITS + 1;
 118
 119   /* FIXME: We can skip 5 too, filling with a 5-part pattern. */
 120   MPN_ZERO (bit_array, limbs);
 121   bit_array[0] = SIEVE_SEED;
 122
 123   if ((bits + 1) % GMP_LIMB_BITS != 0)
 124     bit_array[limbs-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
 125
 126   if (n > SEED_LIMIT) {
 127     mp_limb_t mask, index, i;
 128
 129     ASSERT (n > 49);
 130
 131     mask = 1;
 132     index = 0;
 133     i = 1;
 134     do {
 135       if ((bit_array[index] & mask) == 0)
 136         {
 137           mp_size_t step, lindex;
 138           mp_limb_t lmask;
 139           unsigned  maskrot;
 140
 141           step = id_to_n(i);
 142 /*        lindex = n_to_bit(id_to_n(i)*id_to_n(i)); */
 143           lindex = i*(step+1)-1+(-(i&1)&(i+1));
 144 /*        lindex = i*(step+1+(i&1))-1+(i&1); */
 145           if (lindex > bits)
 146             break;
 147
 148           step <<= 1;
 149           maskrot = step % GMP_LIMB_BITS;
 150
 151           lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
 152           do {
 153             bit_array[lindex / GMP_LIMB_BITS] |= lmask;
 154             lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
 155             lindex += step;
 156           } while (lindex <= bits);
 157
 158 /*        lindex = n_to_bit(id_to_n(i)*bit_to_n(i)); */
 159           lindex = i*(i*3+6)+(i&1);
 160
 161           lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
 162           for ( ; lindex <= bits; lindex += step) {
 163             bit_array[lindex / GMP_LIMB_BITS] |= lmask;
 164             lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
 165           };
 166         }
 167       mask = mask << 1 | mask >> (GMP_LIMB_BITS-1);
 168       index += mask & 1;
 169       i++;
 170     } while (1);
 171   }
 172 }
 173
 174 static void
 175 block_resieve (mp_ptr bit_array, mp_size_t limbs, mp_limb_t offset,
 176                       mp_srcptr sieve, mp_limb_t sieve_bits)
 177 {
 178   mp_size_t bits, step;
 179
 180   ASSERT (limbs > 0);
 181
 182   bits = limbs * GMP_LIMB_BITS - 1;
 183
 184   /* FIXME: We can skip 5 too, filling with a 5-part pattern. */
 185   MPN_ZERO (bit_array, limbs);
 186
 187   LOOP_ON_SIEVE_BEGIN(step,0,sieve_bits,0,sieve);
 188   {
 189     mp_size_t lindex;
 190     mp_limb_t lmask;
 191     unsigned  maskrot;
 192
 193 /*  lindex = n_to_bit(id_to_n(i)*id_to_n(i)); */
 194     lindex = __i*(step+1)-1+(-(__i&1)&(__i+1));
 195 /*  lindex = __i*(step+1+(__i&1))-1+(__i&1); */
 196     if (lindex > bits + offset)
 197       break;
 198
 199     step <<= 1;
 200     maskrot = step % GMP_LIMB_BITS;
 201
 202     if (lindex < offset)
 203       lindex += step * ((offset - lindex - 1) / step + 1);
 204
 205     lindex -= offset;
 206
 207     lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
 208     for ( ; lindex <= bits; lindex += step) {
 209       bit_array[lindex / GMP_LIMB_BITS] |= lmask;
 210       lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
 211     };
 212
 213 /*  lindex = n_to_bit(id_to_n(i)*bit_to_n(i)); */
 214     lindex = __i*(__i*3+6)+(__i&1);
 215     if (lindex > bits + offset)
 216       continue;
 217
 218     if (lindex < offset)
 219       lindex += step * ((offset - lindex - 1) / step + 1);
 220
 221     lindex -= offset;
 222
 223     lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
 224     for ( ; lindex <= bits; lindex += step) {
 225       bit_array[lindex / GMP_LIMB_BITS] |= lmask;
 226       lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
 227     };
 228   }
 229   LOOP_ON_SIEVE_END;
 230 }
 231
 232 #define BLOCK_SIZE 2048
 233
 234 /* Fills bit_array with the characteristic function of composite
 235    numbers up to the parameter n. I.e. a bit set to "1" represent a
 236    composite, a "0" represent a prime.
 237
 238    The primesieve_size(n) limbs pointed to by bit_array are
 239    overwritten. The returned value counts prime integers in the
 240    interval [4, n]. Note that n > 4.
 241
 242    Even numbers and multiples of 3 are excluded "a priori", only
 243    numbers equivalent to +/- 1 mod 6 have their bit in the array.
 244
 245    Once sieved, if the bit b is ZERO it represent a prime, the
 246    represented prime is bit_to_n(b), if the LSbit is bit 0, or
 247    id_to_n(b), if you call "1" the first bit.
 248  */
 249
 250 mp_limb_t
 251 gmp_primesieve (mp_ptr bit_array, mp_limb_t n)
 252 {
 253   mp_size_t size;
 254   mp_limb_t bits;
 255
 256   ASSERT (n > 4);
 257
 258   bits = n_to_bit(n);
 259   size = bits / GMP_LIMB_BITS + 1;
 260
 261   if (size > BLOCK_SIZE * 2) {
 262     mp_size_t off;
 263     off = BLOCK_SIZE + (size % BLOCK_SIZE);
 264     first_block_primesieve (bit_array, id_to_n (off * GMP_LIMB_BITS));
 265     for ( ; off < size; off += BLOCK_SIZE)
 266       block_resieve (bit_array + off, BLOCK_SIZE, off * GMP_LIMB_BITS, bit_array, off * GMP_LIMB_BITS - 1);
 267   } else {
 268     first_block_primesieve (bit_array, n);
 269   }
 270
 271   if ((bits + 1) % GMP_LIMB_BITS != 0)
 272     bit_array[size-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
 273
 274
 275   return size * GMP_LIMB_BITS - mpn_popcount (bit_array, size);
 276 }
 277
 278 #undef BLOCK_SIZE
 279 #undef SEED_LIMIT
 280 #undef SIEVE_SEED
 281 #undef LOOP_ON_SIEVE_END
 282 #undef LOOP_ON_SIEVE_STOP
 283 #undef LOOP_ON_SIEVE_BEGIN
 284 #undef LOOP_ON_SIEVE_CONTINUE