|
8 | 8 | #include <vector>
|
9 | 9 | #include <bit>
|
10 | 10 | namespace cp_algo::math {
|
| 11 | + std::vector<int64_t> factorize(int64_t m); |
| 12 | + |
| 13 | + int64_t euler_phi(int64_t m) { |
| 14 | + auto primes = factorize(m); |
| 15 | + auto [from, to] = std::ranges::unique(primes); |
| 16 | + primes.erase(from, to); |
| 17 | + int64_t ans = m; |
| 18 | + for(auto it: primes) { |
| 19 | + ans -= ans / it; |
| 20 | + } |
| 21 | + return ans; |
| 22 | + } |
| 23 | + template<modint_type base> |
| 24 | + int64_t period(base x) { |
| 25 | + auto ans = euler_phi(base::mod()); |
| 26 | + base x0 = bpow(x, ans); |
| 27 | + for(auto t: factorize(ans)) { |
| 28 | + while(ans % t == 0 && x0 * bpow(x, ans / t) == x0) { |
| 29 | + ans /= t; |
| 30 | + } |
| 31 | + } |
| 32 | + return ans; |
| 33 | + } |
| 34 | + // Find min non-negative x s.t. a*b^x = c (mod m) |
| 35 | + std::optional<uint64_t> discrete_log(int64_t b, int64_t c, uint64_t m, int64_t a = 1) { |
| 36 | + if(std::abs(a - c) % m == 0) { |
| 37 | + return 0; |
| 38 | + } |
| 39 | + if(std::gcd(a, m) != std::gcd(a * b, m)) { |
| 40 | + auto res = discrete_log(b, c, m, a * b % m); |
| 41 | + return res ? std::optional(*res + 1) : res; |
| 42 | + } |
| 43 | + // a * b^x is periodic here |
| 44 | + using base = dynamic_modint; |
| 45 | + return base::with_mod(m, [&]() -> std::optional<uint64_t> { |
| 46 | + size_t sqrtmod = std::max<size_t>(1, std::sqrt(m) / 2); |
| 47 | + std::unordered_map<int64_t, int> small; |
| 48 | + base cur = a; |
| 49 | + for(size_t i = 0; i < sqrtmod; i++) { |
| 50 | + small[cur.getr()] = i; |
| 51 | + cur *= b; |
| 52 | + } |
| 53 | + base step = bpow(base(b), sqrtmod); |
| 54 | + cur = 1; |
| 55 | + for(size_t k = 0; k < m; k += sqrtmod) { |
| 56 | + auto it = small.find((base(c) * cur).getr()); |
| 57 | + if(it != end(small)) { |
| 58 | + auto cand = base::with_mod(period(base(b)), [&](){ |
| 59 | + return base(it->second - k); |
| 60 | + }).getr(); |
| 61 | + if(base(a) * bpow(base(b), cand) == base(c)) { |
| 62 | + return cand; |
| 63 | + } else { |
| 64 | + return std::nullopt; |
| 65 | + } |
| 66 | + } |
| 67 | + cur *= step; |
| 68 | + } |
| 69 | + return std::nullopt; |
| 70 | + }); |
| 71 | + } |
11 | 72 | // https://en.wikipedia.org/wiki/Berlekamp-Rabin_algorithm
|
12 | 73 | template<modint_type base>
|
13 | 74 | std::optional<base> sqrt(base b) {
|
@@ -98,17 +159,8 @@ namespace cp_algo::math {
|
98 | 159 | int64_t primitive_root(int64_t p) {
|
99 | 160 | using base = dynamic_modint;
|
100 | 161 | return base::with_mod(p, [p](){
|
101 |
| - auto fact = factorize(p - 1); |
102 |
| - auto is_primitive_root = [&](base x) { |
103 |
| - for(auto t: fact) { |
104 |
| - if(bpow(x, (p - 1) / t) == 1) { |
105 |
| - return false; |
106 |
| - } |
107 |
| - } |
108 |
| - return true; |
109 |
| - }; |
110 | 162 | base t = 1;
|
111 |
| - while(!is_primitive_root(t)) { |
| 163 | + while(period(t) != p - 1) { |
112 | 164 | t = random::rng();
|
113 | 165 | }
|
114 | 166 | return t.getr();
|
|
0 commit comments