Count Square-free integers
(Math/countsquarefree.hpp)

• View this file on GitHub
• Last update: 2024-04-26 03:18:17+09:00
• Include: #include "Math/countsquarefree.hpp"

Depends on

• Kth Root(Integer) (Math/kthroot.hpp)

Code

#pragma once
#include "Math/kthroot.hpp"

ll CountSquarefree(ll n) {
    if (n <= 3)
        return n;
    const int I = Kthroot(n, 5);
    const int D = sqrtl(n / I);
    vector<int> prime(D + 1, 1), mobius(D + 1);
    mobius[1] = 1;
    rep(i, 2, D + 1) if (prime[i]) {
        for (int j = D / i; j >= 1; j--) {
            if (j >= i and prime[j])
                prime[i * j] = 0;
            if (mobius[j])
                mobius[i * j] = -mobius[j];
        }
    }

    ll ret = 0;
    rep(i, 1, D + 1) ret += mobius[i] * (n / (ll(i) * i));

    auto mertens = mobius;
    rep(i, 0, D) mertens[i + 1] += mertens[i];

    vector<ll> Mxi(I);
    for (int i = I - 1; i >= 1; i--) {
        const int xi = sqrtl(n / i);
        const int sqxi = sqrtl(xi);
        Mxi[i] = 1;
        rep(j, 2, sqxi + 1) {
            if (xi / j <= D)
                Mxi[i] -= mertens[xi / j];
            else
                Mxi[i] -= Mxi[i * j * j];
        }
        rep(j, 1, xi / (sqxi + 1) + 1) Mxi[i] -=
            (xi / j - xi / (j + 1)) * mertens[j];
        ret += Mxi[i];
    }
    ret -= ll(I - 1) * mertens[D];
    return ret;
}

/**
 * @brief Count Square-free integers
 */

#line 2 "Math/kthroot.hpp"

uint64_t Kthroot(uint64_t k, uint64_t a) {
    assert(k >= 1);
    if (a == 0 || a == 1 || k == 1) return a;
    if (k >= 64) return 1;
    if (k == 2) return sqrtl(a);
    if (a == uint64_t(-1)) --a;
    struct S {
        uint64_t v;
        S& operator*=(const S& o) {
            v = v <= uint64_t(-1) / o.v ? v * o.v : uint64_t(-1);
            return *this;
        }
    };
    auto power = [&](S x, ll n) -> S {
        S v{1};
        while (n) {
            if (n & 1) v *= x;
            x *= x;
            n /= 2;
        }
        return v;
    };
    uint64_t res = pow(a, nextafter(1 / double(k), 0));
    while (power(S{res + 1}, k).v <= a) ++res;
    return res;
}

/**
 * @brief Kth Root(Integer)
*/
#line 3 "Math/countsquarefree.hpp"

ll CountSquarefree(ll n) {
    if (n <= 3)
        return n;
    const int I = Kthroot(n, 5);
    const int D = sqrtl(n / I);
    vector<int> prime(D + 1, 1), mobius(D + 1);
    mobius[1] = 1;
    rep(i, 2, D + 1) if (prime[i]) {
        for (int j = D / i; j >= 1; j--) {
            if (j >= i and prime[j])
                prime[i * j] = 0;
            if (mobius[j])
                mobius[i * j] = -mobius[j];
        }
    }

    ll ret = 0;
    rep(i, 1, D + 1) ret += mobius[i] * (n / (ll(i) * i));

    auto mertens = mobius;
    rep(i, 0, D) mertens[i + 1] += mertens[i];

    vector<ll> Mxi(I);
    for (int i = I - 1; i >= 1; i--) {
        const int xi = sqrtl(n / i);
        const int sqxi = sqrtl(xi);
        Mxi[i] = 1;
        rep(j, 2, sqxi + 1) {
            if (xi / j <= D)
                Mxi[i] -= mertens[xi / j];
            else
                Mxi[i] -= Mxi[i * j * j];
        }
        rep(j, 1, xi / (sqxi + 1) + 1) Mxi[i] -=
            (xi / j - xi / (j + 1)) * mertens[j];
        ret += Mxi[i];
    }
    ret -= ll(I - 1) * mertens[D];
    return ret;
}

/**
 * @brief Count Square-free integers
 */

Back to top page