-CP_ALGO_BIT_PRAGMA_PUSHnamespace cp_algo::math{using cp_algo::structures::dynamic_bit_array;constexpr size_t base_threshold=1<<21;constexpr auto to_ord(auto x){return x/2;}constexpr auto to_val(auto x){return 2*x+1;}const auto base_prime_bits=[](){dynamic_bit_array prime(base_threshold);prime.set_all();prime.reset(to_ord(1));for(size_t i=3;to_ord(i*i)<base_threshold;i+=2){if(prime[to_ord(i)]){for(size_t j=i*i;to_ord(j)<base_threshold;j+=2*i){prime.reset(to_ord(j));}}}return prime;}();const auto base_primes=[](){cp_algo::big_vector<uint32_t>primes;for(uint32_t i=3;to_ord(i)<base_threshold;i+=2){if(base_prime_bits[to_ord(i)]){primes.push_back(i);}}return primes;}();constexpr size_t sqrt_threshold=1<<16;const auto sqrt_primes=std::span(begin(base_primes),std::ranges::upper_bound(base_primes,sqrt_threshold));constexpr size_t max_wheel_size=std::min<size_t>(base_threshold,1<<21);struct wheel_t{dynamic_bit_array mask;uint32_t product;};auto make_wheel(std::vector<uint32_t>primes,uint32_t product){assert(product%(2*dynamic_bit_array::width)==0);wheel_t wheel;wheel.product=product;wheel.mask.resize(product/2);wheel.mask.set_all();for(auto p:primes){for(size_t j=to_ord(p);j<wheel.mask.size();j+=p){wheel.mask.reset(j);}}return wheel;}auto medium_primes=sqrt_primes;const auto wheels=[](){uint32_t product=2*dynamic_bit_array::width;std::vector<uint32_t>current;std::vector<wheel_t>wheels;for(size_t i=0;i<size(sqrt_primes);i++){uint32_t p=sqrt_primes[i];if(product*p>max_wheel_size){if(size(current)==1){medium_primes=sqrt_primes.subspan(i-1);return wheels;}wheels.push_back(make_wheel(current,product));current={p};product=2*dynamic_bit_array::width*p;}else{current.push_back(p);product*=p;}}assert(false);}();void sieve_dense(auto&prime,uint32_t l,uint32_t r,wheel_t const&wheel){if(l>=r)return;const auto width=(uint32_t)dynamic_bit_array::width;uint32_t wl=l/width;uint32_t wr=(r+width-1)/width;uint32_t N=(uint32_t)wheel.mask.words;for(uint32_t i=wl;i<wr;i+=N){auto block=std::min(N,wr-i);for(uint32_t j=0;j<block;j++){prime.word(i+j)&=wheel.mask.word(j);}}}constexpr auto add210=[](){std::array<uint8_t,210>add;auto good=[&](int x){return x%2&&x%3&&x%5&&x%7;};for(int i=0;i<210;i++){add[i]=1;while(!good(i+add[i])){add[i]++;}}return add;}();constexpr uint8_t gap210[]={5,1,2,1,2,3,1,3,2,1,2,3,3,1,3,2,1,3,2,3,4,2,1,2,1,2,4,3,2,3,1,2,3,1,3,3,2,1,2,3,1,3,2,1,2,1,5,1};constexpr auto state210=[](){std::array<uint8_t,210>state;int idx=0;for(int i=0;i<210;i++){if(i%2&&i%3&&i%5&&i%7){state[i]=uint8_t(idx++);}else{state[i]=-1;}}return state;}();template<class BitArray>void sieve210(BitArray&prime,uint32_t l,uint32_t r,uint32_t p,uint8_t state){if(l>=r)return;while(l<r){prime.reset(l);l+=p*gap210[state];state=state==47?0:state+1;}}dynamic_bit_array odd_sieve(uint32_t N){N++;dynamic_bit_array prime(N/2);prime.set_all();for(size_t i=0;i<std::min(prime.words,base_prime_bits.words);i++){prime.word(i)=base_prime_bits.word(i);}cp_algo::checkpoint("init");static constexpr uint32_t dense_block=1<<24;for(uint32_t start=base_threshold;start<N;start+=dense_block){uint32_t r=std::min(start+dense_block,N);for(auto const&wheel:wheels){auto l=start/wheel.product*wheel.product;sieve_dense(prime,to_ord(l),to_ord(r),wheel);}}cp_algo::checkpoint("dense sieve");static constexpr uint32_t sparse_block=1<<22;for(uint32_t start=base_threshold;start<N;start+=sparse_block){uint32_t r=std::min(start+sparse_block,N);for(auto p:medium_primes){if(p*p>=r)break;auto k=std::max(start/p,p);if(state210[k%210]==0xFF){k+=add210[k%210];}sieve210(prime,to_ord(k*p),to_ord(r),p,state210[k%210]);}}cp_algo::checkpoint("sparse sieve");return prime;}}
0 commit comments