Count Square Free
(number-theory/count-square-free.hpp)

• View this file on GitHub
• Last update: 2026-07-25 02:01:37+09:00
• Include: #include "number-theory/count-square-free.hpp"

$N$ 以下の無平方数の個数を $O(N^{2/5}\log\log N)$ 時間で求める．

アルゴリズム

参考：Counting square free numbers - Blog of smsxgz

$N$ 以下の無平方数の個数を $S(N)$ とする．

$~O(\sqrt{N})$

素因数についての包除で次を得る． \(S(N)=\sum_{d}\mu(d)\left\lfloor\frac{N}{d^2}\right\rfloor\) $n$ 以下に対するMobius関数の列挙は篩で $~O(n)$ なのでこれは $~O(\sqrt{N})$ で計算できる．

高速化

$D$ を十分大きくとると，$\left\lfloor\frac{N}{d^2}\right\rfloor$ が $d\gt D$ の範囲であまり変化しなくなる．

$S_1(N)=\sum_{1\leq d\leq D}\mu(d)\left\lfloor\frac{N}{d^2}\right\rfloor,S_2(N)=\sum_{d\gt D}\mu(d)\left\lfloor\frac{N}{d^2}\right\rfloor$ とする．

$S_1$ は直接計算もできるが，$S_2$ でも使えるので $\mu$ 関数の値を列挙してそこから計算するようにする． $\mu(1),\dots,\mu(D)$ の列挙は篩を考えると空間 $O(D)$，時間 $O(D\log\log D)$ で計算できる．

$S_2$ を考察． \(\begin{align*} S_2(N) &=\sum_{d\gt D}\mu(d)\sum_i1_{i=\left\lfloor\frac{N}{d^2}\right\rfloor}i\\ &=\sum_i i\sum_{d\gt D}\mu(d)1_{i\leq\frac{N}{d^2}\lt i+1}\\ &=\sum_i i\sum_{d\gt D,\left\lfloor\sqrt{\frac{N}{i+1}}\right\rfloor\lt d\leq\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor}\mu(d) \end{align*}\)

$x_i=\sqrt{\frac{N}{i}}$ とおく．正整数 $I$ をとり $D=\lfloor x_I\rfloor$ とすれば，Mertens関数 $M(x)=\sum_{i=1}^{\lfloor x\rfloor}\mu(i)$ を用いて次のように表せる． \(S_2(N) =\sum_{i=1}^{I-1}i(M(x_i)-M(x_{i+1})) =\sum_{i=1}^{I-1}M(x_i)-(I-1)M(x_I)\)

ここでメビウス反転公式から $\sum_{n=1}^{x}M\left(\frac{x}{n}\right)=1$ である． 変形すれば $M(x)=1-\sum_{n=2}^{x}M\left(\frac{x}{n}\right)$．

特に $M(x/n),n\geq 2$ の値がわかっているとき $M(x)$ が $O(\sqrt{x})$ で計算できる．

$M(x_I),M(x_{I-1}),\dots,M(x_1)$ をこの順に計算していくことを考える．

$M(x_k)$ を計算するには $M(x_k/i),2\leq i\leq x_k$ の値が必要．

• $x_k/i\lt D$ のとき：$S_1$ の計算で $\mu(j),j\leq D$ を列挙しているので累積和をとればよい．
• $x_k/i\geq D$ のとき：$x_k/i=x_{ki^2}\geq D=x_I$ より $ki^2\leq I$ であり $M(x_{ki^2})$ の値はすでに求めてある．

$M(x_k/i)$ は整数部分 $\lfloor x_k/i\rfloor$ で決まるので $M(x_k)$ は $O(\sqrt{x_k})$ で計算できる．

計算量を解析する．

$S_1$ の計算量は $O(D\log\log D)=O(\sqrt{N/I}\log\log(N/I))$．

$S_2$ の計算量は，次が成り立つから $O(N^{1/4}I^{3/4})$． \(\sum_{k=1}^{I}\sqrt{x_k} =\sum_{k=1}^{I}\frac{N^{1/4}}{k^{1/4}} =O(N^{1/4}I^{3/4})\)

