nekolib/math/mod_tetration.rs
1//! tetration。
2
3use super::euler_phi;
4
5use euler_phi::EulerPhi;
6
7/// tetration。
8///
9/// ${}^b a\\bmod n$ を求める。${}^\\bullet \\bullet$ は次のように定義される。
10/// $$ {}^b a = \\begin{cases}
11/// 1, & \\text{if } b = 0; \\\\
12/// a^{\\left({}^{(b-1)} a\\right)}, & \\text{if } b \\ge 1.
13/// \\end{cases} $$
14/// 直感的に書けば、$\\underbrace{a^{(a^{(\\cdots ^a)})}}\_{b\\text{ many}} \\bmod n$ である。
15///
16/// # Idea
17/// 大変大きくなりうる $z$ に対して $a^z\\bmod n$ を求めることを考える。
18/// このとき、[`dlog`] の [Idea] と同じ議論から、ある $\\mu$, $\\lambda$ が存在して
19/// $z\\lt\\mu$ または $z=\\mu+q\\cdot\\lambda+r$ とでき、後者のとき
20/// $a^z\\equiv a^{\\mu+r}\\pmod{n}$ が成り立つ。
21///
22/// ここで、$\\mu\\le\\log\_2(n)$, $\\lambda\\sqsubseteq\\varphi(n)$ である。
23/// さらに、任意の $n$ に対して $\\log\_2(n)\\le\\varphi(n)$ なので、$z\\ge\\varphi(n)$
24/// ならば $z\\ge\\mu$ とわかる。よって、以下のようにできる。
25/// $$ \\begin{aligned}
26/// a^z \\equiv \\begin{cases}
27/// a^z, & \\text{if }z\\lt \\varphi(n); \\\\
28/// a^{(z\\bmod\\varphi(n))+\\varphi(n)}, & \\text{otherwise}.
29/// \\end{cases}
30/// \\end{aligned} $$
31///
32/// 直感的には、指数部が $\\varphi(n)$ 以上であればすでに周期の中に入っており、入った後は
33/// $\\varphi(n)$ を法として合同かつ $\\varphi(n)$ 以上の値さえ得られれば十分ということである。
34///
35/// [`dlog`]: fn.dlog.html
36/// [Idea]: fn.dlog.html#idea
37///
38/// ## When $b$ is large
39/// 前述のように、${}^b a\\bmod{n}$ を求める際に ${}^{b-1} a$ を法 ${\\varphi(n)}$ で考える。
40/// その次は $\\varphi(\\varphi(n)), \\varphi(\\varphi(\\varphi(n))), \\dots$ と続く。
41/// $\\varphi^\\star(n)$ 段では考えるべき法が $1$ となり、それより上の段のことは無視できる。
42///
43/// そこで、$\\varphi^\\star(n)$ を考える。奇素数 $p$ に対して $\\varphi(p^e)=p^{e-1}\\cdot(p-1)$
44/// が偶数であることと、$\\varphi(2^e)=2^{e-1}$ であることから、$\\varphi(\\varphi(n))\\lt n/2$
45/// が成り立つ。すなわち、$\\varphi^\\star(n)\\le 2\\log(n)$ である[^1]。
46///
47/// [^1]: ゆるゆるな bound である。実際にはもっと速く減りそう。
48///
49/// よって、$b\\ge 2\\log(n)$ であれば ${}^{b+1} a\\equiv {}^b a\\pmod{n}$ となる。
50///
51/// ## Common bugs
52/// 繰り返し二乗法で $\\varphi(n)$ 以上か判断しつつ $a^z\\bmod\\varphi(n)$ を求める際、
53/// $a^{2^\\bullet}$ が $\\varphi(n)$ 以上になった時点で $a^z\\ge\\varphi(n)$
54/// と判断してしまうと、誤検出してしまう場合がある。
55/// ```ignore
56/// fn mod_pow(mut a: u64, mut b: u64, n: u64) -> u64 {
57/// let mut res = 1 % n;
58/// let mut large = false;
59/// while b > 0 {
60/// if b & 1 == 1 {
61/// res *= a;
62/// if res >= n { res %= n; large = true; }
63/// }
64/// a *= a;
65/// if a >= n { a %= n; large = true; } // !
66/// b >>= 1;
67/// }
68/// if large { res + n } else { res }
69/// }
70/// ```
71/// 最後のループで初めて `a >= n` になると、`res < n` なのに `res + n` が返ってしまう。
72/// このバグや、あるいはそもそも $z\\ge\\varphi(n)$ を仮定していることなどにより、
73/// ${}^3 2\\bmod 32 = 0$ としてしまうコードがたくさん提出されている。
74/// [$\\bullet$](https://judge.yosupo.jp/submission/4054)
75/// [$\\bullet$](https://judge.yosupo.jp/submission/4564)
76/// [$\\bullet$](https://judge.yosupo.jp/submission/12501)
77/// [$\\bullet$](https://judge.yosupo.jp/submission/18734)
78/// [$\\bullet$](https://judge.yosupo.jp/submission/23725)
79/// [$\\bullet$](https://judge.yosupo.jp/submission/25108)
80/// [$\\bullet$](https://judge.yosupo.jp/submission/28794)
81/// [$\\bullet$](https://judge.yosupo.jp/submission/36536)
82/// [$\\bullet$](https://judge.yosupo.jp/submission/38102)
83/// [$\\bullet$](https://judge.yosupo.jp/submission/39646)
84/// [$\\bullet$](https://judge.yosupo.jp/submission/40708)
85/// [$\\bullet$](https://judge.yosupo.jp/submission/42416)
86///
87/// # Complexity
88/// $O(\\sqrt{n})$ time.
89///
90/// 律速は、$\\varphi(n), \\varphi(\\varphi(n)), \\dots$ を求めるパートであり、
91/// $O(\\sum\_{i=0}^{\\varphi^\\star(n)} \\sqrt{n/2^i}) = O(\\sqrt{n})$ である。
92///
93/// # Examples
94/// ```
95/// use nekolib::math::ModTetration;
96///
97/// let n = 10_u64.pow(9);
98///
99/// assert_eq!(0_u64.mod_tetration(0, n), 1);
100/// assert_eq!(0_u64.mod_tetration(1, n), 0);
101/// assert_eq!(0_u64.mod_tetration(2, n), 1);
102/// assert_eq!(0_u64.mod_tetration(3, n), 0);
103///
104/// assert_eq!(1_u64.mod_tetration(0, n), 1);
105/// assert_eq!(1_u64.mod_tetration(1, n), 1);
106///
107/// assert_eq!(2_u64.mod_tetration(0, n), 1);
108/// assert_eq!(2_u64.mod_tetration(1, n), 2);
109/// assert_eq!(2_u64.mod_tetration(2, n), 4);
110/// assert_eq!(2_u64.mod_tetration(3, n), 16);
111/// assert_eq!(2_u64.mod_tetration(4, n), 65536);
112///
113/// assert_eq!(3_u64.mod_tetration(9, n), 64_195_387);
114/// assert_eq!(3_u64.mod_tetration(10, n), 464_195_387);
115/// assert_eq!(3_u64.mod_tetration(11, n), 464_195_387);
116/// assert_eq!(3_u64.mod_tetration(99, n), 464_195_387);
117/// ```
118///
119/// # Notations
120/// ${}^b a$ は $a\\uparrow\\uparrow b$ (Knuth's up-arrow notation) や
121/// $a\\to b\\to 2$ (Conway chained arrow notation) などとも表記される。
122pub trait ModTetration {
123 fn mod_tetration(self, b: Self, n: Self) -> Self;
124}
125
126trait ModTetrationInternal {
127 fn rec(self, b: Self, n: Self) -> Self;
128 fn mod_pow(self, b: Self, n: Self, large: bool) -> Self;
129}
130
131macro_rules! impl_uint {
132 ($t:ty) => {
133 impl ModTetration for $t {
134 fn mod_tetration(self, b: Self, n: Self) -> Self {
135 let a = self;
136 match (a, b, n) {
137 (.., 1) => return 0,
138 (_, 0, _) => return 1,
139 (_, 1, _) => return a % n,
140 _ => match a.rec(b, n) {
141 z if z >= n => z - n,
142 z => z,
143 }
144 }
145 }
146 }
147 impl ModTetrationInternal for $t {
148 fn rec(self, b: Self, n: Self) -> Self {
149 let a = self;
150 match (a, b, n) {
151 (0, ..) => return 1 - b % 2,
152 (1, ..) => return 1,
153 (.., 1) => return 1,
154 (_, 1, _) => return a,
155 _ => {
156 let phi = n.euler_phi();
157 let res = a.rec(b - 1, phi);
158 (a % n).mod_pow(res, n, res >= phi || a >= phi)
159 }
160 }
161 }
162 fn mod_pow(self, mut b: Self, n: Self, mut large: bool) -> Self {
163 let mut a = self;
164 let mut res = 1;
165 let mut large_buf = false;
166 while b > 0 {
167 if b & 1 == 1 {
168 res *= a;
169 large |= large_buf;
170 if res >= n {
171 res %= n;
172 large = true;
173 }
174 }
175 a *= a;
176 if a >= n {
177 a %= n;
178 large_buf = true;
179 }
180 b >>= 1;
181 }
182 if large { res + n } else { res }
183 }
184 }
185 };
186 ( $($t:ty)* ) => { $(impl_uint!($t);)* };
187}
188
189impl_uint!(u8 u16 u32 u64 u128 usize);
190
191#[test]
192fn test() {
193 for n in 1_u64..100000 {
194 if 2.mod_tetration(2, n) != 4 % n {
195 eprintln!("{:?}", (2, 2, n));
196 }
197 if 2.mod_tetration(3, n) != 16 % n {
198 eprintln!("{:?}", (2, 3, n));
199 }
200 if 2.mod_tetration(4, n) != 65536 % n {
201 eprintln!("{:?}", (2, 4, n));
202 }
203 if 2.mod_tetration(5, n) != 2_u64.mod_pow(65536, n, true) - n {
204 eprintln!("{:?}", (2, 5, n));
205 }
206
207 if 3.mod_tetration(2, n) != 27 % n {
208 eprintln!("{:?}", (3, 2, n));
209 }
210 if 3.mod_tetration(3, n) != 7_625_597_484_987 % n {
211 eprintln!("{:?}", (3, 3, n));
212 }
213 }
214}