線形篩 (osa-k素因数分解)
(Src/Number/LinearSieve.hpp)

• View this file on GitHub
• Last update: 2025-05-29 15:37:34+09:00
• Include: #include "Src/Number/LinearSieve.hpp"

概要

ある正整数 $n$ 以下の任意の正整数 $x$ について、適切な前計算の後

• 素数判定
• 素因数分解
• 約数列挙

を提供することができます。いずれも試し割りのアルゴリズムより高速です。

ライブラリの使い方

コンストラクタ

LinearSieve(u32 n)

$n$ を設定します。

制約: $n$ は正整数である

計算量: $O(n)$

size

usize size() const noexcept

コンストラクタで設定した $n$ を返します。

isPrime

bool isPrime(u32 x) const noexcept

$x$ が素数か判定します。

制約: $x$ は $n$ 以下の正整数である

計算量: $O(1)$

operator[]

bool operator[](u32 x) const noexcept

$x$ が素数か判定します。

制約: $x$ は $n$ 以下の正整数である

計算量: $O(1)$

divisor

template <std::integral T = u32>
std::vector<T> divisor(u32 x) const

$x$ の約数を列挙する。ソートされていないことに注意

制約: $x$ は $n$ 以下の正整数である

計算量: $x$ の約数の個数に比例する

factorize

std::vector factorize(u32 x) const

$x$ を素因数分解したものを返す。

PrimeFactorとは以下のような構造のclassである。

class PrimeFactor {
public:
    u32 factor(); // 素因数を返す O(1)
    u32 exponent(); // 指数を返す O(1)
};

制約: $x$ は $n$ 以下の正整数である

計算量: $O(\log x)$

参考

• 線形篩

Depends on

• Src/Number/PrimeFactor.hpp
• 標準データ型のエイリアス (Src/Template/TypeAlias.hpp)

Verified with

• Test/AtCoder/abc170_d.test.cpp
• Test/AtCoder/abc177_e.test.cpp
• ABC185-E Adjacent GCD (約数・倍数関係の高速ゼータ・メビウス変換) (Test/AtCoder/arc185_e.test.cpp)
• Test/CF/CF902-B.test.cpp

Code

#pragma once

#include "../Template/TypeAlias.hpp"
#include "./PrimeFactor.hpp"

#include <concepts>
#include <vector>
#include <utility>
#include <cassert>

namespace zawa {

class LinearSieve {
public:

    using V = u32;
    using F = PrimeFactor<V>;

private:

    std::vector<V> primes_;
    std::vector<V> lpf_;

public:

    explicit LinearSieve(V n) : primes_{}, lpf_(n + 1) {
        for (V i{2} ; i <= n ; i++) {
            if (!lpf_[i]) {
                lpf_[i] = i;
                primes_.emplace_back(i);
            }
            for (V p : primes_) {
                if (static_cast<u64>(p) * i > n) break;
                lpf_[p * i] = p;
            }
        }
    }

    V size() const noexcept {
        return static_cast<V>(lpf_.size()) - 1;
    }

    bool isPrime(V x) const noexcept {
        assert(x < lpf_.size());
        return lpf_[x] == x;
    }

    bool operator[](V x) const noexcept {
        assert(x < lpf_.size());
        return lpf_[x] == x;
    }

    std::vector<V> primes() const {
        return primes_;
    }

    // @note: response array is not sorted.
    template <std::integral T = V>
    std::vector<T> divisor(V x) const {
        assert(0u < x and x < lpf_.size());
        std::vector<T> res(1, 1u);
        while (x > 1) {
            V factor{lpf_[x]};
            u32 exponent{};
            while (lpf_[x] == factor) {
                exponent++;
                x /= lpf_[x];
            }
            usize line{res.size()};
            V now{1};
            for (u32 i{} ; i < exponent ; i++) {
                now *= factor;
                for (usize j{} ; j < line ; j++) {
                    res.emplace_back(static_cast<T>(res[j] * now));
                }
            }
        }
        return res;
    }

    std::vector<F> factorize(V x) const {
        assert(0u < x and x < lpf_.size());
        std::vector<F> res;
        while (x > 1) {
            V factor{lpf_[x]};
            u32 exponent{};
            while (lpf_[x] == factor) {
                exponent++;
                x /= lpf_[x];
            }
            res.emplace_back(factor, exponent);
        }
        return res;
    }