ここで $I=N^\alpha$ とおけば上の計算量はそれぞれ $O(N^{(1-\alpha)/2}\log\log N)$ および $O(N^{1/4+3/4\alpha})$ になる． $\alpha=1/5$ とすれば全体の計算量は $O(N^{2/5}\log\log N)$ となった（このとき $D=N^{2/5}$）．

Depends on

• math/util.hpp
• Mobius Function (number-theory/mobius-function.hpp)
• 素数篩 (number-theory/prime-sieve.hpp)

Verified with

• verify/number-theory/LC_counting_squarefrees.test.cpp

Code

#pragma once

#include "number-theory/mobius-function.hpp"
#include "math/util.hpp"

long long CountSquareFree(long long n) {
  using ll = long long;
  if (n <= 0) return 0;
  ll ret = 0;
  if (n < 100) {
    auto mu = MobiusFunction::table((int)n);
    for (int i = 1; i <= n; i++) ret += n / i / i * mu[i];
  } else {
    ll l = Math::floor_root(n, 5);
    ll d = Math::floor_root(n / l, 2);
    auto mu = MobiusFunction::table(d);
    for (int i = 1; i <= d; i++) ret += n / i / i * mu[i];
    for (int i = 1; i <= d; i++) mu[i] += mu[i - 1];  // mu 自体はもういらないので使い回し
    vector<ll> a(l + 1);
    for (int k = l; k > 0; k--) {
      ll x = Math::floor_root(n / k, 2);
      ll sq = Math::floor_root(x, 2);
      ll m = x / (sq + 1);
      ll v = 1;
      for (ll i = 1; i <= m; i++) v -= (x / i - x / (i + 1)) * mu[i];
      for (ll i = 2; i <= sq; i++) v -= x / i <= d ? mu[x / i] : a[k * i * i];
      a[k] = v;
      if (k < l) ret += a[k];
    }
    ret -= (l - 1) * a[l];
  }
  return ret;
}

/**
 * @brief Count Square Free
 * @docs docs/number-theory/count-square-free.md
 */

#line 2 "number-theory/count-square-free.hpp"

#line 2 "number-theory/mobius-function.hpp"

#line 2 "math/util.hpp"

