nekolib/math/mod_recip_table_.rs
1//! 素数 $m$ を法とした逆元のテーブル。
2
3/// 素数 $m$ を法とした逆元のテーブル。
4///
5/// $\\gdef\\recip#1#2{#1\_{(#2)}\^{-1}}$
6/// 以下、$i\^{-1}\\bmod j$ を $\\recip{i}{j}$ と書く。
7///
8/// 次で定められる \\(a = (a\_0, a\_1, \\dots, a\_n)\\) を返す。
9/// \\[
10/// a\_i = \\begin{cases}
11/// \\recip{i}{m}, & \\text{if }\\recip{i}{m}\\text{ exists}; \\\\
12/// 0, & \\text{otherwise}.
13/// \\end{cases}
14/// \\]
15/// Note: $\\recip{i}{m}\\ne 0$。
16///
17/// # Idea
18///
19/// 次のことが成り立つ:
20/// - $\\recip{0}{m}$ は存在しない。
21/// - $\\recip{1}{m} = 1$ である。
22///
23/// $\\recip{(i+m)}{m} = \\recip{i}{m}$ なので、$i\\lt m$ について考える。
24///
25/// 各 $j$ ($1\\le j\\lt i$) については $\\recip{j}{m}$
26/// が得られているとする。ただし、$m$ が素数なので、これらは常に存在する。
27///
28/// $m = q\\cdot i+r$ ($0\\le r\\lt i$) とおき、$\\recip{i}{m}$ を辺々掛けると
29/// $$
30/// \\begin{aligned}
31/// m\\cdot\\recip{i}{m} &= q\\cdot i\\cdot\\recip{i}{m} + r\\cdot\\recip{i}{m} \\\\
32/// 0 &\\equiv q + r\\cdot\\recip{i}{m} \\pmod{m}.
33/// \\end{aligned}
34/// $$
35/// よって、
36/// $$ \\recip{i}{m} \\equiv -q\\cdot\\recip{r}{m} \\pmod{m}. $$
37/// $q\\gt 0$ と $\\recip{r}{m}\\gt 0$ であることから、
38/// $$ \\recip{i}{m} = m - (q\\cdot\\recip{r}{m} \\bmod m) $$
39/// となる。
40///
41/// $r\\lt i$ より $\\recip{r}{m}$
42/// は既に得られているため、$\\recip{i}{m}$ が順次 $O(1)$ time で求まる。
43///
44/// # Notes
45/// 一般のケースでは $\\gcd(r, m)=\\gcd(i, m)$ とは限らないので注意。
46/// たとえば、$\\recip{3}{8}$ を求める際、$8 = 3\\cdot 2 + 2$ なので
47/// $\\recip{2}{8}$ が必要になるが、これは存在しない。
48///
49/// 一般のケースで必要になる場合は、線形篩を用いる方法がある。See:
50/// [`LinearSieve::recips`](../struct.LinearSieve.html#method.recips)
51///
52/// # Complexity
53/// $O(n)$ time。
54///
55/// # Examples
56/// ```
57/// use nekolib::math::mod_recip_table_prime;
58///
59/// assert_eq!(mod_recip_table_prime(2, 3), [0, 1, 2]);
60/// assert_eq!(mod_recip_table_prime(4, 5), [0, 1, 3, 2, 4]);
61/// assert_eq!(mod_recip_table_prime(9, 5), [0, 1, 3, 2, 4].repeat(2));
62/// ```
63pub fn mod_recip_table_prime(n: u64, m: u64) -> Vec<u64> {
64 let mut dp = vec![0; n as usize + 1];
65 if 1 <= n {
66 dp[1] = 1;
67 }
68 for i in 2..=n.min(m - 1) {
69 let (q, r) = (m / i, m % i);
70 if dp[r as usize] > 0 {
71 dp[i as usize] = m - q * dp[r as usize] % m;
72 }
73 }
74 for i in m..=n {
75 dp[i as usize] = dp[(i - m) as usize];
76 }
77 dp
78}
79
80#[test]
81fn test() {
82 use factors_dup::FactorsDup;
83 use gcd::Gcd;
84 for m in (2_u64..=1000).filter(|m| m.factors_dup().count() == 1) {
85 let n = m - 1;
86 let actual = mod_recip_table_prime(n, m);
87 for i in 0..=n {
88 let recip = actual[i as usize];
89 if recip == 0 {
90 assert_ne!(i.gcd(m), 1);
91 } else {
92 assert!(recip < m);
93 assert_eq!(i * recip % m, 1);
94 }
95 }
96 }
97}