    i32 mobius(V x) const {
        assert(0u < x and x < lpf_.size());
        i32 res = 1;
        while (x > 1u) {
            V factor = lpf_[x];
            u32 exp = 0;
            while (lpf_[x] == factor) {
                x /= factor;
                exp++;
            }
            if (exp >= 2u) return 0;
            res *= -1;
        }
        return res;
    }
};

} // namespace zawa

#line 2 "Src/Number/LinearSieve.hpp"

#line 2 "Src/Template/TypeAlias.hpp"

#include <cstdint>
#include <cstddef>

namespace zawa {

using i16 = std::int16_t;
using i32 = std::int32_t;
using i64 = std::int64_t;
using i128 = __int128_t;

using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

using usize = std::size_t;

} // namespace zawa
#line 2 "Src/Number/PrimeFactor.hpp"

#line 4 "Src/Number/PrimeFactor.hpp"

#include <type_traits>

namespace zawa {

template <class T>
class PrimeFactor {
    static_assert(std::is_unsigned_v<T>, "T must be unsigned integer");
private:
    T factor_{};
    u32 exponent_{};

public: 
    PrimeFactor() = default;

    PrimeFactor(T factor, u32 exponent) : factor_{factor}, exponent_{exponent} {}

    inline T factor() const noexcept {
        return factor_;
    }

    inline u32 exponent() const noexcept {
        return exponent_;
    }

    friend bool operator<(const PrimeFactor<T>& lhs, const PrimeFactor<T>& rhs) {
        return lhs.factor() != rhs.factor() ?
            lhs.factor() < rhs.factor() :
            lhs.exponent() < rhs.exponent();
    }

    friend bool operator>(const PrimeFactor<T>& lhs, const PrimeFactor<T>& rhs) {
        return lhs.factor() != rhs.factor() ?
            lhs.factor() > rhs.factor() :
            lhs.exponent() > rhs.exponent();
    }
};

} // namespace zawa
#line 5 "Src/Number/LinearSieve.hpp"

#include <concepts>
#include <vector>
#include <utility>
#include <cassert>

namespace zawa {

class LinearSieve {
public:

    using V = u32;
    using F = PrimeFactor<V>;

private:

    std::vector<V> primes_;
    std::vector<V> lpf_;

public:

    explicit LinearSieve(V n) : primes_{}, lpf_(n + 1) {
        for (V i{2} ; i <= n ; i++) {
            if (!lpf_[i]) {
                lpf_[i] = i;
                primes_.emplace_back(i);
            }
            for (V p : primes_) {
                if (static_cast<u64>(p) * i > n) break;
                lpf_[p * i] = p;
            }
        }
    }

    V size() const noexcept {
        return static_cast<V>(lpf_.size()) - 1;
    }

    bool isPrime(V x) const noexcept {
        assert(x < lpf_.size());
        return lpf_[x] == x;
    }

    bool operator[](V x) const noexcept {
        assert(x < lpf_.size());
        return lpf_[x] == x;
    }

    std::vector<V> primes() const {
        return primes_;
    }

    // @note: response array is not sorted.
    template <std::integral T = V>
    std::vector<T> divisor(V x) const {
        assert(0u < x and x < lpf_.size());
        std::vector<T> res(1, 1u);
        while (x > 1) {
            V factor{lpf_[x]};
            u32 exponent{};
            while (lpf_[x] == factor) {
                exponent++;
                x /= lpf_[x];
            }
            usize line{res.size()};
            V now{1};
            for (u32 i{} ; i < exponent ; i++) {
                now *= factor;
                for (usize j{} ; j < line ; j++) {
                    res.emplace_back(static_cast<T>(res[j] * now));
                }
            }
        }
        return res;
    }

    std::vector<F> factorize(V x) const {
        assert(0u < x and x < lpf_.size());
        std::vector<F> res;
        while (x > 1) {
            V factor{lpf_[x]};
            u32 exponent{};
            while (lpf_[x] == factor) {
                exponent++;
                x /= lpf_[x];
            }
            res.emplace_back(factor, exponent);
        }
        return res;
    }

    i32 mobius(V x) const {
        assert(0u < x and x < lpf_.size());
        i32 res = 1;
        while (x > 1u) {
            V factor = lpf_[x];
            u32 exp = 0;
            while (lpf_[x] == factor) {
                x /= factor;
                exp++;
            }
            if (exp >= 2u) return 0;
            res *= -1;
        }
        return res;
    }
};

} // namespace zawa

Back to top page