Struct nekolib::math::linear_sieve::LinearSieve

pub struct LinearSieve { /* private fields */ }

線形篩。

Idea

各整数の最小素因数を求めつつ、素数のリストを作成していく。

$\gdef\lpf#1{\operatorname{lpf}(#1)}$ $i$ の最小素因数を $\mathrm{lpf}(i)$ と書く。$j \le \lpf{i}$ なる各素数 $j$ に対して、$\lpf{i\times j} = j$ とわかる。 素因数分解の一意性から、各整数の最小素因数の更新は一回ずつしか行われず、線形時間で構築できる。

また、$\lpf{i}$ が $\lpf{i / {\lpf{i}}}$ と等しいかで場合分けしながら DP することで、各 $i$ が $\lpf{i}$ で何回割り切れるかも求められる。 なお、これは DP せず各 $i$ に対して愚直に計算しても $O(n)$ になることが示せるらしい。

Complexity

演算	時間計算量
new	$\Theta(n)$
is_prime	$\Theta(1)$
least_factor	$\Theta(1)$
factors_dup	$\Theta(1)$ delay
factors	$\Theta(1)$ delay
euler_phi	$\Theta(\omega(n))$
euler_phi_star	$O(\omega(n)\log(n))$
divisors	$O(\sigma(n))$
divisors_count	$O(\omega(n))$
divisors_sum	$O(\omega(n))$
primes	$\Theta(1)$ delay
recips	$O(n)$

$n$ の素因数の個数を $\omega(n)$ とすると、以下の式が成り立つらしい。 $$ \sum_{i\le n} \omega(i) = n\ln(\ln(n)) + O(n). $$

また、$n$ 以下の素数の個数を $\pi(n)$ とすると、以下の式が成り立つ（素数定理）。 $$ \pi(n) = \frac{n}{\ln(n)} + O{\left(\frac{n}{\ln(n)^2}\right)}. $$

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert!(!sieve.is_prime(1));
assert!(sieve.is_prime(2));
assert!(sieve.is_prime(23));
assert!(!sieve.is_prime(24));

assert_eq!(sieve.least_factor(1), None);
assert_eq!(sieve.least_factor(3), Some(3));
assert_eq!(sieve.least_factor(24), Some(2));

assert_eq!(sieve.factors_dup(1).next(), None);
assert_eq!(sieve.factors_dup(60).collect::<Vec<_>>(), vec![2, 2, 3, 5]);

assert_eq!(
    sieve.primes().take(10).collect::<Vec<_>>(),
    vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
);

References

• https://cp-algorithms.com/algebra/prime-sieve-linear.html
• https://twitter.com/hidesugar2/status/1431243186362458114
• https://maspypy.com/%E7%B4%A0%E6%95%B0%E3%81%AB%E9%96%A2%E3%81%99%E3%82%8B%E4%B8%8A%E3%81%8B%E3%82%89%E3%81%AE%E8%A9%95%E4%BE%A1%EF%BC%88%E5%88%9D%E7%AD%89%E7%9A%84%E3%81%AA%E8%A8%BC%E6%98%8E%EF%BC%89

Implementations

impl LinearSieve

pub fn new(n: usize) -> Self

$n$ 以下の自然数に対する篩を用意する。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);

pub fn is_prime(&self, n: usize) -> bool

$n$ が素数であれば true を返す。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert!(!sieve.is_prime(1));
assert!(sieve.is_prime(2));
assert!(sieve.is_prime(23));
assert!(!sieve.is_prime(24));

pub fn least_factor(&self, n: usize) -> Option<usize>

$n$ の最小素因数を返す。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.least_factor(1), None);
assert_eq!(sieve.least_factor(3), Some(3));
assert_eq!(sieve.least_factor(24), Some(2));

pub fn factors_dup(&self, n: usize) -> impl Iterator<Item = usize> + '_

$n$ の素因数を列挙する。重複あり。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.factors_dup(1).next(), None);
assert_eq!(sieve.factors_dup(60).collect::<Vec<_>>(), vec![2, 2, 3, 5]);

pub fn factors(&self, n: usize) -> impl Iterator<Item = (usize, u32)> + '_

$n$ を素因数分解する。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.factors(1).next(), None);
assert_eq!(
    sieve.factors(60).collect::<Vec<_>>(),
    vec![(2, 2), (3, 1), (5, 1)]
);

pub fn euler_phi(&self, n: usize) -> usize

$\phi(n)$ を求める。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.euler_phi(1), 1);
assert_eq!(sieve.euler_phi(35), 24);
assert_eq!(sieve.euler_phi(60), 16);

pub fn euler_phi_star(&self, n: usize) -> usize

$\phi^\star(n)$ を求める。

$n$ に $\phi$ を繰り返し適用して $1$ にするために必要な最小回数である。 $n\gt 1$ に対して $\phi(\phi(n)) \le n/2$ なので、 $\phi^\star(n) = O(\log(n))$ である。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.euler_phi_star(1), 0);
assert_eq!(sieve.euler_phi_star(35), 5);
assert_eq!(sieve.euler_phi_star(60), 5);

