基础数论复习
这个复习没有什么顺序的... 想到什么复习什么而已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\) 。
这是很容易证明的：
\[\begin{align} x^2 \equiv 1 \pmod p \\ x^2 -1 \equiv 0 \pmod p \\ (x-1) (x+1) \equiv 0 \pmod p \end{align} \]又 \(\because\) \(p\) 为奇素数，有且仅有 \(1,p\) 两个因子。
\(\therefore\) 只有两解 \(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} \bmod 341 = 1\) 。由于 \(340\) 是偶数，因此我们检查 \(2^{170}\) ，得到 \(2^{170} \bmod 341=1\) ，满足二次探测定理。同时由于 \(170\) 还是偶数，因此我们进一步检查\(2^{85} \bmod 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\) 个候选人围成一个圈，依次编号为 \(0\dots 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) \bmod 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) \bmod 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} \rfloor\) 。
由于此时起点为 \(8\) ，则等价于：
2 3 4 - 5 6 7 - 0 1因此我们仍然可以从 \(f[8]\) 的结果来推导出 \(f[10]\) 的结果。
我们需要根据 \(f[8]\) 的值进行分类讨论。假设 \(f[8]=s\) ，则根据 \(s\) 和 \(n \bmod k\) 的大小关系有两种情况：
\[\begin{cases} s' = s - n \bmod k + n & (s< n \bmod k) \\ \displaystyle s' = s - n \bmod k + \lfloor \frac{(s - n \bmod k)}{k-1} \rfloor & (s \ge n \bmod 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} \rfloor = 0\) 。
因此当 \(n\) 小于 \(k\) 时，就只能采用第一种递推的算法来计算了。
每次我们将 \(\displaystyle n \to n - \lfloor \frac{n}{k} \rfloor\) ，我们可以近似看做把 \(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 \bmod b)\) 。
这个是比较好证明的：
假设 \(a = k*b+r\) ，有 \(r = a \bmod b\) 。不妨设 \(d\) 为 \(a\) 和 \(b\) 的一个任意一个公约数，则有 \(a \equiv b \equiv 0 \pmod d\) 。
由于同余的性质 \(a-kb \equiv r \equiv 0 \pmod d\) 因此 \(d\) 是 \(b\) 和 \(a \bmod 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 \bot b\) （互质）求 \((x, y)\) 的解。
接下来我们利用前面的朴素欧几里得定律推一波式子。
不难发现此时 \(x\) 变成了 \(y\) ， \(y\) 变成了 \(\displaystyle x - \lfloor \frac{a}{b} \rfloor ~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 \bmod b + b) \bmod 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) \times p ^ {k - 1}\) 。
证明：不难发现所有 \(p\) 的倍数都与 \(x\) 不互质，其他所有数都与它互质。
\(p\) 的倍数刚好有 \(p^{k-1}\) 个（包括了 \(x\) 本身）。
那么就有 \(\varphi(x) = p^k - p^{k-1} = (p - 1) \times p^{k - 1}\) 。
-
若 \(p, q\) 互质，则有 \(\varphi(p \times q) = \varphi(p) \times \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} \times {p_2}^{k_2} \times \cdots \times {p_{n} }^{k_n} = \prod _{i=1}^{n} {p_i}^{k_i}\) 。
\[\displaystyle \varphi(x) = x \times (1 - \frac{1}{p_1}) \times (1 - \frac{1}{p_2}) \times \cdots \times (1 - \frac{1}{p_n}) = x~\prod_{i=1}^{n} (1-\frac{1}{p_i}) \]证明：
我们考虑用前几条定理一起来证明。
\[\begin{align} \varphi(x) &= \prod_{i=1}^{n} \varphi(p_i^{k_i}) \\ &= \prod_{i=1}^{n} (p_i-1)\times {p_i}^{k_i-1}\\ &=\prod_{i=1}^{n} {p_i}^{k_i} \times(1 - \frac{1}{p_i})\\ &=x~ \prod_{i=1}^{n} (1- \frac{1}{p_i}) \end{align} \] -
若 \(p\) 为 \(x\) 的约数（ \(p\) 为质数, \(x\) 为任意正整数），我们有 $\varphi(x \times p) = \varphi(x) \times p $ 。
证明：
我们利用之前的第 \(4\) 个性质来证明就行了。
不难发现 \(\displaystyle \prod _{i=1}^{n} (1 - \frac{1}{p_i})\) 是不会变的，前面的那个 \(x\) 会变成 \(x \times p\) 。
所以乘 \(p\) 就行了。
-
若 \(p\) 不是 \(x\) 的约数（ \(p\) 为质数, \(x\) 为任意正整数），我们有 \(\varphi(x \times p) = \varphi(x) \times (p - 1)\) 。
证明 ：\(p, x\) 互质，由于 \(\varphi\) 积性直接得到。
求欧拉函数
求解单个欧拉函数
我们考虑质因数分解，然后直接利用之前的性质 \(4\) 来求解。
快速分解的话可以用前面讲的 \(Pollard~Rho\) 算法。
求解一系列欧拉函数
首先需要学习 线性筛 ，我们将其改一些地方。
质数 \(p\) 的 \(\varphi(p) = p-1\) 。
对于枚举一个数 \(i\) 和另外一个质数 \(p\) 的积 \(x\) 。
我们在线性筛，把合数筛掉会分两种情况。
- \(p\) 不是 \(i\) 的一个约数，由于性质 \(5\) 就有 \(\varphi(x) = \varphi(i) \times (p - 1)\) 。
- \(p\) 是 \(i\) 的一个约数，此时我们会跳出循环，由于性质 \(6\) 有 \(\varphi(x) = \varphi(i) \times 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 \bmod 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 \times M_i \equiv 1 \pmod {m_i}\) ）。
不难发现这个 \(t_i\) 是一定存在的，因为由于前面要满足两两互质，那么 \(M_i \bot m_i\) ，所以必有逆元。
那么在模 \(M\) 意义下，方程 \((S)\) 有且仅有一个解 ： \(\displaystyle x = (\sum_{i=1}^{n} r_i t_i M_i) \bmod 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 \times 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} \bmod {{p_i}^{k_i}}\) 的值。
然后我们可以考虑用前文提到的 \(CRT\) 来进行合并（需要用扩欧求逆元）。
问题就转化成如何求 \(\displaystyle {n \choose m} \bmod {{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} \rfloor !\) ，那么很容易得到如下递推式：
不难发现这个求单个的复杂度是 \(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} \rfloor\) 的部分递归处理就行了。
然后就可以快速处理出与 \(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\n", ExLucas(n, m));
return 0;
}
