Linear Sieve
(math/linear_sieve.hpp)

• View this file on GitHub
• Last update: 2023-12-28 15:52:36+09:00
• Include: #include "math/linear_sieve.hpp"
• Link: https://37zigen.com/linear-sieve/
• Link: https://atcoder.jp/contests/abc162/submissions/25095562

説明

線形篩。線形で、 $[1, N]$ の整数の最小の素因数を求める。

prime_table()

$[1, N]$ に含まれる素数を列挙する。 $O(N)$

factorize(int n)

$n$ を素因数分解する。 $O(\log N)$

divisors(int n)

$n$ の約数を列挙する。

fast_zeta(std::vector f)

$f$ を約数に関してゼータ変換（Multiple）する。 $O(N\log \log N)$

fast_mobius(std::vector F)

$F$ を約数に関してメビウス変換（Multiple）する。 $O(N\log \log N)$

pow_table(int k)

$[1, N]$ に含まれる $i$ に対して $i^k$ を計算する。 $O(N \frac{\log k}{\log N})$

inv_table()

$[1, N]$ に含まれる $i$ に対して $i^{-1} \mod m$ を計算する。 $O(N \frac{\log m}{\log N})$

Depends on

• template/int_alias.hpp

Required by

• Dirichlet Convolution (convolution/dirichlet_convolution.hpp)
• Dirichlet Series (math/DirichletSeries.hpp)

Verified with

• test/math/Sum_of_Totient_Function.test.cpp
• test/yuki/yuki_2747.test.cpp

Code

#pragma once

#include "../template/int_alias.hpp"


/*
    reference: https://37zigen.com/linear-sieve/
    verify:    https://atcoder.jp/contests/abc162/submissions/25095562
*/

#include <cassert>

#include <vector>


namespace ebi {

struct linear_sieve {
  private:
    using u64 = std::uint64_t;
    int n;
    std::vector<int> sieve;
    std::vector<int> prime;

  public:
    linear_sieve(int _n) : n(_n), sieve(std::vector<int>(_n + 1, -1)) {
        for (int i = 2; i <= n; i++) {
            if (sieve[i] < 0) {
                sieve[i] = i;
                prime.emplace_back(i);
            }
            for (auto p : prime) {
                if (u64(p) * u64(i) > u64(n) || p > sieve[i]) break;
                sieve[p * i] = p;
            }
        }
    }

    std::vector<int> prime_table() const {
        return prime;
    }

    std::vector<std::pair<int, int>> prime_power_table(int m) const {
        assert(m <= n);
        std::vector<std::pair<int, int>> table(m + 1, {1, 1});
        for (int i = 2; i <= m; i++) {
            int p = sieve[i];
            table[i] = {p, p};
            if (sieve[i / p] == p) {
                table[i] = table[i / p];
                table[i].second *= p;
            }
        }
        return table;
    }

    std::vector<std::pair<int, int>> factorize(int x) {
        assert(x <= n);
        std::vector<std::pair<int, int>> res;
        while (x > 1) {
            int p = sieve[x];
            int exp = 0;
            if (p < 0) {
                res.emplace_back(x, 1);
                break;
            }
            while (sieve[x] == p) {
                x /= p;
                exp++;
            }
            res.emplace_back(p, exp);
        }
        return res;
    }

    std::vector<int> divisors(int x) {
        assert(x <= n);
        std::vector<int> res;
        res.emplace_back(1);
        auto pf = factorize(x);
        for (auto p : pf) {
            int sz = (int)res.size();
            for (int i = 0; i < sz; i++) {
                int ret = 1;
                for (int j = 0; j < p.second; j++) {
                    ret *= p.first;
                    res.emplace_back(res[i] * ret);
                }
            }
        }
        return res;
    }

    template <class T> std::vector<T> fast_zeta(const std::vector<T> &f) {
        std::vector<T> F = f;
        int sz = f.size();
        assert(sz <= n + 1);
        for (int i = 2; i < sz; i++) {
            if (sieve[i] != i) continue;
            for (int j = (sz - 1) / i; j >= 1; j--) {
                F[j] += F[j * i];
            }
        }
        return F;
    }

    template <class T> std::vector<T> fast_mobius(const std::vector<T> &F) {
        std::vector<T> f = F;
        int sz = F.size();
        assert(sz <= n + 1);
        for (int i = 2; i < sz; i++) {
            if (sieve[i] != i) continue;
            for (int j = 1; j * i < sz; j++) {
                f[j] -= f[j * i];
            }
        }
        return f;
    }

