Cod sursa(job #2929181)

Utilizator lucametehauDart Monkey lucametehau Data 25 octombrie 2022 00:45:06
Problema Ciurul lui Eratosthenes Scor 100
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 7.89 kb
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <vector>
#include <map>
#include <unordered_map>
#include <set>
#include <cstring>
#include <chrono>
#include <cassert>
#include <bitset>
#include <stack>
#include <queue>
#include <iomanip>
#include <omp.h>
#include <random>
//#include <ext/pb_ds/assoc_container.hpp>

#pragma GCC optimize("Ofast")
//#pragma GCC optimize("Ofast,unroll-loops")
//#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native")
#define x first
#define y second
#define ld long double
#define ll long long
#define ull unsigned long long
#define us unsigned short
#define lsb(x) ((x) & (-(x)))
#define pii pair <int, int>

using namespace std;

#ifdef _MSC_VER
#  include <intrin.h>
#  define __builtin_popcount __popcnt
#  define __builtin_popcountll __popcnt64
#endif
//using namespace __gnu_pbds;

mt19937 gen(time(0));
uniform_int_distribution <int64_t> rng;

ifstream in("ciur.in");
ofstream out("ciur.out");

const int N = (int)1e4; // precalculated values
const int N4 = (int)1e3;
const ll INF = (ll)1e18;

ll n;

int hashHits, hashHitsP;

int cntPrimes;

bitset <N + 5> viz;
int prime[10005];
int pi[N + 5];

int lp[N4 + 5], mu[N4 + 5];

ll get_pi(ll n);
ll P(int niv, ll n, int a);

void sieve() {
    for (int i = 4; i <= N; i += 2)
        viz[i] = 1;
    for (int i = 3; i * i <= N; i++) {
        if (!viz[i]) {
            for (int j = i * i; j <= N; j += 2 * i)
                viz[j] = 1;
        }
    }

    for (int i = 1; i <= N4; i++)
        lp[i] = i, mu[i] = 1;

    for (int i = 2; i <= N4; i++) {
        if (lp[i] == i) {
            for (int j = i; j <= N4; j += i)
                lp[j] = min(lp[j], i), mu[j] *= -1;
            if (1LL * i * i > N4)
                continue;
            for (int j = i * i; j <= N4; j += i * i)
                mu[j] = 0;
        }
    }

    cntPrimes = 0;

    for (int i = 2; i <= N; i++) {
        pi[i] = pi[i - 1];
        if (!viz[i])
            prime[++cntPrimes] = i, pi[i]++;
    }
}

struct SimpleDS {
    static const int SIZE = 32;
    static const int BUCKETS = 10000 / SIZE;
    uint32_t v[BUCKETS];
    uint32_t allMask[SIZE];

    SimpleDS() {
        allMask[0] = 1;
        for (int i = 1; i < SIZE; i++)
            allMask[i] = (allMask[i - 1] << 1) | 1;
        for (int i = 0; i < BUCKETS; i++)
            v[i] = allMask[SIZE - 1];
    }

    void init(int lg) {
        for (int i = 0; i <= lg / SIZE; i++)
            v[i] = allMask[SIZE - 1];
    }

    void setZero(int ind) {
        int b = ind / SIZE, p = ind % SIZE;
        v[b] &= ~(1u << p);
    }

    int count(int poz) {
        int ans = 0;
        for (int i = 0; i < poz / SIZE; i++)
            ans += __builtin_popcount(v[i]);
        ans += __builtin_popcount(v[poz / SIZE] & allMask[poz % SIZE]);
        return ans - 1;
    }

    int count(int l, int r) {
        if (l > r)
            return 0;
        int bl = l / SIZE, br = r / SIZE;
        if (bl == br)
            return __builtin_popcount(v[bl] & (l % SIZE == 0 ? allMask[SIZE - 1] : ~allMask[l % SIZE - 1]) & allMask[r % SIZE]);
        int ans = __builtin_popcount(v[bl] & (l % SIZE == 0 ? allMask[SIZE - 1] : ~allMask[l % SIZE - 1]));
        for (int i = bl + 1; i < br; i++)
            ans += __builtin_popcount(v[i]);
        ans += __builtin_popcount(v[br] & allMask[r % SIZE]);

        return ans;
    }
} simple_ds;

int f(int n, int m) {
    if (m == 0)
        return n;
    if (n == 0)
        return 0;
    return f(n, m - 1) - f(n / prime[m], m - 1);
}

ll phi(ll n, int m) {
    int a = pi[m];
    ll s1 = 0, s2 = 0;

    m = prime[a];

    for (int i = 1; i <= m; i++) {
        s1 += n / i * mu[i];
    }

    int bucket_size = sqrt(n / m);
    vector <int> phi(a + 1);
    vector <ll> lst(a + 1);

    fill(phi.begin(), phi.end(), 0);
    for (int i = 1; i <= a; i++)
        lst[i] = prime[i];

    int b = pi[int(sqrt(m))];

    for (ll start = 1; start <= n / m; start += bucket_size) {
        ll end = min(start + bucket_size, n / m + 1);

        int lg = end - start;

        simple_ds.init(lg);

        for (int i = 1; i <= b; i++) {
            ll l = max<ll>(m / prime[i], n / (end * prime[i]));
            ll r = min<ll>(m, n / (start * prime[i]));

            if (prime[i] >= r) {
                break;
            }

            int ind = 1;
            int value = 0;

            for (ll x = r; x > l; x--) {
                if (mu[x] && lp[x] > prime[i]) {
                    ll val = n / (1LL * x * prime[i]);

                    value += simple_ds.count(ind, val - start + 1);
                    ind = val - start + 2;

                    s2 -= (phi[i] + value) * mu[x];
                }
            }

            //cout << cnt << "\n";

            phi[i] += value + simple_ds.count(ind, lg);

            for (; lst[i] < end; lst[i] += (1 + (i > 1)) * prime[i]) {
                simple_ds.setZero(lst[i] - start + 1);
            }
        }

        for (int i = b + 1; i < a; i++) {
            int ind = pi[min<ll>(m, n / (start * prime[i]))];
            if (prime[i] >= prime[ind])
                break;

            ll l = max<ll>(prime[i], n / (end * prime[i]));
            int id = 1, value = 0;
                
            for (int x = ind; prime[x] > l; x--) {
                ll val = n / (1LL * prime[x] * prime[i]);

                value += simple_ds.count(id, val - start + 1);
                id = val - start + 2;

                s2 += phi[i] + value;
            }


            phi[i] += value + simple_ds.count(id, lg);

            for (; lst[i] < end; lst[i] += (1 + (i > 1)) * prime[i]) {
                simple_ds.setZero(lst[i] - start + 1);
            }
        }
    }

    return s1 + s2;
}

bitset <1000005> ciur;

ll cnt_primes(ll l, ll r) {
    //simple_ds.init(r - l + 1);
    for (int i = 0; i <= r - l; i++)
        ciur[i] = 1;

    ll cnt = r - l + 1;

    for (int i = 1; 1LL * prime[i] * prime[i] <= r; i++) {
        int p = prime[i];
        for (ll j = max<ll>(1LL * p * p, (l + p - 1) / p * p); j <= r; j += p) {
            cnt -= ciur[j - l];
            ciur[j - l] = 0;
        }
    }

    return cnt;
}

ll P2(ll n, int m) {
    int b = get_pi(sqrt(n));
    ld t1 = clock();
    ll ans = 0;
    int a = pi[m];

    ll sum_pi = get_pi(n / prime[b]);
    ans += sum_pi;

    for (int i = b - 1; i > a; i--) {
        sum_pi += cnt_primes(n / prime[i + 1] + 1, n / prime[i]);

        ans += sum_pi;
    }

    //cout << (clock() - t1) / CLOCKS_PER_SEC << "s for P2\n";

    ans -= 1LL * (b - a) * (b + a - 1) / 2;

    return ans;
}

ll get_pi(ll n) { // number of primes <= n
    if (n <= N)
        return pi[n];
    int root = round(pow(n, 1.0 / 3)) * log(log(n)), a = pi[root];

    ll ans = phi(n, root) + a - 1;
    //ll ans = f(n, a) + a - 1;
    ans -= P2(n, root);
    return ans;
}

ll get_prime(ll n) { // n-th prime number
    if (n <= cntPrimes)
        return prime[n];
    ll st = n * log(n) + n * log(log(n)) - n, dr = n * log(n) + n * log(log(n)), mid;
    while (st <= dr) {
        mid = (st + dr) >> 1;
        //reset();
        if (get_pi(mid) < n)
            st = mid + 1;
        else
            dr = mid - 1;
    }
    return st;
}

int main() {
    in >> n;
    ld start = clock();
    sieve();
    //cout << (clock() - start) / CLOCKS_PER_SEC << "s time taken for precalc\n";
    //precalc();
    out << get_pi(n) << "\n";
    ld end = clock();
    //cout << (end - start) / CLOCKS_PER_SEC << "s time taken\n";

    //cout << "hits: " << hashHits << "\n";
    //cout << "hits for p: " << hashHitsP << "\n";
    //cout << small_calls << "/" << calls << " " << 100.0 * small_calls / calls << "% small calls\n";
    return 0;
}