同余系基本定理2
command_block
2020-04-04 12:57:04
上回请看 [同余系基本定理1](https://www.luogu.com.cn/blog/command-block/tong-yu-xi1)
默认大家已经对同余系有一定的了解,而且能够熟练运用基础算法。
- ## [P1495 【模板】中国剩余定理(CRT)](https://www.luogu.com.cn/problem/P1495)
有若干个同余方程$\begin{cases}x=c_1\pmod {p_1}\\x=c_2\pmod {p_2}\\...\\x=c_m\pmod {p_m}\end{cases}$
满足$p_{1...m}$两两互质,求$x$的最小正整数解。
考虑把这些方程两两**合并**。
先看两个方程$\begin{cases}x=c_1\pmod{p_1}\\x=c_2\pmod {p_2}\end{cases}$
可以构造$x=c_1p_2+c_2p_1$,这样在模$p_1$时显出$c_1p_2$,模$p_2$时显出$c_2p_1$.
问题在于,我们希望得到$c_1$而非$c_1p_2$。这好办,再乘上$p_2$在模$p_1$时的逆即可。
合并完之后的模数是$p_1p_2$,由于互质所以总能求逆。
边界方程可以视作$x=0\pmod 1$.
复杂度$O(m\log p)$.
```cpp
#include<cstdio>
#define ll long long
using namespace std;
void exgcd(ll a,ll b,ll &x,ll &y){
if (b==0){x=1;y=0;return ;}
exgcd(b,a%b,y,x);y-=(a/b)*x;
}
ll inv(ll a,ll m){
ll x,y;exgcd(a,m,x,y);
return (x%m+m)%m;
}
int n;
ll c,p,ret,mp;
int main()
{
scanf("%d",&n);
ret=0;mp=1;
for (int i=1;i<=n;i++){
scanf("%lld%lld",&p,&c);
ll sp=mp*p;
ret=(ret*inv(p,mp)%sp*p+mp*inv(mp,p)%sp*c)%sp;
mp=sp;
}printf("%lld",ret);
return 0;
}
```
- **例题** : [P3868 [TJOI2009]猜数字](https://www.luogu.com.cn/problem/P3868)
- ## [P3807 【模板】卢卡斯定理](https://www.luogu.com.cn/problem/P3807)
$\color{blue}\text{定理:}$ 当$p$是**素数**时
$$\dbinom{n}{m}\bmod p=\dbinom{\lfloor n/p\rfloor}{\lfloor m/p\rfloor}\dbinom{n\bmod p}{m\bmod p}\bmod p$$
$\color{blue}\text{证明:}$ (可不掌握)
注意到$[x^m](1+x)^n=\dbinom{n}{m}$
构造$(1+x)^p=1+\binom{p}{1}x+\binom{p}{2}x+...+\binom{p}{p}x^p$
- 对于$1\leq i\leq p-1$,总有$\binom{p}{i}=0\pmod p$.
原因是$\dbinom{p}{i}=\dfrac{p!}{(p-i)!i!}$,而$p$是素数,所以 $(p-i)!$ 和 $i!$ 都不含因子$p$,无法抵消$p!$中的恰一个因子$p$.
那么可以得到$(1+x)^p=1+x^p\pmod p$
能进一步推出$(1+x)^{p^k}=(1+x^p)^{p^{k-1}}=...=1+x^{p^k}$
我们把$n,m$都表示成两个**p进制数**,即$\begin{cases}n=\sum\limits_{i=0}a_ip^i\\m=\sum\limits_{i=0}b_ip^i\end{cases}$
然后考虑$(1+x)^n=\prod\limits_{i=0}(1+x)^{a_ip^i}$
$(1+x)^n=\prod\limits_{i=0}(1+x^{p^k})^{a^i}\pmod p$
提取第$m$项系数。
$[x^m](1+x)^n=[x^m]\prod\limits_{i=0}(1+x^{p^k})^{a^i}\pmod p$
由于乘积项的次数**不交**,可以把$m$分解。
“不交”的意思是某两位之间不会互相影响。因为将一个 $p$ 进制数的后 $n$ 位加起来之后肯定不足第 $n+1$ 位的一个 $1$ 大。
$\dbinom{n}{m}=\prod\limits_{i=0}[x^{b_ip^k}](1+x^{p^k})^{a^i}\pmod p$
$\dbinom{n}{m}=\prod\limits_{i=0}[x^{b_i}](1+x)^{a^i}\pmod p$
$\dbinom{n}{m}=\prod\limits_{i=0}\dbinom{a_i}{b_i}\pmod p$
即 : 把 $n,m$ 按照 $p$ 进制拆位,每一位对应求组合数乘积即可。
```cpp
#include<cstdio>
#define MaxN 100500
#define ll long long
using namespace std;
int mod;
ll powM(ll a,int t=mod-2){
ll ret=1;
while(t){
if (t&1)ret=ret*a%mod;
a=a*a%mod;t>>=1;
}return ret;
}
ll fac[MaxN],ifac[MaxN];
ll sC(int n,int m){
if (n<m)return 0;
return fac[n]*ifac[m]%mod*ifac[n-m]%mod;
}
ll C(ll n,ll m){
if (!m)return 1;
return sC(n%mod,m%mod)*C(n/mod,m/mod)%mod;
}
void Init()
{
fac[0]=1;
for (int i=1;i<mod;i++)
fac[i]=fac[i-1]*i%mod;
ifac[mod-1]=powM(fac[mod-1]);
for (int i=mod-1;i;i--)
ifac[i-1]=ifac[i]*i%mod;
}
void solve()
{
ll n,m;
scanf("%lld%lld%d",&n,&m,&mod);
m+=n;Init();
printf("%lld\n",C(m,n));
}
int main()
{
int T;scanf("%d",&T);
while(T--)solve();
return 0;
}
```
- **例题** [P4345 [SHOI2015]超能粒子炮·改](https://www.luogu.com.cn/problem/P4345)
**题意** : 多次给出 $n,k$ ,求 $\sum\limits_{i=0}^k\dbinom{n}{i}\pmod {2333}$
$T\leq 10^5;n,k\leq 10^{18}$,时限$\texttt{1s}$.
设 $p=2333$ ,显然这是个素数。
根据卢卡斯定理有 $\dbinom{n}{m}=\dbinom{\lfloor n/p\rfloor}{\lfloor m/p\rfloor}\dbinom{n\bmod p}{m\bmod p}$
$ANS=\sum\limits_{i=0}^k\dbinom{n}{i}$
$=\sum\limits_{i=0}^k\dbinom{\lfloor n/p\rfloor}{\lfloor i/p\rfloor}\dbinom{n\bmod p}{i\bmod p}$
前半部分是块状变化的,我们分别枚举$\lfloor i/p\rfloor$ 和 $i\bmod p$。注意最后多出来的一个散块。
$=\sum\limits_{i=0}^{\lfloor k/p\rfloor-1}\dbinom{\lfloor n/p\rfloor}{i}\sum\limits_{j=0}^{p-1}\dbinom{n\bmod p}{j}+\dbinom{\lfloor n/p\rfloor}{\lfloor k/p\rfloor}\sum\limits_{j=0}^{k\bmod p}\dbinom{n\bmod p}{j}$
注意到$\sum\limits_{i=0}^{\lfloor k/p\rfloor-1}\dbinom{\lfloor n/p\rfloor}{i},\quad\sum\limits_{j=0}^{p-1}\dbinom{n\bmod p}{j},\quad\sum\limits_{j=0}^{k\bmod p}\dbinom{n\bmod p}{j}$均是子问题。
设 $f(n,k)=\sum\limits_{i=0}^k\dbinom{n}{i}$
则有 $f(n,k)=f(\lfloor n/p\rfloor,\lfloor n/k\rfloor-1)f(n\bmod p,p-1)+\dbinom{\lfloor n/p\rfloor}{\lfloor k/p\rfloor}f(n\bmod p,k\bmod p)$
先预处理 $p\times p$ 的 $f$,这用组合数递推不难做到 $O(p^2)$.
然后就是求解 $\dbinom{\lfloor n/p\rfloor}{\lfloor k/p\rfloor}$ ,使用`Lucas`定理即可做到 $O(\log_p n)$
观察参数 $n$ 的变化 : 不断除以 $p$。
然后由于 $f(n,k)$ 在 $k>n$ 时等于 $f(n,n)$ ,当 $n\leq p$ 时就已经达到边界,只需要递归 $O(\log_p n)$ 次。
复杂度$O(T\log_p^2n+p^2)$.
[提交记录](https://www.luogu.com.cn/record/33866714)
- **周边** : Kummer定理
$\color{blue}\text{定理:}$ $\dbinom{n+m}{m}$ 中素因子 $p$ 的个数 $=n,m$ 在 $p$ 进制下做加法的**进位次数**。
注意到 $n!$ 中 $p$ 的次数是 $\sum\limits_{i=0}\lfloor n/p^i\rfloor$ ,证明比较显然。
根据 $\dbinom{n+m}{m}=\dfrac{(n+m)!}{n!m!}$
则因子 $p$ 的个数是 $\sum\limits_{i=0}\lfloor (n+m)/p^i\rfloor-\lfloor n/p^i\rfloor-\lfloor m/p^i\rfloor$
当某一位上 $\lfloor n/p^i\rfloor$ 相当于在 $p$ 进制表示下去掉后 $i$ 位。
那么只有 $n+m$ 在这里进位了才能逃脱被去掉的命运而恰好多$1$.
- **例题** : [CF582D Number of Binominal Coefficients](https://www.luogu.com.cn/problem/CF582D)
**题意** : 给出参数 $lim,p,a$ ,求有多少个组合数 $\dbinom{n}{k}$ 是 $p^a$ 的倍数,且 $0\leq k \leq n \leq lim$.
$p,a\leq 10^9,lim\leq 10^{1000}$,答案对$10^9+7$取模 ,时限$\texttt{4s}$.
根据 Kummer 定理,能将问题转化成这样的形式:
求两个数 $x,y$ 使得 $x+y\leq lim$ ,而且会在 $p$ 进制下产生多于 $a$ 个进位。
这看起来就很数位`DP`的样子。
设$f[i][j][0/1][0/1]$表示考虑了前$i$位,还需要产生$j$个进位,是否卡上界,这一位是否进位。
注意这里是不需要记录前导 $0$ 的,而且由于进位会影响上一位,还要钦定是否进位。
转移的时候分几类讨论 : (下一位进/不进位)(这一位进/不进位),分别计算可能的数对个数即可。
复杂度$O(\log^2lim)$
```cpp
```
- **中段综合测试** : [P2480 [SDOI2010]古代猪文](https://www.luogu.com.cn/problem/P2480)
**题意** : 给出$n,g$,求$g^{\small\sum\limits_{k|n}\binom{n}{k}}$,对$999911659$取模。
$n,g\leq 10^9$,时限$\texttt{1s}$.
首先,模数是个质数,根据费马小定理,指数要对 $mod-1$ 取模。
则对 $999911658=2\times 3\times 4679\times 35617$ 取模。
观察到没有重复的素因子,这比较和善。
求出分别模 $2,3,4679,35617$ 的结果,然后使用中国剩余定理合并即可(甚至不需要`Exgcd`)。
求组合数可以使用Lucas定理。
复杂度是$O(d(n)\log p+T(2+3+4679+35617))$
- ## 扩展欧拉定理(降幂塔)
上集我们已经证明了 : $[a\perp m]\ \Longrightarrow\ a^{\varphi(m)}=1\pmod m$
现在我们来看扩展形式, $\color{blue}\text{定理:}$
$$a^n=\begin{cases}a^{n\bmod\varphi(m)}&(a\perp m)\\a^n&(a\not\perp m,\ n< \varphi(m))\\a^{(n\bmod \varphi(m))+\varphi(m)}&(a\not\perp m,\ n\geq\varphi(m))\end{cases}\large\pmod m$$
第一种情况已被证明,第二种情况是显然的,现在我们来证明第三种情况。
说人话就是 : 当 $a\not\perp m$ 时,第一个 $\varphi(m)$ 不是循环,从第二个才开始进入。
不妨先考虑 $a=p^k$ 的情况。
此时若有 $p\not\perp m$ ,则 $m$ 可以表示成 $s*p^c$ ,且 $s\perp p$。
那么根据普通欧拉定理有 $p^{\varphi(s)}=1\pmod s$ 。 同时,由于 $\varphi$ 函数的积性,有 $\varphi(s)|\varphi(m)$。
于是有 $p^{\varphi(m)}=1\pmod s$
两边同时乘以 $p^c$ ,能够得到 $p^{\varphi(m)+c}=p^c\pmod {m}$
注意到 $\varphi(m)\geq \varphi(p^c)\geq p^{c-1}(p-1)\geq c$
所以有 $p^c=p^{c+\varphi(m)}=p^{c\bmod\varphi(m)+\varphi(m)}\pmod m$ ,我们证明了对于 $p^c$ 满足结论。
结论同时也对于 $p^{kc}$ 成立,可以归纳证明。
若有 $2c\geq n\geq \varphi(m) \geq c$ ,能推知 $n-c\leq \varphi(m)$
则 $p^n=p^{c}p^{n-c}=p^{c\bmod\varphi(m)+\varphi(m)+(n-c)\bmod\varphi(m)}$
有可能会出现 $c|n$ ,此时已经证毕。否则有 $c\bmod\varphi(m)+(n-c)\bmod\varphi(m)=n\bmod\varphi(m)$
然后对 $n$ 的上界归纳即可。
现在我们对任意的 $p^n$ 证明了 $p^n=p^{(n\bmod\varphi(m))+\varphi(m)}$
然后使用唯一分解定理即可把结论扩展至全体正整数。
- **例题** : [P4139 上帝与集合的正确用法](https://www.luogu.com.cn/problem/P4139)
**题意** : 求 $2^{2^{2^{2^{\dots}}}}\bmod m$ , $\ p\leq 10^7$
设$G(p)=2^{2^{2^{2^{\dots}}}}\bmod m$
根据扩展欧拉定理能得到 $G(p)=2^{G(\varphi(p))+\varphi(p)}$
边界就是 $2^{2^{2^{2^{\dots}}}}\bmod 1=0$
取多少次 $\varphi(p)$ 会得到 $1$ 呢?
注意到欧拉函数的其中一个定义式 : $\varphi(m)=m\prod\limits_{p|m}\frac{p-1}{p}$
当 $m$ 无素因子 $2$ 时,一个奇素因子一定会产生因子 $2$ ,因为 $p-1$ 是偶数。
当 $m$ 有素因子 $2$ 时,值至少减半。
综上,经过 $O(\log p)$ 次递归就能得到 $1$。
线性筛出欧拉函数复杂度是 $O(p+T\log^2p)$ 的。
[提交记录](https://www.luogu.com.cn/record/33847265)
- **附加题** :
- [P4861 按钮](https://www.luogu.com.cn/problem/P4861)
若 $\gcd(k,m)>1$ ,显然无解。
否则满足经典欧拉定理,则有 $k^{\varphi(m)}=1\pmod m$ ,但是 $\varphi(m)$ 并不是答案,可能存在更小的循环节。
可能的循环节一定是 $\varphi(m)$ 的约数,否则能导出矛盾。
于是只需要判定 $\varphi(m)$ 的约数即可,复杂度为 $O(\sqrt{m}\log m)$。
[提交记录] ()
- [P3934 Nephren Ruq Insania](https://www.luogu.com.cn/problem/P3934)
考虑乘方降幂塔,在 $O(\log p)$ 层之后模数就会变为 $1$,所以我们只需要暴力取出 $l$ 后的若干项即可。
上一题中,由于 $2^{2^{2^{2^{\dots}}}}$ 恒大于 $\varphi(m)$ ,所以总可以使用第三种情况。
现在由于可能产生第二种情况,我们计算时需要按照如下原则:
如果算的的数大于 $m$ 则替换成 $(a\bmod m)+m$,否则不变。
[提交记录](https://www.luogu.com.cn/record/33849679)
- [CF906D Power Tower](https://www.luogu.com.cn/problem/CF906D) : 上一题的静态版,但数据较强。
- [P3747 [六省联考2017]相逢是问候](https://www.luogu.com.cn/problem/P3747)
仍然考虑结论 : 降幂塔不超过 $O(\log p)$ 层。某个位置被赋值 $O(\log p)$ 次之后就必然变为 $c^{c^{c^{\dots}}}\bmod p$。
考虑把不再变化的一整个区间冻结,可以证明只会产生 $O(n\log p)$ 次修改叶子的操作。
每次单点修改需要计算一次降幂塔,朴素实现复杂度是$O(\log p^2)$的,无法通过。
模数只有 $O(\log p)$ 个,可以分别于处理 $c$ 的光速幂。
总复杂度为$O(\sqrt{p}\log p+m\log ^2p)$。
[提交记录](https://www.luogu.com.cn/record/33862724)
- [P5240 Derivation](https://www.luogu.com.cn/problem/P5240?contestId=14994)
- ## [P4777 【模板】扩展中国剩余定理（EXCRT）](https://www.luogu.com.cn/problem/P4777)
仍然是若干个方程组$\begin{cases}x=c_1\pmod {p_1}\\x=c_2\pmod {p_2}\\...\\x=c_m\pmod {p_m}\end{cases}$
求$x$的最小解,但不保证$p_{1...m}$互质。
仍然考虑两两合并$\begin{cases}x=c_1\pmod{p_1}\\x=c_2\pmod {p_2}\end{cases}$
原先我们要直接构造某个模数在另一个模数下的逆,现在似乎不太行了。
第一个方程的通解为$c_1+kp_1$,现在就是求一个 $k$ 使得$c_1+kp_1=c_2\pmod{p_2}$
移项有$kp_1=c_2-c_1\pmod{p_2}$,熟练的同学已经知道怎么做了。
等价于$kp_1+yp_2=c_2-c_1$
根据斐蜀定理,方程有解的充要条件是$\gcd(p_1,p_2)|(c_2-c_1)$.
如果可解,先扩欧求出$kp_1+yp_2=\gcd(p_1,p_2)$,然后乘以$\dfrac{c_2-c_1}{\gcd(p_1,p_2)}$即可。
注意,新的模数是$lcm(p_1,p_2)$。此外还要注意防止溢出。
复杂度仍然是$O(m\log p)$.
```cpp
#include<cstdio>
#define ll long long
using namespace std;
inline ll mul(ll a,ll b,ll m){
ll d=((long double)a/m*b+0.5);
ll r=a*b-d*m;
return r<0?r+m:r;
}
ll gcd(ll a,ll b)
{return !b ? a : gcd(b,a%b);}
void exgcd(ll a,ll b,ll &x,ll &y){
if (b==0){x=1;y=0;return ;}
exgcd(b,a%b,y,x);y-=(a/b)*x;
}
int n;
ll c1,p1,c2,p2;
int main()
{
scanf("%d",&n);
c1=0;p1=1;
for (int i=1;i<=n;i++){
scanf("%lld%lld",&p2,&c2);
ll pd=gcd(p1,p2),sp=p2/pd*p1,k,y;
c2=(c2-c1%p2+p2)%p2;
//if (c2%pd) No Sol;
exgcd(p1,p2,k,y);
k=mul(k,(c2/pd),p2);
c1=(c1+mul(p1,k,sp))%sp;
p1=sp;
}printf("%lld",c1);
return 0;
}
```
- **例题** : [P4774 [NOI2018]屠龙勇士](https://www.luogu.com.cn/problem/P4774)
发现不会做这道题是写本文的动力之一。
首先容易发现,屠龙的规则和顺序严格确定,那么如果我们能成功,每一轮使用的剑是相同的。拿`multiset`维护一下便可。
设杀第$i$条龙的剑攻击力为$k_i$.
击杀每条龙的条件很像取模,我们能列出式子$k_ix=a_i\pmod {p_i}$
问题在于,可能有$a_i>p_i$,此时可能在比$0$大但为$p_i$倍数的血量停下来,所以$x$满足上述方程并不是充要条件。
考虑对每条龙记录砍到负数的最小次数的最大值,最后$x$在通解中取大一点就好了。
接下来的问题就是解方程组$\begin{cases}k_1x=c_1\pmod {p_1}\\k_2x=c_2\pmod {p_2}\\...\\k_mx=c_m\pmod {p_m}\end{cases}$
考虑将$kx=c\pmod p$变为形如$x=c\pmod p$的经典形式。
先特判$0$ : 当$p|k$且$p|c$时方程恒成立,不考虑。当$p|k$但$p\!\!\not|\ c$时无解。
由于$k,p$不一定互质,我们不能直接通过逆元来转化。
不过,我们仍然能尝试求解$kx+py=c$,无解则原方程组无解。
否则由斐蜀定理得$\gcd(k,p)|c$,设$d=\gcd(k,p)$。
则有$\dfrac{k}{d}x+\dfrac{p}{d}y=\dfrac{c}{d}$.
我们对$\dfrac{p}{d}$取模即得到$\dfrac{k}{d}x=\dfrac{c}{d}\pmod {\dfrac{p}{d}}$
由于$\dfrac{k}{d}\perp\dfrac{p}{d}$,我们就可以求逆元变为经典形式了。
剩下的就是一个扩展中国剩余定理。
```cpp
#include<algorithm>
#include<cstdio>
#include<set>
#define Itor set<ll>::iterator
#define ll long long
#define MaxN 100500
using namespace std;
const int mod=998244353;
ll mul(ll a,ll b,ll m)
{return (((a>>20)*b%m<<20)+(a-(a>>20<<20))*b)%m;}
ll gcd(ll a,ll b)
{return !b ? a : gcd(b,a%b);}
ll lcm(ll a,ll b)
{return a/gcd(a,b)*b;}
void exgcd(ll a,ll b,ll &x,ll &y){
if (!b){x=1;y=0;return ;}
exgcd(b,a%b,y,x);y-=(a/b)*x;
}
ll inv(ll a,ll m){
ll x,y;exgcd(a,m,x,y);
return (x+m)%m;
}
bool fl;
void merge(ll &c1,ll &p1,ll c2,ll p2)
{
ll c=(c2-c1%p2+p2)%p2,d=gcd(p1,p2);
if (c%d){puts("-1");fl=1;return ;}
ll k,y;exgcd(p1,p2,k,y);
ll p=lcm(p1,p2);
c1=(c1+mul(mul(k,p1,p),c/d,p))%p;
p1=p;
}
multiset<ll> s;
ll a[MaxN],k[MaxN],p[MaxN],t[MaxN]
,mx,tc,tp;
void calc(ll k,ll c,ll p)
{
mx=max(mx,(c-1)/k+1);
if (k%p==0){
if (c%p==0)return ;
else {puts("-1");fl=1;return ;}
}
ll d=gcd(k,p);
if (c%d){puts("-1");fl=1;return ;}
k/=d;p/=d;c/=d;
c=mul(c,inv(k,p),p);
merge(tc,tp,c,p);
}
int n,m;
void solve()
{
scanf("%d%d",&n,&m);
for (int i=1;i<=n;i++)scanf("%lld",&a[i]);
for (int i=1;i<=n;i++)scanf("%lld",&p[i]);
for (int i=1;i<=n;i++)scanf("%lld",&t[i]);
s.clear();
for (int i=1;i<=m;i++){
ll x;scanf("%lld",&x);
s.insert(x);
}
for (int i=1;i<=n;i++){
Itor it=s.upper_bound(a[i]);
if (it!=s.begin())it--;
k[i]=*it;s.erase(it);
s.insert(t[i]);
}
mx=0;tp=1;tc=0;fl=0;
for (int i=1;i<=n&&!fl;i++)
calc(k[i],a[i],p[i]);
if (fl)return ;
ll tk=mx/tp;
while(tk*tp+tc<mx)tk++;
printf("%lld\n",tk*tp+tc);
}
int main()
{
int T;scanf("%d",&T);
while(T--)solve();
return 0;
}
```
- ## [P4720 【模板】扩展卢卡斯](https://www.luogu.com.cn/problem/P4720)
跟普通卢卡斯定理关系不大 (
我们把模数分解成$\prod p_i^{c_i}$的形式,对于每个$p^c$单独求解,然后使用普通中国剩余定理合并即可。
考虑$\dbinom{n}{m}=\dfrac{n!}{m!(n-m)!}\pmod{p^c}$
但这毫无作用,往往有 $p|(n!)$ ,这样就无法求逆了。
考虑将$p$的幂次提尽,设$v(n)$为$n$中含有$p$的幂次。
则有$\dfrac{n!}{p^{v(n)}}≠0\pmod{p^c}$,这样就可以求逆了。设$r(n)=\dfrac{n!}{p^{v(n)}}$
我们计算$\dfrac{r(n)}{r(m)r(n-m)}p^{v(n)-v(m)-v(n-m)}$即可得到$\dbinom{n}{m}$。
上文介绍了$v(n)=\sum\limits_{i=0}\lfloor n/p^i\rfloor$ ,现在问题变为求$r(n)$.
对于$n!=1*2*3*...*n$,我们可以把其中为$p$的倍数的项提出来。
变为$\Big[1*2*...*(p-1)\Big]*p*\Big[(p+1)(p+2)...(p+p-1)\Big]*2p*...$
对于中间的每一块,$\prod_{i=1}^{p-1}(kp+i)\pmod {p^c}$,可以变为$[0,p^c)$中的某一段连乘积,维护阶乘及其逆元即可。
这样的段会有$O(n/p)$个,但是由于$\bmod\ p^c$,我们可以批量计算。
具体见代码,注意散块的边界情况。这样计算一次的复杂度是$O(p^c)$的。
对于那些为$p$的倍数的位置,共$\lfloor n/p\rfloor$个,统一除以$p$之后,又变成了$1*2*3*...$
于是就变为求解$r(\lfloor n/p\rfloor)$.
预处理复杂度为$O(\sum p_i^{c_i})$,回答一次的复杂度为$O(\log^2 p)$.
( 如果采用欧拉定理+光速幂可以达到$O(\sqrt{p}+\sum p_i^{c_i})$预处理$O(\log p)$回答 )
```cpp
#include<algorithm>
#include<cstdio>
#define MaxN 1005000
#define ll long long
using namespace std;
void exgcd(ll a,ll b,ll &x,ll &y){
if (b==0){x=1;y=0;return ;}
exgcd(b,a%b,y,x);y-=(a/b)*x;
}
ll inv(ll a,ll m){
ll x,y;exgcd(a,m,x,y);
return (x%m+m)%m;
}
ll powM(ll a,ll t,int mod)
{
ll ret=1;
while(t){
if (t&1)ret=ret*a%mod;
a=a*a%mod;t>>=1;
}return ret;
}
ll v(ll n,int p){
ll ret=0;
while(n)ret+=(n/=p);
return ret;
}
ll sav[MaxN];
ll r(ll n,int p,int mod){
if (n==0)return 1;
return powM(sav[mod-1],n/mod,mod)*sav[n%mod]%mod*r(n/p,p,mod)%mod;
}
ll C(ll n,ll m,int p,int mod)
{
int c=v(n,p)-v(m,p)-v(n-m,p);
ll ans=powM(p,c,mod);
if (!ans)return 0;
sav[0]=1;
for (int i=1;i<mod;i++)
if (i%p==0)sav[i]=sav[i-1];
else sav[i]=sav[i-1]*i%mod;
ans=ans*r(n,p,mod)%mod;
ans=ans*inv(r(m,p,mod),mod)%mod;
ans=ans*inv(r(n-m,p,mod),mod)%mod;
return ans;
}
ll ret,mp=1;
void add(ll c,ll p)
{
ll sp=mp*p;
ret=(ret*inv(p,mp)%sp*p+mp*inv(mp,p)%sp*c)%sp;
mp=sp;
}
ll n,m;int p;
int main()
{
scanf("%lld%lld%lld",&n,&m,&p);
int d=2;
while(d*d<=p){
if (p%d==0){
int mod=1;
while(p%d==0){
p/=d;
mod*=d;
}add(C(n,m,d,mod),mod);
}d++;
}if (p>1)add(C(n,m,p,p),p);
printf("%lld",ret);
return 0;
}
```
- **例题①** : [P2183 [国家集训队]礼物](https://www.luogu.com.cn/problem/P2183)
**题意** : 有$n$个不同的小球,投入$m$个盒子里面,钦定第$i$个盒子要投$w_i$个球,问方案数。
对 $mod$ 取模,保证 $mod\leq 10^9$ 且分解之后的每个 $p^c\leq 10^5$
$m\leq 5,n\leq 10^9$ ,时限$\texttt{1s}$.
先判掉无解的情况。设$S[k]=\sum\limits_{i=1}^kw_i$
不难发现题目叫我们求的是$\prod\limits_{i=1}^m\dbinom{n-S[i-1]}{w_i}$
也就是在固定模数下求解$m$次组合数,套用`EXLucas`即可。
理论复杂度为$O(\sqrt{p}+\sum p^c+m\log^2p)$,但是$O\Big(m(\sqrt{p}+\sum p^c+\log^2p)\Big)$也能过我就偷懒了。
[评测记录](https://www.luogu.com.cn/record/32493149)
- **例题②** : [P3726 [AH2017/HNOI2017]抛硬币](https://www.luogu.com.cn/problem/P3726)
[题解Link](https://www.luogu.com.cn/blog/command-block/solution-p3726)