Verify/LC_counting_primes.test.cpp

View this file on GitHub
Last update: 2024-04-26 03:18:17+09:00
Problem: https://judge.yosupo.jp/problem/counting_primes

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/counting_primes"

#include "Template/template.hpp"

#include "Math/primesum.hpp"


ll F(ll x){return x;}

int main(){
    ll n;
    cin>>n;
    PrimeSum<ll,F> pc(n);
    cout<<pc[n]<<'\n';
    return 0;
}

#line 1 "Verify/LC_counting_primes.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/counting_primes"

#line 1 "Template/template.hpp"
#include <bits/stdc++.h>

using namespace std;

#define rep(i, a, b) for (int i = (int)(a); i < (int)(b); i++)
#define rrep(i, a, b) for (int i = (int)(b-1); i >= (int)(a); i--)
#define ALL(v) (v).begin(), (v).end()
#define UNIQUE(v) sort(ALL(v)), (v).erase(unique(ALL(v)), (v).end())
#define SZ(v) (int)v.size()
#define MIN(v) *min_element(ALL(v))
#define MAX(v) *max_element(ALL(v))
#define LB(v, x) int(lower_bound(ALL(v), (x)) - (v).begin())
#define UB(v, x) int(upper_bound(ALL(v), (x)) - (v).begin())

using uint = unsigned int;
using ll = long long int;
using ull = unsigned long long;
using i128 = __int128_t;
using u128 = __uint128_t;
const int inf = 0x3fffffff;
const ll INF = 0x1fffffffffffffff;

template <typename T> inline bool chmax(T &a, T b) {
    if (a < b) {
        a = b;
        return 1;
    }
    return 0;
}
template <typename T> inline bool chmin(T &a, T b) {
    if (a > b) {
        a = b;
        return 1;
    }
    return 0;
}
template <typename T, typename U> T ceil(T x, U y) {
    assert(y != 0);
    if (y < 0)
        x = -x, y = -y;
    return (x > 0 ? (x + y - 1) / y : x / y);
}
template <typename T, typename U> T floor(T x, U y) {
    assert(y != 0);
    if (y < 0)
        x = -x, y = -y;
    return (x > 0 ? x / y : (x - y + 1) / y);
}
template <typename T> int popcnt(T x) {
    return __builtin_popcountll(x);
}
template <typename T> int topbit(T x) {
    return (x == 0 ? -1 : 63 - __builtin_clzll(x));
}
template <typename T> int lowbit(T x) {
    return (x == 0 ? -1 : __builtin_ctzll(x));
}

#ifdef LOCAL
#define show(...) _show(0, #__VA_ARGS__, __VA_ARGS__)
#else
#define show(...) true
#endif
template <typename T> void _show(int i, T name) {
    cerr << '\n';
}
template <typename T1, typename T2, typename... T3>
void _show(int i, const T1 &a, const T2 &b, const T3 &...c) {
    for (; a[i] != ',' && a[i] != '\0'; i++)
        cerr << a[i];
    cerr << ":" << b << " ";
    _show(i + 1, a, c...);
}
template <class T, class U>
ostream &operator<<(ostream &os, const pair<T, U> &p) {
    os << "P(" << p.first << ", " << p.second << ")";
    return os;
}
template <typename T, template <class> class C>
ostream &operator<<(ostream &os, const C<T> &v) {
    os << "[";
    for (auto d : v)
        os << d << ", ";
    os << "]";
    return os;
}
#line 2 "Math/sieve.hpp"

template<int L=50101010>vector<int> sieve(int N){
    bitset<L> isp;
    int n,sq=ceil(sqrt(N));
    for(int z=1;z<=5;z+=4){
        for(int y=z;y<=sq;y+=6){
            for(int x=1;x<=sq and (n=4*x*x+y*y)<=N;++x){
                isp[n].flip();
            }
            for(int x=y+1;x<=sq and (n=3*x*x-y*y)<=N;x+=2){
                isp[n].flip();
            }
        }
    }
    for(int z=2;z<=4;z+=2){
        for(int y=z;y<=sq;y+=6){
            for (int x=1;x<=sq and (n=3*x*x+y*y)<=N;x+=2){
                isp[n].flip();
            }
            for(int x=y+1;x<=sq and (n=3*x*x-y*y)<=N;x+=2){
                isp[n].flip();
            }
        }
    }
    for(int y=3;y<=sq;y+=6){
        for(int z=1;z<=2;++z){
            for(int x=z;x<=sq and (n=4*x*x+y*y)<=N;x+=3){
                isp[n].flip();
            }
        }
    }
    for(int n=5;n<=sq;++n)if(isp[n]){
        for(int k=n*n;k<=N;k+=n*n){
            isp[k]=false;
        }
    }
    isp[2]=isp[3]=true;

    vector<int> ret;
    for(int i=2;i<=N;i++)if(isp[i]){
        ret.push_back(i);
    }
    return ret;
}

/**
 * @brief Prime Sieve
 */
#line 3 "Math/primesum.hpp"

template<typename T,T (*F)(ll)>struct PrimeSum{
    ll N,SQ;
    vector<T> lo,hi;
    PrimeSum(ll n=0):N(n),SQ(sqrtl(N)),lo(SQ+1),hi(SQ+1){
        rep(i,1,SQ+1){
            lo[i]=F(i)-1;
            hi[i]=F(N/i)-1;
        }
        auto ps=sieve(SQ);
        for(auto& p:ps){
            ll q=ll(p)*p;
            if(q>N)break;
            T sub=lo[p-1],fp=lo[p]-lo[p-1];
            ll L=min(SQ,N/q),M=SQ/p;
            rep(i,1,M+1)hi[i]-=fp*(hi[i*p]-sub);
            rep(i,M+1,L+1)hi[i]-=fp*(lo[double(N)/(i*p)]-sub);
            for(int i=SQ;i>=q;i--)lo[i]-=fp*(lo[double(i)/p]-sub);
        }
    }
    T operator[](ll x) {
        return (x<=SQ?lo[x]:hi[N/x]);
    }
};

/**
 * @brief Prime Sum
 * @docs docs/primesum.md
 */
#line 5 "Verify/LC_counting_primes.test.cpp"

ll F(ll x){return x;}

int main(){
    ll n;
    cin>>n;
    PrimeSum<ll,F> pc(n);
    cout<<pc[n]<<'\n';
    return 0;
}