这个复习没有什么顺序的... 想到什么复习什么而已qwq
Miller-Rabin 质数测试
问题描述
测试一个正整数 (p) 是否为素数，需要较快且准确的测试方法。((p le 10^{18}))
算法解决
费马小定理
首先需要了解 费马小定理 。
费马小定理：对于质数 (p) 和任意整数 (a) ，有 (a^p equiv a pmod p) 。
反之，若满足 (a^p equiv a pmod p)，(p) 也有很大概率为质数。
将两边同时约去一个 (a) ，则有 (a^{p-1} equiv 1 pmod p) 。
也即是说：假设我们要测试 (n) 是否为质数。我们可以随机选取一个数 (a) ，然后计算 (a^{n-1} pmod n) ，如果结果不为 (1) ，我们可以 (100\%) 断定 (n) 不是质数。否则我们再随机选取一个新的数 (a) 进行测试。如此反复多次，如果每次结果都是 (1) ，我们就假定 (n) 是质数。
二次探测定理
该测试被称为 (Fermat) 测试。(Miller) 和 (Rabin) 在 (Fermat) 测试上，建立了 (Miller-Rabin) 质数测试算法。与 (Fermat) 测试相比，增加了一个 二次探测定理 。
二次探测定理：如果 (p) 是奇素数，则 (x^2 equiv 1 pmod p) 的解为 (x equiv 1) 或 (x equiv p - 1 pmod p) 。
这是很容易证明的：
[egin{align} x^2 equiv 1 pmod p \ x^2 -1 equiv 0 pmod p \ (x-1) (x+1) equiv 0 pmod p end{align} ]又 (ecause) (p) 为奇素数，有且仅有 (1,p) 两个因子。
( herefore) 只有两解 (x equiv 1) 或 (x equiv p - 1 pmod p) 。
如果 (a^{n-1} equiv 1 pmod n) 成立，(Miller-Rabin) 算法不是立即找另一个 (a) 进行测试，而是看 (n-1) 是不是偶数。如果 (n-1) 是偶数，另 (displaystyle u=frac{n-1}{2}) ，并检查是否满足二次探测定理即 (a^u equiv 1) 或 (a^u equiv n - 1 pmod n) 。
假设 (n=341) ，我们选取的 (a=2) 。则第一次测试时，(2^{340} mod 341 = 1) 。由于 (340) 是偶数，因此我们检查 (2^{170}) ，得到 (2^{170} mod 341=1) ，满足二次探测定理。同时由于 (170) 还是偶数，因此我们进一步检查(2^{85} mod 341=32) 。此时不满足二次探测定理，因此可以判定 (341) 不为质数。
将这两条定理合起来，也就是最常见的 (Miller-Rabin) 测试 。
单次测试失败的几率为 (displaystyle frac{1}{4}) 。那么假设我们做了 (times) 次测试，那么我们总错误的几率为 (displaystyle (frac{1}{4})^{times}) 如果 (times) 为 (20) 次，那么错误概率约为 (0.00000000000090949470) 也就是意味着不可能发生的事件。
单次时间复杂度是 (O(log n)) ，总复杂度就是 (O(times log n)) 了。注意 long long 会爆，用 __int128 就行了。
代码实现
inline bool Miller_Rabin(ll p) {
if (p <= 2) return p == 2;
if (!(p & 1)) return false;
ll u = p - 1; int power = 0;
for (; !(u & 1); u >>= 1) ++ power;
For (i, 1, times) {
ll a = randint(2, p - 1), x = fpm(a, u, p), y;
for (int i = 1; i <= power; ++ i, x = y) {
if ((y = mul(x, x, p)) == 1 && x != 1 && x != p - 1) return false;
}
if (x != 1) return false;
}
return true;
}
Pollard-Rho 大合数分解
问题描述
分解一个合数 (n) ，时间效率要求较高。（比试除法 (O(sqrt n)) 要快许多）
算法解决
生日悖论
在一个班级里，假设有 (60) 人，所有人生日不同概率是多少？
依次按人考虑，第一个人有 (1) 的概率，第二个人要与第一个人不同就是 (displaystyle 1 - frac{1}{365}) ，第三个人与前两人不同，那么就是 (displaystyle 1 - frac{2}{365}) 。那么第 (i) 个人就是 (displaystyle 1 - frac{i-1}{365}) 。
那么概率就是
假设我们带入 (n) 等于 (60) ，那么这个概率就是 (0.0058) ，也就是说基本上不可能发生。
好像是因为和大众的常识有些违背，所以称作生日悖论。
假设一年有 (N) 天，那么只要 (n ge sqrt N) 存在相同生日的概率就已经大于 (50 \%) 了。
利用伪随机数判断
本段参考了 Doggu 大佬的博客 。
我们考虑利用之前那个性质来判断，生成了许多随机数，然后两两枚举判断。
假设枚举的两个数为 (i, j) ，假设 (i > j) ，那么判断一下 (gcd(i - j, n)) 是多少，如果非 (1) ，那么我们就已经找出了一个 (n) 的因子了，这样的概率随着枚举的组数变多会越来越快。
如果我们一开始把这些用来判断的数都造出来这样，这样就很坑了。不仅内存有点恶心，而且其实效率也不够高。
考虑优化，我们用一个伪随机数去判断，不断生成两个随机数，然后像流水线一样逐个判断就行了。
有一个随机函数似乎很优秀，大部分人都用的这个：
(a) 在单次判断中是一个常数。
然后一般这种简单的随机数都会出现循环节，我们判断一下循环节就行了。
但为了节省内存，需要接下来这种算法。
Floyd 判圈法
两个人同时在一个环形跑道上起跑，一个人是另外一个人的两倍速度，那么这个人绝对会追上另外一个人。追上时，意味着其中一个人已经跑了一圈了。
然后就可以用这个算法把退出条件变松，只要有环，我们就直接退出就行了。
如果这次失败了，那么继续找另外一个 (a) 就行了。
用 Miller-Rabin 优化
这个算法对于因子多，因子值小的数 (n) 是较为优异的。
也就是对于它的一个小因子 (p) ， (Pollard-Rho) 期望在 (O(sqrt p)) 时间内找出来。
但对于 (n) 是一个质数的情况，就没有什么较好的方法快速判断了，那么复杂度就是 (O(sqrt n)) 的。
此时就退化成暴力试除法了。
此时我们考虑用前面学习的 (Miller-Rabin) 算法来优化这个过程就行了。
此时把这两个算法结合在一起，就能解决这道题啦qwq
代码实现
我似乎要特判因子为 (2) 的数，那个似乎一直在里面死循环qwq
namespace Pollard_Rho {
ll c, Mod;
ll f(ll x) { return (mul(x, x, Mod) + c) % Mod; }
ll find(ll x) {
if (!(x & 1)) return 2; Mod = x;
ll a = randint(2, x - 1), b = a; c = randint(2, x - 1);
do {
a = f(a); b = f(f(b));
ll p = __gcd(abs(a - b), x);
if (p > 1) return p;
} while (b != a);
return find(x);
}
void ReSolve(ll x, vector<ll> &factor) {
if (x <= 1) return;
if (Miller_Rabin(x)) { factor.pb(x); return; }
ll fac = find(x); ReSolve(fac, factor); ReSolve(x / fac, factor);
}
};
约瑟夫问题
问题描述
首先 (n) 个候选人围成一个圈，依次编号为 (0dots n-1) 。然后随机抽选一个数 (k) ，并 (0) 号候选人开始按从 (1) 到 (k) 的顺序依次报数，(n-1)号候选人报数之后，又再次从 (0) 开始。当有人报到 (k) 时，这个人被淘汰，从圈里出去。下一个人从 (1) 开始重新报数。
现在询问：最后一个被淘汰在哪个位置。 ((n le 10^9, k le 1000)) 。
算法解决
暴力模拟
最直观的解法是用循环链表模拟报数、淘汰的过程，复杂度是 (O(nk)) 。
递推一
令 (f[n]) 表示当有 (n) 个候选人时，最后当选者的编号。 （注意是有 (n) 个人，而不是 (n) 轮）
我们考虑用数学归纳法证明：
显然当只有 (1) 个候选人时，该候选人就是当选者，并且他的编号为 (0) 。
假设我们已经求解出了 (f[n - 1]) ，并且保证 (f[n - 1]) 的值是正确的。
现在先将 (n) 个人按照编号进行排序：
0 1 2 3 ... n-1那么第一次被淘汰的人编号一定是 (k-1) (假设 (k < n) ，若 (k > n) 则为 ((k-1) mod n) )。将被选中的人标记为 (underline{#}) ：
0 1 2 3 ... k-2 # k k+1 k+2 ... n-1第二轮报数时，起点为 (k) 这个候选人。并且只剩下 (n-1) 个选手。假如此时把 (k) 看作 (0') ，(k+1) 看作 (1') ... 则对应有：
0 1 2 3 ... k-2 # k k+1 k+2 ... n-1 n-k' ... n-2' 0' 1' 2' ... n-k-1'此时在 (0',1',...,n-2') 上再进行一次 (k) 报数的选择。而 (f[n-1]) 的值已经求得，因此我们可以直接求得当选者的编号 (s') 。
但是，该编号 (s') 是在 (n-1) 个候选人报数时的编号，并不等于 (n) 个人时的编号，所以我们还需要将(s') 转换为对应的 (s) 。通过观察，(s) 和 (s') 编号相对偏移了 (k) ，又因为是在环中，因此得到 (s = (s'+k) mod n) ，即 (f[n] = (f[n - 1] + k) pmod n) 。
然后就证明完毕了。
这个时间复杂度就是 (O(n)) 的了，算比较优秀了，但对于这题来说还是不行。
递推二
因此当 (n) 小于 (k) 时，就只能采用第一种递推的算法来计算了。
那么我们就可以用第二种递推，解决的思路仍然和上面相同，而区别在于我们每次减少的 (n) 的规模不再是 (1) 。
同样用一个例子来说明，初始 (n=10,k=4) :
初始序列：
0 1 2 3 4 5 6 7 8 9当 (7) 号进行过报数之后：
0 1 2 - 4 5 6 - 8 9而对于任意一个 (n,k) 来说，退出的候选人数量为 (displaystyle lfloor frac{n}{k} floor) 。
由于此时起点为 (8) ，则等价于：
2 3 4 - 5 6 7 - 0 1因此我们仍然可以从 (f[8]) 的结果来推导出 (f[10]) 的结果。
我们需要根据 (f[8]) 的值进行分类讨论。假设 (f[8]=s) ，则根据 (s) 和 (n mod k) 的大小关系有两种情况：
[egin{cases} s' = s - n mod k + n & (s< n mod k) \ displaystyle s' = s - n mod k + lfloor frac{(s - n mod k)}{k-1} floor & (s ge n mod k) end{cases} ]上面那种就是最后剩的的几个数，对于上面例子就是 (s={0, 1}) ，(s' = {8, 9}) 。
下面那种就是其余的几个数，对于上面例子就是 (s = {5,6,7}) ，(s' = {4,5,6}) 。
这两种就是先定位好初始位置，然后再加上前面的空隙，得到新的位置。
由于我们不断的在减小 (n) 的规模，最后一定会将 (n) 减少到小于 (k) ，此时 (displaystyle lfloor frac{n}{k} floor = 0) 。
因此当 (n) 小于 (k) 时，就只能采用第一种递推的算法来计算了。
每次我们将 (displaystyle n o n - lfloor frac{n}{k} floor) ，我们可以近似看做把 (n) 乘上 (displaystyle 1-frac{1}{k}) 。
按这样分析的话，复杂度就是 (displaystyle O(-log_{frac{k-1}{k}}n+k)) 的。这个对于 (k) 很小的时候会特别快。
代码实现
int Josephus(int n, int k) {
if (n == 1) return 0;
int res = 0;
if (n < k) {
For (i, 2, n)
res = (res + k) % i;
return res;
}
res = Josephus(n - n / k, k);
if (res < n % k)
res = res - n % k + n;
else
res = res - n % k + (res - n % k) / (k - 1);
return res;
}
扩展欧几里得
问题描述
对于三个自然数 (a,b,c) ，求解 (ax+by = c) 的 ((x, y)) 的整数解。
算法解决
首先我们要判断是否存在解，对于这个这个存在整数解的充分条件是 (gcd(a, b)~ |~c) 。
也就是说 (c) 为 (a,b) 最大公因数的一个倍数。
朴素欧几里得
对于求解 (gcd(a, b)) 我们需要用 朴素欧几里得定理 。
(gcd(a,b) = gcd(b, a mod b)) 。
这个是比较好证明的：
假设 (a = k*b+r) ，有 (r = a mod b) 。不妨设 (d) 为 (a) 和 (b) 的一个任意一个公约数，则有 (a equiv b equiv 0 pmod d) 。
由于同余的性质 (a-kb equiv r equiv 0 pmod d) 因此 (d) 是 (b) 和 (a mod b) 的公约数。
然后 (forall~ d~|gcd(a, b)) 都满足这个性质，所以这个定理成立啦。
所以我们就可以得到算 (gcd) 的一个简单函数啦。
int gcd(int a, int b) {
return !b ? a : gcd(b, a % b);
}
这个复杂度是 (O(log)) 的，十分迅速。
然后判定是否有解后，我们需要在这个基础上求一组解 ((x, y)) ，由于 (a,b,c) 都是 (gcd(a,b)) 的倍数。
对于 (a, b) 有负数的情况，我们需要将他们其中一个负数加上另外一个数直到非负。（由于前面朴素欧几里得定理是不会影响的）两个负数，直接将整个式子反号，然后放到 (c) 上就行了。
我们将它们都除以 (gcd(a, b)) ，不影响后面的计算。
也就是我们先求对于 (ax + by = 1) 且 (a ot b) （互质）求 ((x, y)) 的解。
接下来我们利用前面的朴素欧几里得定律推一波式子。
不难发现此时 (x) 变成了 (y) ， (y) 变成了 (displaystyle x - lfloor frac{a}{b} floor ~y) ，利用这个性质，我们可以递归的去求解 ((x,y)) 。
边界条件其实和前面朴素欧几里得是一样的 (b=0) 的时候，我们有 (a=1,ax+by=1) 那么此时 (x = 1, y = 0) 。
这样做完的话我们用 (O(log)) 的时间就会得到一组 ((x, y)) 的特殊解。
解系扩展
但常常我们要求对于 (x ~or~ y) 的最小非负整数解，这个的话我们需要将单个 ((x, y)) 扩展成一个解系。
如果学过不定方程的话，就可以轻易得到这个解系的，在此不过多赘述了。
要记住，((x, y)) 都需要乘上 (c) 。
然后我们直接令 (x) 为 ((x mod b + b) mod b) 就行了。
代码实现
超简短版本：
递归调用的话， (y=x', x = y') ，只需要修改 (y) 就行了。（是不是很好背啊）
void Exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) x = 1, y = 0;
else Exgcd(b, a % b, y, x), y -= a / b * x;
}
欧拉函数
定义
我们定义 (varphi (x)) 为 小于 (x) 的正整数中与 (x) 互质的数的个数，称作欧拉函数。数学方式表达就是
但需要注意，我们定义 (varphi(1) = 1) 。
性质
-
若 (x) 为质数，(varphi(x) = x - 1) 。
证明：这个很显然了，因为除了质数本身的数都与它互质。
-
若 (x = p^k) （ (p) 为质数， (x) 为单个质数的整数幂），则 (varphi(x) = (p - 1) imes p ^ {k - 1}) 。
证明：不难发现所有 (p) 的倍数都与 (x) 不互质，其他所有数都与它互质。
(p) 的倍数刚好有 (p^{k-1}) 个（包括了 (x) 本身）。
那么就有 (varphi(x) = p^k - p^{k-1} = (p - 1) imes p^{k - 1}) 。
-
若 (p, q) 互质，则有 (varphi(p imes q) = varphi(p) imes varphi(q)) ，也就是欧拉函数是个积性函数。
证明 ：
如果 (a) 与 (p) 互质 ((a<p)) ， (b) 与 (q) 互质 ((b<q)) ， (c) 与 (pq) 互质 ((c<pq)) ，则 (c) 与数对 ((a,b)) 是一一对应关系。由于 (a) 的值有 (varphi (p)) 种可能， (b) 的值有 (varphi(q)) 种可能，则数对 ((a,b)) 有 (φ(p)φ(q)) 种可能，而 (c) 的值有 (φ(pq)) 种可能，所以 (φ(pq)) 就等于 (φ(p)φ(q)) 。
具体来说这一条需要 中国剩余定理 以及 完全剩余系 才能证明，感性理解一下它的思路吧。
（后面再填） -
对于一个正整数 (x) 的质数幂分解 (displaystyle x = {p_1}^{k_1} imes {p_2}^{k_2} imes cdots imes {p_{n} }^{k_n} = prod _{i=1}^{n} {p_i}^{k_i}) 。
[displaystyle varphi(x) = x imes (1 - frac{1}{p_1}) imes (1 - frac{1}{p_2}) imes cdots imes (1 - frac{1}{p_n}) = x~prod_{i=1}^{n} (1-frac{1}{p_i}) ]证明：
我们考虑用前几条定理一起来证明。
[egin{align} varphi(x) &= prod_{i=1}^{n} varphi(p_i^{k_i}) \ &= prod_{i=1}^{n} (p_i-1) imes {p_i}^{k_i-1}\ &=prod_{i=1}^{n} {p_i}^{k_i} imes(1 - frac{1}{p_i})\ &=x~ prod_{i=1}^{n} (1- frac{1}{p_i}) end{align} ] -
若 (p) 为 (x) 的约数（ (p) 为质数, (x) 为任意正整数），我们有 $varphi(x imes p) = varphi(x) imes p $ 。
证明：
我们利用之前的第 (4) 个性质来证明就行了。
不难发现 (displaystyle prod _{i=1}^{n} (1 - frac{1}{p_i})) 是不会变的，前面的那个 (x) 会变成 (x imes p) 。
所以乘 (p) 就行了。
-
若 (p) 不是 (x) 的约数（ (p) 为质数, (x) 为任意正整数），我们有 (varphi(x imes p) = varphi(x) imes (p - 1)) 。
证明 ：(p, x) 互质，由于 (varphi) 积性直接得到。
求欧拉函数
求解单个欧拉函数
我们考虑质因数分解，然后直接利用之前的性质 (4) 来求解。
快速分解的话可以用前面讲的 (Pollard~Rho) 算法。
求解一系列欧拉函数
首先需要学习 线性筛 ，我们将其改一些地方。
质数 (p) 的 (varphi(p) = p-1) 。
对于枚举一个数 (i) 和另外一个质数 (p) 的积 (x) 。
我们在线性筛，把合数筛掉会分两种情况。
- (p) 不是 (i) 的一个约数，由于性质 (5) 就有 (varphi(x) = varphi(i) imes (p - 1)) 。
- (p) 是 (i) 的一个约数，此时我们会跳出循环，由于性质 (6) 有 (varphi(x) = varphi(i) imes p) 。
代码实现
bitset<N> is_prime; int prime[N], phi[N], cnt = 0;
void Init(int maxn) {
is_prime.set(); is_prime[0] = is_prime[1] = false; phi[1] = 1;
For (i, 2, maxn) {
if (is_prime[i]) phi[i] = i - 1, prime[++ cnt] = i;
For (j, 1, cnt) {
int res = i * prime[j];
if (res > maxn) break;
is_prime[res] = false;
if (i % prime[j]) phi[res] = phi[i] * (prime[j] - 1);
else { phi[res] = phi[i] * prime[j]; break; }
}
}
}
扩展欧拉定理
证明看 这里 ，有点难证，记结论算啦qwq
这个常常可以用于降幂这种操作。它的一个扩展就是最前面的那个 费马小定理 。
求解模线性方程组
问题描述
给定了 (n) 组除数 (m_i) 和余数 (r_i) ，通过这 (n) 组 ((m_i,r_i)) 求解一个 (x) ，使得 (x mod m_i = r_i) 这就是 模线性方程组 。
数学形式表达就是 ：
求解一个 (x) 满足上列所有方程。
算法解决
中国剩余定理
中国剩余定理提供了一个较为通用的解决方法。（似乎是唯一一个以中国来命名的定理）
如果 (m_1, m_2, dots , m_n) 两两互质，则对于任意的整数 (r_1, r_2 , dots , r_n) 方程组 ((S)) 有解，并且可以用如下方法来构造：
令 (displaystyle M = prod_{i=1} ^ {n} m_i) ，并设 (displaystyle M_i = frac{M}{m_i}) ，令 (t_i) 为 (M_i) 在模 (m_i) 意义下的逆元（也就是 (t_i imes M_i equiv 1 pmod {m_i}) ）。
不难发现这个 (t_i) 是一定存在的，因为由于前面要满足两两互质，那么 (M_i ot m_i) ，所以必有逆元。
那么在模 (M) 意义下，方程 ((S)) 有且仅有一个解 ： (displaystyle x = (sum_{i=1}^{n} r_i t_i M_i) mod M) 。
不难发现这个特解是一定满足每一个方程的，只有一个解的证明可以考虑特解 (x) 对于第 (i) 个方程，下个继续满足条件的 (x) 为 (x + m_i) 。那么对于整体来说，下个满足条件的数就是 (x + mathrm{lcm} (m_1, m_2, dots, m_n) = x + M) 。
严谨数学证明请见 百度百科 。
利用扩欧合并方程组
我们考虑合并方程组，比如从 (n=2) 开始递推。
也就等价于
此处 (k_1, k_2 in mathbb{N}) 。联立后就得到了如下一个方程：
我们令 (a = m_1, b = m_2, c = r_2 - r_1) 就变成了 (ax+by=c) 的形式，用之前讲过的 扩展欧几里得 ，可以求解。
首先先判断有无解。假设存在解，我们先解出一组 ((k_1,k_2)) ，然后带入解出 (x=x_0) 的一个特解。
我们将这个可以扩展成一个解系：
由于前面不定方程的结论， (k_1) 与其相邻解的间距为 (displaystyle frac{m_2}{gcd(m_1,m_2)}) ，又有 (x=m_1 imes k_1 + r_1) 。所以 (x) 与其相邻解的距离为 (displaystyle frac{m_1 m_2}{gcd(m_1,m_2)} = mathrm{lcm}(m_1,m_2)) 。
所以我们令 (M=mathrm{lcm}(m_1, m_2), R = x_0) 则又有新的模方程 (x equiv R pmod M) 。
然后我们就将两个方程合并成一个了，只要不断重复这个过程就能做完了。
这个比 中国剩余定理 优秀的地方就在于它这个不需要要求 (m) 两两互质，并且可以较容易地判断无解的情况。
代码实现
注意有些细节，比如求两个数的 (mathrm{gcd}) 的时候，一定要先除再乘，防止溢出！！
ll mod[N], rest[N];
ll Solve() {
For (i, 1, n - 1) {
ll a = mod[i], b = mod[i + 1], c = rest[i + 1] - rest[i], gcd = __gcd(a, b), k1, k2;
if (c % gcd) return - 1;
a /= gcd; b /= gcd; c /= gcd;
Exgcd(a, b, k1, k2);
k1 = ((k1 * c) % b + b) % b;
mod[i + 1] = mod[i] / __gcd(mod[i], mod[i + 1]) * mod[i + 1] ;
rest[i + 1] = (mod[i] * k1 % mod[i + 1] + rest[i]) % mod[i + 1];
}
return rest[n];
}
扩展卢卡斯
问题描述
求
(1 le m le n le 10^{18}, 2 le p le 10^6) 不保证 (p) 为质数。
问题求解
虽说是扩展卢卡斯，但是和卢卡斯定理没有半点关系。
用了几个性质，首先可以考虑 (p) 分解质因数。
假设分解后
我们可以对于每个 (p_i) 单独考虑，假设我们求出了 (displaystyle {n choose m} mod {{p_i}^{k_i}}) 的值。
然后我们可以考虑用前文提到的 (CRT) 来进行合并（需要用扩欧求逆元）。
问题就转化成如何求 (displaystyle {n choose m} mod {{p_i}^{k_i}}) 了。
首先我们考虑它有多少个关于 (p_i) 的因子（也就是多少次方）。
我们令 (f(n)) 为 (n!) 中包含 (p_i) 的因子数量，那么 (displaystyle {n choose m} = frac{n!}{m!(n - m)!}) ，所以因子数量就为 (f(n) - f(m) - f(n - m)) 。
那如何求 (f(n)) 呢。
我们举出一个很简单的例子来讨论。
(n=19,p_i=3,k_i=2) 时:
不难发现每次把 (p_i) 倍数提出来的东西，就是 (displaystyle lfloor frac{n}{p_i} floor !) ，那么很容易得到如下递推式：
不难发现这个求单个的复杂度是 (O(log n)) 的，十分优秀。
然后考虑剩下与 (p_i) 互不相关的如何求，不难发现在 ({p_i}^{k_i}) 意义下会存在循环节（比如前面的 ((1∗2∗4∗5∗7∗8)^2) ，可能最后多出来一个或者多个不存在循环节中的数。
不难发现循环节长度是 (< {p_i}^{k_i}) ，因为同余的在 (> {p_i}^{k_i}) 之后马上出现，然后把多余的 (displaystyle lfloor frac{n}{p_i} floor) 的部分递归处理就行了。
然后就可以快速处理出与 (p_i) 无关的了，最后合并一下就行了。
简介算法流程：
- 对于每个 ({p_i}^{k_i}) 单独考虑。
- 用 (f(n)) 处理出有关于 (p_i) 的因子个数。
- 然后递归处理出无关 (p_i) 的组合数的答案。
- 把这两个乘起来合并即可。
- 最后用 (CRT) 合并所有解即可。
瞎分析的复杂度（ Great_Influence 博客上抄的）。（不可靠）
代码解决
#include <bits/stdc++.h>
#define For(i, l, r) for(register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for(register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
#define DEBUG(...) fprintf(stderr, __VA_ARGS__)
using namespace std;
typedef long long ll;
template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }
void File() {
#ifdef zjp_shadow
freopen ("P4720.in", "r", stdin);
freopen ("P4720.out", "w", stdout);
#endif
}
ll n, m, p;
void Exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) x = 1, y = 0;
else Exgcd(b, a % b, y, x), y -= a / b * x;
}
inline ll fpm(ll x, ll power, ll Mod) {
ll res = 1;
for (; power; power >>= 1, (x *= x) %= Mod)
if (power & 1) (res *= x) %= Mod;
return res;
}
inline ll fac(ll n, ll pi, ll pk) {
if (!n) return 1;
ll res = 1;
For (i, 2, pk) if (i % pi) (res *= i) %= pk;
res = fpm(res, n / pk, pk);
For (i, 2, n % pk) if (i % pi) (res *= i) %= pk;
return res * fac(n / pi, pi, pk) % pk;
}
inline ll Inv(ll n, ll Mod) {
ll x, y; Exgcd(n, Mod, x, y);
return (x % Mod + Mod) % Mod;
}
inline ll CRT(ll b, ll Mod) {
return b * Inv(p / Mod, Mod) % p * (p / Mod) % p;
}
inline ll factor(ll x, ll Mod) {
return x ? factor(x / Mod, Mod) + (x / Mod) : 0;
}
inline ll Comb(ll n, ll m, ll pi, ll pk) {
ll k = factor(n, pi) - factor(m, pi) - factor(n - m, pi);
if (!fpm(pi, k, pk)) return 0;
return fac(n, pi, pk) * Inv(fac(m, pi, pk), pk) % pk * Inv(fac(n - m, pi, pk), pk) % pk * fpm(pi, k, pk) % pk;
}
inline ll ExLucas(ll n, ll m) {
ll res = 0, tmp = p;
For (i, 2, sqrt(p + .5)) if (!(tmp % i)) {
ll pk = 1; while (!(tmp % i)) pk *= i, tmp /= i;
(res += CRT(Comb(n, m, i, pk), pk)) %= p;
}
if (tmp > 1) (res += CRT(Comb(n, m, tmp, tmp), tmp)) %= p;
return res;
}
int main () {
File();
cin >> n >> m >> p;
printf ("%lld
", ExLucas(n, m));
return 0;
}