assert_eq!(sieve.euler_phi(35), 24);
assert_eq!(sieve.euler_phi(24), 8);
assert_eq!(sieve.euler_phi(8), 4);
assert_eq!(sieve.euler_phi(4), 2);
assert_eq!(sieve.euler_phi(2), 1);

assert_eq!(sieve.euler_phi(60), 16);
assert_eq!(sieve.euler_phi(16), 8);

pub fn divisors( &self, n: usize ) -> impl Iterator<Item = usize> + DoubleEndedIterator

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

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.divisors(1).next(), Some(1));
assert_eq!(sieve.divisors(1).nth(1), None);
assert_eq!(
    sieve.divisors(60).collect::<Vec<_>>(),
    vec![1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, 60]
);

pub fn divisors_count(&self, n: usize) -> usize

$n$ の約数の個数を返す。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.divisors_count(1), 1);
assert_eq!(sieve.divisors_count(60), sieve.divisors(60).count());

pub fn divisors_sum(&self, n: usize) -> usize

$n$ の約数の総和を返す。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.divisors_sum(1), 1);
assert_eq!(sieve.divisors_sum(8), 15);
assert_eq!(sieve.divisors_sum(60), 168);

pub fn primes(&self) -> impl Iterator<Item = usize> + DoubleEndedIterator + '_

素数を列挙する。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(
    sieve.primes().take(10).collect::<Vec<_>>(),
    vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
);

pub fn recips(&self, n: usize, m: usize) -> Vec<usize>

法 $m$ での逆元を返す。

$\gdef\recip#1#2{#1^{-1}_{(#2)}}$ $i^{-1}\bmod j$ を $\recip{i}{j}$ と書く。 次で定められる $a = (a_0, a_1, \dots, a_n)$ を返す。 $$ a_i = \begin{cases} \recip{i}{m}, & \text{if }\recip{i}{m}\text{ exists}; \\ 0, & \text{otherwise}. \end{cases} $$ Note: $\recip{i}{m}\ne 0$。

Idea

$$ \recip{i}{m} \equiv \recip{\lpf{i}}{m}\cdot\recip{(i/{\lpf{i}})}{m}\pmod{m} $$ に基づく。素数は $O(n/\log(n))$ 個しかないため、互除法で愚直に求めても全体では $O(n)$ 時間となる。

See also

https://37zigen.com/linear-sieve/#i-4

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(60);
assert_eq!(sieve.recips(0, 1), [0]);
assert_eq!(sieve.recips(4, 5), [0, 1, 3, 2, 4]);
assert_eq!(sieve.recips(9, 10), [0, 1, 0, 7, 0, 0, 0, 3, 0, 9]);

pub fn dp<T>( &self, zero: T, one: T, eq: impl Fn(&T, usize) -> T, gt: impl Fn(&T, usize) -> T ) -> Vec<T>

最小素因数を用いて DP を行う。

関数 $f$ であって、任意の $i\gt 1$ に対して $$ f(i) = \begin{cases} g_{=}(f(i/j), \lpf{i}), & \text{if }\lpf{i} = \lpf{i/j}; \\ g_{\gt}(f(i/j), \lpf{i}), & \text{if }\lpf{i} \gt \lpf{i/j} \\ \end{cases} $$ となるものを計算する。ただし $j = \lpf{i}$ とする。

(zero, one, eq, gt) はそれぞれ $(f(0), f(1), g_{=}, g_{\gt})$ である。

Examples

use nekolib::math::LinearSieve;

let sieve = LinearSieve::new(10);

// Moebius mu function
let mu = sieve.dp(0, 1, |&x, _| 0, |&x, _| -x);
assert_eq!(mu, [0, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1]);

// Euler phi function
let phi = sieve.dp(0, 1, |&x, p| x * p, |&x, p| x * (p - 1));
assert_eq!(phi, [0, 1, 1, 2, 2, 4, 2, 6, 4, 6, 4]);

// the number of distinct prime factors
let omega = sieve.dp(0, 0, |&x, _| x, |&x, _| x + 1);
assert_eq!(omega, [0, 0, 1, 1, 1, 1, 2, 1, 1, 1, 2]);

// the number of prime factors
let cap_omega = sieve.dp(0, 0, |&x, _| x + 1, |&x, _| x + 1);
assert_eq!(cap_omega, [0, 0, 1, 1, 2, 1, 2, 1, 3, 2, 2]);

// sum of divisors
let eq = |&(prod, sum, pow): &_, p| (prod, sum + pow * p, pow * p);
let gt = |&(prod, sum, _): &_, p| (prod * sum, 1 + p, p);
let sigma: Vec<_> = sieve.dp((1, 1, 1), (1, 1, 1), eq, gt)
    .into_iter().map(|(prod, sum, _)| prod * sum)
    .collect();
assert_eq!(sigma, [1, 1, 3, 4, 7, 6, 12, 8, 15, 13, 18]);