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