1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109 | #pragma once
#include <cassert><--- Include file: not found. Please note: Cppcheck does not need standard library headers to get proper results.
#include <vector><--- Include file: not found. Please note: Cppcheck does not need standard library headers to get proper results.
#include <set><--- Include file: not found. Please note: Cppcheck does not need standard library headers to get proper results.
#include <map><--- Include file:
#include <utility>
#include <algorithm><--- Include file: not found. Please note: Cppcheck does not need standard library headers to get proper results.
#include <random>
#include <ranges><--- Include file: not found. Please note: Cppcheck does not need standard library headers to get proper results.
#include "snippet/internal/types.hpp"
#include "internal/dev_env.hpp"
#include "numeric/internal/primality_test.hpp"
#include "numeric/modular/modint_interface.hpp"
#include "numeric/bit.hpp"
#include "random/engine.hpp"
namespace uni {
namespace internal {
//Thanks to: https://github.com/NyaanNyaan/library/blob/master/prime/fast-factorize.hpp
namespace fast_factorize_impl {
namespace internal {
// Pollard's rho algorithm
template<modint_family Mint, class T>
T find_factor(const T n) noexcept(NO_EXCEPT) {
if(~n & 1) return 2;
if(is_prime<Mint>(n)) return n;
assert(static_cast<u64>(Mint::mod()) == n);
Mint rr;
auto f = [&](const Mint& x) noexcept(NO_EXCEPT) { return x * x + rr; };
static random_engine_64bit rand(std::random_device{}());
auto rand_ = [&]() noexcept(NO_EXCEPT) { return rand() % (n - 2) + 2; };
while(true) {
Mint x, y, ys, q = Mint::one;
rr = rand_(), y = rand_();
T g = 1;
constexpr int m = 128;
for(int r = 1; g == 1; r <<= 1) {
x = y;
for(int i = 0; i < r; ++i) y = f(y);
for(int k = 0; g == 1 && k < r; k += m) {
ys = y;
for(int i = 0; i < m && i < r - k; ++i) q *= x - (y = f(y));
g = uni::binary_gcd(q.val(), n);
}
}
if(g == n) {
do {
g = uni::binary_gcd((x - (ys = f(ys))).val(), n);
} while(g == 1);
}
if(g != n) return g;
}
assert(false);
}
template<modint_family Small, modint_family Large, class R>
void factorize(const u64 n, R *const res) noexcept(NO_EXCEPT) {
if(n <= 1) return;
u64 p;
if constexpr(std::same_as<Small, Large>) p = find_factor<Small, typename Small::value_type>(n);
else {
if(n <= Small::max()) p = find_factor<Small, typename Small::value_type>(n);
else p = find_factor<Large, typename Large::value_type>(n);
}
if(p == static_cast<u64>(n)) {
res->emplace_back(p);
return;
}
factorize<Small, Large>(p, res);
factorize<Small, Large>(n / p, res);
}
} // namespace internal
} // namespace fast_factorize_impl
using fast_factorize_impl::internal::factorize;
} // namespace internal
} // namespace uni
|