1 /// 2 /// @file Wheel.hpp 3 /// @brief Wheel factorization is used to skip multiles of 4 /// small primes in the sieve of Eratosthenes. 5 /// 6 /// Copyright (C) 2020 Kim Walisch, <kim.walisch@gmail.com> 7 /// 8 /// This file is distributed under the BSD License. See the COPYING 9 /// file in the top level directory. 10 /// 11 12 #ifndef WHEEL_HPP 13 #define WHEEL_HPP 14 15 #include <stdint.h> 16 #include <algorithm> 17 #include <cassert> 18 19 namespace primesieve { 20 21 /// The WheelInit data structure is used to calculate the 22 /// first multiple >= start of each sieving prime 23 /// 24 struct WheelInit 25 { 26 uint8_t nextMultipleFactor; 27 uint8_t wheelIndex; 28 }; 29 30 extern const WheelInit wheel30Init[30]; 31 extern const WheelInit wheel210Init[210]; 32 33 /// The WheelElement data structure is used to skip multiples 34 /// of small primes using wheel factorization 35 /// 36 struct WheelElement 37 { 38 /// Bitmask used to unset the bit corresponding to the current 39 /// multiple of a SievingPrime object 40 uint8_t unsetBit; 41 /// Factor used to calculate the next multiple of a sieving prime 42 /// that is not divisible by any of the wheel factors 43 uint8_t nextMultipleFactor; 44 /// Overflow needed to correct the next multiple index 45 /// (due to sievingPrime = prime / 30) 46 uint8_t correct; 47 /// Used to calculate the next wheel index: 48 /// wheelIndex += next; 49 int8_t next; 50 }; 51 52 extern const WheelElement wheel30[8*8]; 53 extern const WheelElement wheel210[48*8]; 54 55 /// The abstract Wheel class is used skip multiples of small 56 /// primes in the sieve of Eratosthenes. The EratSmall, 57 /// EratMedium and EratBig classes are derived from Wheel. 58 /// 59 template <int MODULO, int SIZE, const WheelElement* WHEEL, const WheelInit* INIT> 60 class Wheel 61 { 62 public: 63 /// Add a new sieving prime to the sieving algorithm. 64 /// Calculate the first multiple > segmentLow of prime and 65 /// the position within the sieve array of that multiple 66 /// and its wheel index. When done store the sieving prime. 67 /// addSievingPrime(uint64_t prime,uint64_t segmentLow)68 void addSievingPrime(uint64_t prime, uint64_t segmentLow) 69 { 70 assert(segmentLow % 30 == 0); 71 72 // This hack is required because in primesieve the 8 73 // bits of each byte (of the sieve array) correspond to 74 // the offsets { 7, 11, 13, 17, 19, 23, 29, 31 }. 75 // So we are looking for: multiples > segmentLow + 6. 76 segmentLow += 6; 77 78 // calculate the first multiple (of prime) > segmentLow 79 uint64_t quotient = (segmentLow / prime) + 1; 80 quotient = std::max(prime, quotient); 81 uint64_t multiple = prime * quotient; 82 // prime not needed for sieving 83 if (multiple > stop_ || 84 multiple < segmentLow) 85 return; 86 87 // calculate the next multiple of prime that is not 88 // divisible by any of the wheel's factors 89 uint64_t nextMultipleFactor = INIT[quotient % MODULO].nextMultipleFactor; 90 uint64_t nextMultiple = prime * nextMultipleFactor; 91 if (nextMultiple > stop_ - multiple) 92 return; 93 94 nextMultiple += multiple - segmentLow; 95 uint64_t multipleIndex = nextMultiple / 30; 96 uint64_t wheelIndex = wheelOffsets_[prime % 30] + INIT[quotient % MODULO].wheelIndex; 97 storeSievingPrime(prime, multipleIndex, wheelIndex); 98 } 99 100 protected: 101 uint64_t stop_ = 0; 102 virtual ~Wheel() = default; 103 virtual void storeSievingPrime(uint64_t, uint64_t, uint64_t) = 0; 104 getMaxFactor()105 static uint64_t getMaxFactor() 106 { 107 return WHEEL[0].nextMultipleFactor; 108 } 109 110 /// Cross-off the current multiple of sievingPrime 111 /// and calculate its next multiple 112 /// unsetBit(uint8_t * sieve,uint64_t sievingPrime,uint64_t * multipleIndex,uint64_t * wheelIndex)113 static void unsetBit(uint8_t* sieve, 114 uint64_t sievingPrime, 115 uint64_t* multipleIndex, 116 uint64_t* wheelIndex) 117 { 118 sieve[*multipleIndex] &= WHEEL[*wheelIndex].unsetBit; 119 *multipleIndex += WHEEL[*wheelIndex].nextMultipleFactor * sievingPrime; 120 *multipleIndex += WHEEL[*wheelIndex].correct; 121 *wheelIndex += WHEEL[*wheelIndex].next; 122 } 123 124 private: 125 static const uint64_t wheelOffsets_[30]; 126 }; 127 128 template <int MODULO, int SIZE, const WheelElement* WHEEL, const WheelInit* INIT> 129 const uint64_t 130 Wheel<MODULO, SIZE, WHEEL, INIT>::wheelOffsets_[30] = 131 { 132 0, SIZE * 7, 0, 0, 0, 0, 133 0, SIZE * 0, 0, 0, 0, SIZE * 1, 134 0, SIZE * 2, 0, 0, 0, SIZE * 3, 135 0, SIZE * 4, 0, 0, 0, SIZE * 5, 136 0, 0, 0, 0, 0, SIZE * 6 137 }; 138 139 /// 3rd wheel, skips multiples of 2, 3 and 5 140 using Wheel30_t = Wheel<30, 8, wheel30, wheel30Init>; 141 142 /// 4th wheel, skips multiples of 2, 3, 5 and 7 143 using Wheel210_t = Wheel<210, 48, wheel210, wheel210Init>; 144 145 } // namespace 146 147 #endif 148