namespace Math {
template <class T>
T safe_mod(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  a %= b;
  return a >= 0 ? a : a + b;
}
template <class T>
T floor(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  return a >= 0 ? a / b : (a + 1) / b - 1;
}
template <class T>
T ceil(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  return a > 0 ? (a - 1) / b + 1 : a / b;
}
long long isqrt(long long n) {
  if (n <= 0) return 0;
  long long x = sqrt(n);
  while ((__int128)(x + 1) * (x + 1) <= n) x++;
  while ((__int128)x * x > n) x--;
  return x;
}
long long floor_root(long long n, int k) {
  assert(n >= 0);
  if (n == 0) return 0;
  assert(k >= 1);
  if (k == 1) return n;
  if (k > 64) return 1;
  long long x = round(pow((long double)n, 1.0L / k));
  auto check = [&](long long a) {
    if (a <= 0) return true;
    __int128_t p = 1;
    for (int i = 0; i < k; ++i)
      if ((p *= a) > n) return false;
    return true;
  };
  while (check(x + 1)) x++;
  while (!check(x)) x--;
  return x;
}
unsigned long long floor_root_unsigned(unsigned long long n, int k) {
  assert(k >= 1);
  if (n <= 1 || k == 1) return n;
  if (k >= 64) return 1;
  int bits = (64 + k - 1) / k;
  unsigned long long ok = 1, ng = min(n, 1ULL << bits);
  auto check = [&](unsigned long long a) {
    __uint128_t p = 1;
    for (int i = 0; i < k; i++) {
      p *= a;
      if (p > n) return false;
    }
    return true;
  };
  while (ok + 1 < ng) {
    unsigned long long mid = ok + (ng - ok) / 2;
    (check(mid) ? ok : ng) = mid;
  }
  return ok;
}
// return g=gcd(a,b)
// a*x+b*y=g
// - b!=0 -> 0<=x<|b|/g
// - b=0  -> ax=g
template <class T>
T ext_gcd(T a, T b, T& x, T& y) {
  T a0 = a, b0 = b;
  bool sgn_a = a < 0, sgn_b = b < 0;
  if (sgn_a) a = -a;
  if (sgn_b) b = -b;
  if (b == 0) {
    x = sgn_a ? -1 : 1;
    y = 0;
    return a;
  }
  T x00 = 1, x01 = 0, x10 = 0, x11 = 1;
  while (b != 0) {
    T q = a / b, r = a - b * q;
    x00 -= q * x01;
    x10 -= q * x11;
    swap(x00, x01);
    swap(x10, x11);
    a = b, b = r;
  }
  x = x00, y = x10;
  if (sgn_a) x = -x;
  if (sgn_b) y = -y;
  if (b0 != 0) {
    a0 /= a, b0 /= a;
    if (b0 < 0) a0 = -a0, b0 = -b0;
    T q = x >= 0 ? x / b0 : (x + 1) / b0 - 1;
    x -= b0 * q;
    y += a0 * q;
  }
  return a;
}
constexpr long long inv_mod(long long x, long long m) {
  x %= m;
  if (x < 0) x += m;
  long long a = m, b = x;
  long long y0 = 0, y1 = 1;
  while (b > 0) {
    long long q = a / b;
    swap(a -= q * b, b);
    swap(y0 -= q * y1, y1);
  }
  if (y0 < 0) y0 += m / a;
  return y0;
}
long long pow_mod(long long x, long long n, long long m) {
  if (m == 1) return 0;
  x = (x % m + m) % m;
  long long y = 1;
  while (n) {
    if (n & 1) y = y * x % m;
    x = x * x % m;
    n >>= 1;
  }
  return y;
}
constexpr long long pow_mod_constexpr(long long x, long long n, int m) {
  if (m == 1) return 0;
  unsigned int _m = (unsigned int)(m);
  unsigned long long r = 1;
  unsigned long long y = x % m;
  if (y >= m) y += m;
  while (n) {
    if (n & 1) r = (r * y) % _m;
    y = (y * y) % _m;
    n >>= 1;
  }
  return r;
}
constexpr bool is_prime_constexpr(int n) {
  if (n <= 1) return false;
  if (n == 2 || n == 7 || n == 61) return true;
  if (n % 2 == 0) return false;
  long long d = n - 1;
  while (d % 2 == 0) d /= 2;
  constexpr long long bases[3] = {2, 7, 61};
  for (long long a : bases) {
    long long t = d;
    long long y = pow_mod_constexpr(a, t, n);
    while (t != n - 1 && y != 1 && y != n - 1) {
      y = y * y % n;
      t <<= 1;
    }
    if (y != n - 1 && t % 2 == 0) {
      return false;
    }
  }
  return true;
}
template <int n>
constexpr bool is_prime = is_prime_constexpr(n);
};  // namespace Math
#line 2 "number-theory/prime-sieve.hpp"

namespace PrimeSieve {
using ll = long long;

vector<int> lpf(int n) {
  assert(n >= 0);
  vector<int> ret(n + 1);
  for (size_t i = 0; i < ret.size(); i++) ret[i] = (int)i;
  for (int p = 2; (ll)p * p <= n; p++) {
    if (ret[p] != p) continue;
    for (ll x = (ll)p * p;; x += p) {
      if (ret[x] == x) ret[x] = p;
      if (n - x < p) break;
    }
  }
  return ret;
}

vector<int> table(int n) {
  assert(n >= 0);
  vector<bool> composite(n + 1, false);
  for (int p = 2; (ll)p * p <= n; p += (p & 1) + 1) {
    if (composite[p]) continue;
    for (ll x = (ll)p * p;; x += p) {
      composite[x] = true;
      if (n - x < p) break;
    }
  }
  vector<int> ps;
  for (int p = 2; p <= n;) {
    if (!composite[p]) ps.push_back(p);
    int step = (p & 1) + 1;
    if (n - p < step) break;
    p += step;
  }
  return ps;
}

vector<vector<pair<ll, int>>> factorize(int n) {
  assert(n >= 0);
  vector<vector<pair<ll, int>>> factors(n + 1);
  auto lp = lpf(n);
  for (int x = 2; x <= n;) {
    int y = x;
    while (y > 1) {
      int p = lp[y], e = 0;
      while (y % p == 0) y /= p, e++;
      factors[x].emplace_back(p, e);
    }
    if (x == n) break;
    x++;
  }
  return factors;
}
};  // namespace PrimeSieve