    template <class modint> std::vector<modint> pow_table(int k) {
        std::vector<modint> table(n + 1, 1);
        table[0] = 0;
        for (int i = 2; i <= n; i++) {
            if (sieve[i] == i) {
                table[i] = modint(i).pow(k);
                continue;
            }
            table[i] = table[sieve[i]] * table[i / sieve[i]];
        }
        return table;
    }

    template <class modint> std::vector<modint> inv_table() {
        return pow_table(modint::mod() - 2);
    }
};

}  // namespace ebi

#line 2 "math/linear_sieve.hpp"

#line 2 "template/int_alias.hpp"

#include <cstdint>

namespace ebi {

using ld = long double;
using std::size_t;
using i8 = std::int8_t;
using u8 = std::uint8_t;
using i16 = std::int16_t;
using u16 = std::uint16_t;
using i32 = std::int32_t;
using u32 = std::uint32_t;
using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

}  // namespace ebi
#line 4 "math/linear_sieve.hpp"

/*
    reference: https://37zigen.com/linear-sieve/
    verify:    https://atcoder.jp/contests/abc162/submissions/25095562
*/

#include <cassert>

#include <vector>


namespace ebi {

struct linear_sieve {
  private:
    using u64 = std::uint64_t;
    int n;
    std::vector<int> sieve;
    std::vector<int> prime;

  public:
    linear_sieve(int _n) : n(_n), sieve(std::vector<int>(_n + 1, -1)) {
        for (int i = 2; i <= n; i++) {
            if (sieve[i] < 0) {
                sieve[i] = i;
                prime.emplace_back(i);
            }
            for (auto p : prime) {
                if (u64(p) * u64(i) > u64(n) || p > sieve[i]) break;
                sieve[p * i] = p;
            }
        }
    }

    std::vector<int> prime_table() const {
        return prime;
    }

    std::vector<std::pair<int, int>> prime_power_table(int m) const {
        assert(m <= n);
        std::vector<std::pair<int, int>> table(m + 1, {1, 1});
        for (int i = 2; i <= m; i++) {
            int p = sieve[i];
            table[i] = {p, p};
            if (sieve[i / p] == p) {
                table[i] = table[i / p];
                table[i].second *= p;
            }
        }
        return table;
    }

    std::vector<std::pair<int, int>> factorize(int x) {
        assert(x <= n);
        std::vector<std::pair<int, int>> res;
        while (x > 1) {
            int p = sieve[x];
            int exp = 0;
            if (p < 0) {
                res.emplace_back(x, 1);
                break;
            }
            while (sieve[x] == p) {
                x /= p;
                exp++;
            }
            res.emplace_back(p, exp);
        }
        return res;
    }

    std::vector<int> divisors(int x) {
        assert(x <= n);
        std::vector<int> res;
        res.emplace_back(1);
        auto pf = factorize(x);
        for (auto p : pf) {
            int sz = (int)res.size();
            for (int i = 0; i < sz; i++) {
                int ret = 1;
                for (int j = 0; j < p.second; j++) {
                    ret *= p.first;
                    res.emplace_back(res[i] * ret);
                }
            }
        }
        return res;
    }

    template <class T> std::vector<T> fast_zeta(const std::vector<T> &f) {
        std::vector<T> F = f;
        int sz = f.size();
        assert(sz <= n + 1);
        for (int i = 2; i < sz; i++) {
            if (sieve[i] != i) continue;
            for (int j = (sz - 1) / i; j >= 1; j--) {
                F[j] += F[j * i];
            }
        }
        return F;
    }

    template <class T> std::vector<T> fast_mobius(const std::vector<T> &F) {
        std::vector<T> f = F;
        int sz = F.size();
        assert(sz <= n + 1);
        for (int i = 2; i < sz; i++) {
            if (sieve[i] != i) continue;
            for (int j = 1; j * i < sz; j++) {
                f[j] -= f[j * i];
            }
        }
        return f;
    }

    template <class modint> std::vector<modint> pow_table(int k) {
        std::vector<modint> table(n + 1, 1);
        table[0] = 0;
        for (int i = 2; i <= n; i++) {
            if (sieve[i] == i) {
                table[i] = modint(i).pow(k);
                continue;
            }
            table[i] = table[sieve[i]] * table[i / sieve[i]];
        }
        return table;
    }

    template <class modint> std::vector<modint> inv_table() {
        return pow_table(modint::mod() - 2);
    }
};

}  // namespace ebi

Back to top page