Skip to content

squarefree_count.hpp

SECTIONMath INCLUDEnoya/squarefree_count.hpp

计算不超过 \(n\) 的无平方因子整数个数，适合 \(n\) 大到不能直接筛。

\[ \displaystyle #\{n\le N:\mu(n)^2=1\} \]

Complexity: Time: O(n^(2/5) log log n). Space: O(n^(2/5)).

AC 记录：counting_squarefrees。

跳到代码 · GitHub ↗

Implementation

当前头文件，省略 include guard；依赖见 #include。

/// @complexity Time: O(n^(2/5) log log n).
/// Space: O(n^(2/5)).

#include <cassert>
#include <cmath>
#include <cstdint>
#include <algorithm>
#include <limits>
#include <vector>

namespace noya {

namespace squarefree_count_detail {

inline std::uint64_t floor_square_root(std::uint64_t val) {
  std::uint64_t rt = std::sqrt(static_cast<long double>(val));
  while ((rt + 1) <= val / (rt + 1)) {
    rt++;
  }
  while (rt > 0 && rt > val / rt) {
    rt--;
  }
  return rt;
}

inline std::uint64_t floor_fifth_root(std::uint64_t val) {
  auto am = [&](std::uint64_t rt) {
    unsigned __int128 prd = 1;
    for (int exp = 0; exp < 5; exp++) {
      prd *= rt;
      if (prd > val) {
        return false;
      }
    }
    return true;
  };
  std::uint64_t low = 0;
  std::uint64_t hig = 65536;
  while (hig - low > 1) {
    std::uint64_t mid = (low + hig) / 2;
    (am(mid) ? low : hig) = mid;
  }
  return low;
}

inline std::vector<std::int8_t> mobius_table(int lmt) {
  std::vector<int> mnp(lmt + 1);
  std::vector<int> ps;
  std::vector<std::int8_t> mu(lmt + 1);
  if (lmt >= 1) {
    mu[1] = 1;
  }
  for (int val = 2; val <= lmt; val++) {
    if (mnp[val] == 0) {
      mnp[val] = val;
      ps.push_back(val);
      mu[val] = -1;
    }
    for (int p : ps) {
      if (p > mnp[val] || val > lmt / p) {
        break;
      }
      mnp[val * p] = p;
      mu[val * p] = val % p == 0 ? 0 : std::int8_t(-mu[val]);
    }
  }
  return mu;
}

} // namespace squarefree_count_detail

/// @brief Count square-free positive integers at most n.  Mobius inversion
/// handles divisors up to `D = sqrt(n / n^(1/5))` directly.  The remaining
/// terms share only `O(n^(1/5))` distinct square-root quotients; their Mertens
/// values are recovered in descending order with the quotient-block
/// recurrence, reducing both time and memory from O(sqrt(n)) to O(n^(2/5)).
inline std::uint64_t squarefree_count(std::uint64_t n) {
  if (n == 0) {
    return 0;
  }
  std::uint64_t spl =
      std::max<std::uint64_t>(1, squarefree_count_detail::floor_fifth_root(n));
  std::uint64_t lim = squarefree_count_detail::floor_square_root(n / spl);
  assert(lim <= std::uint64_t(std::numeric_limits<int>::max()));
  auto mu = squarefree_count_detail::mobius_table(int(lim));
  std::vector<std::int64_t> mer(lim + 1);
  for (std::uint64_t val = 1; val <= lim; val++) {
    mer[val] = mer[val - 1] + mu[val];
  }

  __int128 sm = 0;
  for (std::uint64_t val = 1; val <= lim; val++) {
    sm += __int128(mu[val]) * (n / val / val);
  }

  std::vector<std::int64_t> lmu;
  lmu.reserve(spl > 0 ? spl - 1 : 0);
  __int128 sum = 0;
  for (std::uint64_t idx = spl; idx-- > 1;) {
    std::uint64_t arg = squarefree_count_detail::floor_square_root(n / idx);
    std::uint64_t sr = squarefree_count_detail::floor_square_root(arg);
    std::int64_t val = 1;
    std::uint64_t qb = arg / (sr + 1);
    for (std::uint64_t div = 1; div <= qb; div++) {
      val -= std::int64_t(arg / div - arg / (div + 1)) * mer[div];
    }
    for (std::uint64_t div = 2; div <= sr; div++) {
      std::uint64_t quo = arg / div;
      if (quo <= lim) {
        val -= mer[quo];
      } else {
        std::uint64_t ri = idx * div * div;
        assert(ri < spl);
        val -= lmu[spl - ri - 1];
      }
    }
    lmu.push_back(val);
    sum += val;
  }
  __int128 res = sm + sum - __int128(spl - 1) * mer[lim];
  assert(res >= 0);
  return std::uint64_t(res);
}

} // namespace noya

#ifndef NOYA_SQUAREFREE_COUNT_HPP
#define NOYA_SQUAREFREE_COUNT_HPP 1

/// @complexity Time: O(n^(2/5) log log n).
/// Space: O(n^(2/5)).

#include <cassert>
#include <cmath>
#include <cstdint>
#include <algorithm>
#include <limits>
#include <vector>

namespace noya {

namespace squarefree_count_detail {

inline std::uint64_t floor_square_root(std::uint64_t val) {
  std::uint64_t rt = std::sqrt(static_cast<long double>(val));
  while ((rt + 1) <= val / (rt + 1)) {
    rt++;
  }
  while (rt > 0 && rt > val / rt) {
    rt--;
  }
  return rt;
}

inline std::uint64_t floor_fifth_root(std::uint64_t val) {
  auto am = [&](std::uint64_t rt) {
    unsigned __int128 prd = 1;
    for (int exp = 0; exp < 5; exp++) {
      prd *= rt;
      if (prd > val) {
        return false;
      }
    }
    return true;
  };
  std::uint64_t low = 0;
  std::uint64_t hig = 65536;
  while (hig - low > 1) {
    std::uint64_t mid = (low + hig) / 2;
    (am(mid) ? low : hig) = mid;
  }
  return low;
}

inline std::vector<std::int8_t> mobius_table(int lmt) {
  std::vector<int> mnp(lmt + 1);
  std::vector<int> ps;
  std::vector<std::int8_t> mu(lmt + 1);
  if (lmt >= 1) {
    mu[1] = 1;
  }
  for (int val = 2; val <= lmt; val++) {
    if (mnp[val] == 0) {
      mnp[val] = val;
      ps.push_back(val);
      mu[val] = -1;
    }
    for (int p : ps) {
      if (p > mnp[val] || val > lmt / p) {
        break;
      }
      mnp[val * p] = p;
      mu[val * p] = val % p == 0 ? 0 : std::int8_t(-mu[val]);
    }
  }
  return mu;
}

} // namespace squarefree_count_detail

/// @brief Count square-free positive integers at most n.  Mobius inversion
/// handles divisors up to `D = sqrt(n / n^(1/5))` directly.  The remaining
/// terms share only `O(n^(1/5))` distinct square-root quotients; their Mertens
/// values are recovered in descending order with the quotient-block
/// recurrence, reducing both time and memory from O(sqrt(n)) to O(n^(2/5)).
inline std::uint64_t squarefree_count(std::uint64_t n) {
  if (n == 0) {
    return 0;
  }
  std::uint64_t spl =
      std::max<std::uint64_t>(1, squarefree_count_detail::floor_fifth_root(n));
  std::uint64_t lim = squarefree_count_detail::floor_square_root(n / spl);
  assert(lim <= std::uint64_t(std::numeric_limits<int>::max()));
  auto mu = squarefree_count_detail::mobius_table(int(lim));
  std::vector<std::int64_t> mer(lim + 1);
  for (std::uint64_t val = 1; val <= lim; val++) {
    mer[val] = mer[val - 1] + mu[val];
  }

  __int128 sm = 0;
  for (std::uint64_t val = 1; val <= lim; val++) {
    sm += __int128(mu[val]) * (n / val / val);
  }

  std::vector<std::int64_t> lmu;
  lmu.reserve(spl > 0 ? spl - 1 : 0);
  __int128 sum = 0;
  for (std::uint64_t idx = spl; idx-- > 1;) {
    std::uint64_t arg = squarefree_count_detail::floor_square_root(n / idx);
    std::uint64_t sr = squarefree_count_detail::floor_square_root(arg);
    std::int64_t val = 1;
    std::uint64_t qb = arg / (sr + 1);
    for (std::uint64_t div = 1; div <= qb; div++) {
      val -= std::int64_t(arg / div - arg / (div + 1)) * mer[div];
    }
    for (std::uint64_t div = 2; div <= sr; div++) {
      std::uint64_t quo = arg / div;
      if (quo <= lim) {
        val -= mer[quo];
      } else {
        std::uint64_t ri = idx * div * div;
        assert(ri < spl);
        val -= lmu[spl - ri - 1];
      }
    }
    lmu.push_back(val);
    sum += val;
  }
  __int128 res = sm + sum - __int128(spl - 1) * mer[lim];
  assert(res >= 0);
  return std::uint64_t(res);
}

} // namespace noya

#endif // NOYA_SQUAREFREE_COUNT_HPP

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <limits>
#include <vector>

/// @complexity Time: O(n^(2/5) log log n).
/// Space: O(n^(2/5)).

namespace noya {

namespace squarefree_count_detail {

inline std::uint64_t floor_square_root(std::uint64_t val) {
  std::uint64_t rt = std::sqrt(static_cast<long double>(val));
  while ((rt + 1) <= val / (rt + 1)) {
    rt++;
  }
  while (rt > 0 && rt > val / rt) {
    rt--;
  }
  return rt;
}

inline std::uint64_t floor_fifth_root(std::uint64_t val) {
  auto am = [&](std::uint64_t rt) {
    unsigned __int128 prd = 1;
    for (int exp = 0; exp < 5; exp++) {
      prd *= rt;
      if (prd > val) {
        return false;
      }
    }
    return true;
  };
  std::uint64_t low = 0;
  std::uint64_t hig = 65536;
  while (hig - low > 1) {
    std::uint64_t mid = (low + hig) / 2;
    (am(mid) ? low : hig) = mid;
  }
  return low;
}

inline std::vector<std::int8_t> mobius_table(int lmt) {
  std::vector<int> mnp(lmt + 1);
  std::vector<int> ps;
  std::vector<std::int8_t> mu(lmt + 1);
  if (lmt >= 1) {
    mu[1] = 1;
  }
  for (int val = 2; val <= lmt; val++) {
    if (mnp[val] == 0) {
      mnp[val] = val;
      ps.push_back(val);
      mu[val] = -1;
    }
    for (int p : ps) {
      if (p > mnp[val] || val > lmt / p) {
        break;
      }
      mnp[val * p] = p;
      mu[val * p] = val % p == 0 ? 0 : std::int8_t(-mu[val]);
    }
  }
  return mu;
}

} // namespace squarefree_count_detail

/// @brief Count square-free positive integers at most n.  Mobius inversion
/// handles divisors up to `D = sqrt(n / n^(1/5))` directly.  The remaining
/// terms share only `O(n^(1/5))` distinct square-root quotients; their Mertens
/// values are recovered in descending order with the quotient-block
/// recurrence, reducing both time and memory from O(sqrt(n)) to O(n^(2/5)).
inline std::uint64_t squarefree_count(std::uint64_t n) {
  if (n == 0) {
    return 0;
  }
  std::uint64_t spl =
      std::max<std::uint64_t>(1, squarefree_count_detail::floor_fifth_root(n));
  std::uint64_t lim = squarefree_count_detail::floor_square_root(n / spl);
  assert(lim <= std::uint64_t(std::numeric_limits<int>::max()));
  auto mu = squarefree_count_detail::mobius_table(int(lim));
  std::vector<std::int64_t> mer(lim + 1);
  for (std::uint64_t val = 1; val <= lim; val++) {
    mer[val] = mer[val - 1] + mu[val];
  }

  __int128 sm = 0;
  for (std::uint64_t val = 1; val <= lim; val++) {
    sm += __int128(mu[val]) * (n / val / val);
  }

  std::vector<std::int64_t> lmu;
  lmu.reserve(spl > 0 ? spl - 1 : 0);
  __int128 sum = 0;
  for (std::uint64_t idx = spl; idx-- > 1;) {
    std::uint64_t arg = squarefree_count_detail::floor_square_root(n / idx);
    std::uint64_t sr = squarefree_count_detail::floor_square_root(arg);
    std::int64_t val = 1;
    std::uint64_t qb = arg / (sr + 1);
    for (std::uint64_t div = 1; div <= qb; div++) {
      val -= std::int64_t(arg / div - arg / (div + 1)) * mer[div];
    }
    for (std::uint64_t div = 2; div <= sr; div++) {
      std::uint64_t quo = arg / div;
      if (quo <= lim) {
        val -= mer[quo];
      } else {
        std::uint64_t ri = idx * div * div;
        assert(ri < spl);
        val -= lmu[spl - ri - 1];
      }
    }
    lmu.push_back(val);
    sum += val;
  }
  __int128 res = sm + sum - __int128(spl - 1) * mer[lim];
  assert(res >= 0);
  return std::uint64_t(res);
}

} // namespace noya