/**
 * @brief 素数篩
 * @docs docs/number-theory/prime-sieve.md
 */
#line 5 "number-theory/mobius-function.hpp"

namespace MobiusFunction {
using ll = long long;

vector<int> table(int n) {
  vector<int> mu(n + 1, 1);
  mu[0] = 0;
  auto lpf = PrimeSieve::lpf(n);
  for (int x = 2; x <= n; x++) {
    int p = lpf[x];
    if (x / p % p == 0)
      mu[x] = 0;
    else
      mu[x] = -mu[x / p];
  }
  return mu;
}

ll sum(ll n) {
  if (n <= 0) return 0;

  ll k = ceil(pow((long double)n, 2.0L / 3.0L));
  __int128 n2 = (__int128)n * n;
  auto enough = [&](ll x) { return (__int128)x * x * x >= n2; };
  while (k > 1 && enough(k - 1)) k--;
  while (!enough(k)) k++;
  assert(k <= numeric_limits<int>::max());

  int lim = (int)k;
  vector<ll> small(lim + 1);
  auto mu = table(lim);
  for (int i = 1; i <= lim; i++) small[i] = small[i - 1] + mu[i];

  ll len = n / k + (n % k != 0);
  vector<ll> large(len + 1);
  for (ll i = len; i >= 1; i--) {
    ll x = n / i;
    ll m = Math::isqrt(x);
    ll v = 1;
    for (ll j = 2; j <= m; j++) v -= i * j <= len ? large[i * j] : small[x / j];
    for (ll j = 1; j <= m; j++) v -= (x / j - m) * (small[j] - small[j - 1]);
    large[i] = v;
  }
  return large[1];
}
};  // namespace MobiusFunction

/**
 * @brief Mobius Function
 * @docs docs/number-theory/mobius-function.md
 */
#line 5 "number-theory/count-square-free.hpp"

long long CountSquareFree(long long n) {
  using ll = long long;
  if (n <= 0) return 0;
  ll ret = 0;
  if (n < 100) {
    auto mu = MobiusFunction::table((int)n);
    for (int i = 1; i <= n; i++) ret += n / i / i * mu[i];
  } else {
    ll l = Math::floor_root(n, 5);
    ll d = Math::floor_root(n / l, 2);
    auto mu = MobiusFunction::table(d);
    for (int i = 1; i <= d; i++) ret += n / i / i * mu[i];
    for (int i = 1; i <= d; i++) mu[i] += mu[i - 1];  // mu 自体はもういらないので使い回し
    vector<ll> a(l + 1);
    for (int k = l; k > 0; k--) {
      ll x = Math::floor_root(n / k, 2);
      ll sq = Math::floor_root(x, 2);
      ll m = x / (sq + 1);
      ll v = 1;
      for (ll i = 1; i <= m; i++) v -= (x / i - x / (i + 1)) * mu[i];
      for (ll i = 2; i <= sq; i++) v -= x / i <= d ? mu[x / i] : a[k * i * i];
      a[k] = v;
      if (k < l) ret += a[k];
    }
    ret -= (l - 1) * a[l];
  }
  return ret;
}

/**
 * @brief Count Square Free
 * @docs docs/number-theory/count-square-free.md
 */

Back to top page