中国剩余定理
内容
考虑形如下列形式的方程组：
当 \(m_1,m_2,\dots,m_n\) 互质时，中国剩余定理可以求出满足该方程组的最小正整数解 \(x\)。
过程与实现
设 \(M=\prod\limits_{i=1}^nm_i\)，\(b_i=\frac{M}{m_i}\)，计算每一个 \(b_i\) 在模 \(m_i\) 意义下的逆元 \(b_i^{-1}\)，则 \(x=\sum\limits_{i=1}^na_ib_ib_i^{-1}\pmod M^{①}\)。
void CRT(int a[],int m[],int n){
int M=1,res,b[N],invb[N];
for(int i=1;i<=n;i++) M*=m[i];
for(int i=1;i<=n;i++){
b[i]=M/m[i];
invb[i]=inv(b[i],m[i]);
}
for(int i=1;i<=n;i++)
res=(res+a[i]*b[i]*invb[i])%M;
return res;
}
证明
首先，对于任意 \(j\not =i\)，均有 \(b_j\equiv 0\pmod {m_i}\)，则 \(b_jb_j^{-1}\equiv 0\pmod {m_i}\)，考虑到 \(b_ib_i^{-1}\equiv 1\pmod {m_i}^②\)，故有：
故对于任意 \(i\)，均有 \(x\equiv a_i\pmod {m_i}\)，得证。
扩展中国剩余定理
当 \(m_1,m_2,\dots,m_n\) 不互质时，可以用扩展中国剩余定理计算。
首先考虑两个方程的情况：
设这两个方程分别为 \(x\equiv a_1\pmod {m_1},x\equiv a_2\pmod {m_2}\)，将其改写为不定方程的形式：\(x=pm_1+a_1,x=qm_2+a_2\)，作差，得 \(pm_1-qm_2=a_2-a_1\)，根据裴蜀定理易知，当 \(\gcd(m_1,m_2)\not \mid (a2-a1)\) 时方程无解，否则可以通过扩展欧几里得算法求出一组解 \((p,q)\)，代入 \(x=pm_1+a_1\) 得解 \(x\equiv pm_1+a_1\pmod {\text{lcm}(m1,m2)}\)，这是一个新的方程。
当存在多个方程时，用上面的方法逐一合并即可。
int mode(int x,int p){
return (x%p+p)%p;
}
pair<int,int> merge(int a1,int m1,int a2,int m2){
int x=0,y=0,nx,ny;
int m=mode(a2-a1,m2);
int gcd=exgcd(m1,m2,x,y);
if((a2-a1)%gcd) return {-1,-1};
x=mode(x*(m/gcd),m2/gcd);
nx=a1+m1*x;
ny=m1/gcd*m2;
nx=mode(nx,ny);
return {nx,ny};
}
int exCRT(int a[],int m[],int n){
int ans=a[1],M=m[1];
for(int i=2;i<=n;i++){
pair<int,int> res=merge(ans,M,a[i],m[i]);
ans=res.first,M=res.second;
if(ans==-1) return -1;
}
return ans;
}
注解
①：\(b_i\) 和 \(b_i^{-1}\) 相乘时不需要取模。
②：\(b_i^{-1}\) 可以理解为 \(b_i\) 在模 \(m_i\) 意义下的倒数，两数相乘显然为 \(1\)。
卢卡斯定理
内容
对于质数 \(p\)，有
实现
观察上述式子，左半部分可以迭代求解，右半部分可以通过预处理直接求解，迭代的边界：当 \(m=0\) 时返回 \(1\)。
时间复杂度为 \(O(p+T\log_{p} n)\)，其中 \(T\) 为询问次数。
int Lucas(int n,int m,int p){
if(!m) return 1;
return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p;
}
证明
- 引理 1：
- 证明：
不难发现，\({p\choose n}\bmod p=\frac{p!}{n!(p-n)!}\bmod p\)，当 \(n=0\) 或 \(n=p\) 时上下才能约掉 \(p\)，故 \({p\choose n}\bmod p=[n=0\lor n=p]\)。
则有：
证毕。
考虑二项式 \((1+x)^n\)，有：
对比两边系数 \(x^m\) 的系数，得
其中，有 \(ip+j=m,j<p\)，故 \(i=\lfloor m/p\rfloor,j=m\bmod p\)，故
得证。
扩展卢卡斯定理
过程与证明
当模数不是质数时，就需要用到扩展卢卡斯定理。
扩展卢卡斯定理用于求解 \({n \choose m}\bmod p\)，其中 \(p\) 可能为合数。
- 第一部分：
设 \(p\) 的唯一分解为 \(\prod\limits_{i=1}^kp_i^{\alpha_i}\)，构造如下的同余方程组：
容易发现，只需要求出 \(b_1\sim b_k\) 就可以使用中国剩余定理求出 \({n\choose m}\bmod p^①\)。
- 第二部分：
根据定义，\(b_i={n\choose m}\bmod {p_i^{\alpha_i}}\)，我们就将问题转化为了求 \({n\choose m}\bmod p^{\alpha}\) 的值，其中 \(p\) 是质数。
将组合数展开，得到 \(\frac{n!}{m!(n-m)!}\bmod p^{\alpha}\)，直接计算需要求出分母的逆元，但考虑到逆元可能不存在的情况\(^②\)，将其改写为：
其中，\(x,y,z\) 分别是 \(n!,m!,(n-m)!\) 中因子 \(p\) 的数量。
- 第三部分：
问题转化为了求
的值，先考虑计算 \(n!\bmod p^k\)，根据威尔逊定理的推论\(^③\)，有：
故
则有
其中，\(x,x'\) 分别是 \(n!,\Big(\Big\lfloor\frac{n}{q}\Big\rfloor\Big)\) 中因子 \(p\) 的数量。
式子左半部分可以递归解决，中间和右边的部分利用威尔逊定理的推论解决。
实现
int solve(int n,int x,int P){
if(!n) return 1;
int sum=1;
for(int i=1;i<=P;i++)
if(i%x) sum=sum*i%P;
sum=q_pow(sum,n/P,P);
for(int i=1;i<=n%P;i++)
if(i%x) sum=sum*i%P;
return sum*solve(n/x,x,P)%P;
}
int Crt(int x,int p,int mod){
return inv(p/mod,mod)*(p/mod)*x;
}
int calc(int n,int m,int x,int P){
int res1=solve(n,x,P),res2=solve(m,x,P);
int res3=solve(n-m,x,P),sum=0;
for(int i=n;i;i/=x) sum+=i/x;
for(int i=m;i;i/=x) sum-=i/x;
for(int i=n-m;i;i/=x) sum-=i/x;
return ((q_pow(x,sum,P)*res1%P)*(inv(res2,P)*inv(res3,P)%P))%P;
}
int exLucas(int n,int m,int P){
int ans=0,t=P;
for(int i=2;i*i<=t;i++)
if(t%i==0){
int k=1;
while(t%i==0) k*=i,t/=i;
ans=(ans+Crt(calc(n,m,i,k),P,k))%P;
}
if(t>1) ans=(ans+Crt(calc(n,m,t,t),P,t))%P;
return ans;
}
注解
①：\(p_1^{\alpha_1},p_2^{\alpha_2},\dots,p_k^{\alpha_k}\) 均互质，其最小公倍数就是 \(p\)。
②：根据裴蜀定理，逆元存在的充要条件是 \(\gcd(m!,p^{\alpha})=1\)，不一定成立。
③：详见威尔逊定理。
快速 GCD 与光速幂
快速 GCD
一种基于值域预处理的 \(O(V)\sim O(1)\) GCD 算法。
大致思路是将 \(1\sim V\) 的每个数 \(x\) 分解为三个小于 \(\sqrt x\) 的数，并打一个 \(\sqrt V\times\sqrt V\) 的表后做到 \(O(1)\) 回答。
详见快速 GCD。
void sieve(){
fac[1][0]=fac[1][1]=fac[1][2]=1;
for(int i=2;i<=L;i++){
if(!v[i]){
prime[++cnt]=i;
fac[i][0]=fac[i][1]=1;
fac[i][2]=i;
}
for(int j=1;j<=cnt&&i*prime[j]<=L;j++){
int x=i*prime[j];v[x]=1;
fac[x][0]=fac[i][0]*prime[j];
fac[x][1]=fac[i][1];
fac[x][2]=fac[i][2];
if(fac[x][0]>fac[x][1]) swap(fac[x][0],fac[x][1]);
if(fac[x][1]>fac[x][2]) swap(fac[x][1],fac[x][2]);
if(i%prime[j]==0) break;
}
}
}
void init(){
sieve();
for(int i=0;i<=P;i++)
pre[0][i]=pre[i][0]=i;
for(int i=1;i<=P;i++)
for(int j=1;j<=i;j++)
pre[i][j]=pre[j][i]=pre[j][i%j];
}
int gcd(int a,int b){
int ans=1;
for(int i=0;i<3;i++){
int tmp;
if(fac[a][i]>P){
if(b%fac[a][i]) tmp=1;
else tmp=fac[a][i];
}
else tmp=pre[fac[a][i]][b%fac[a][i]];
b/=tmp;
ans*=tmp;
}
return ans;
}
光速幂
一种在底数固定的情况下的 \(O(\sqrt t)\sim O(1)\) 的求 \(a^b\) 的算法，其中，\(t=\max b\)。
当底数固定时，可以 \(O(\sqrt t)\) 预处理出 \(a^0,a^1,\dots,a^{\sqrt t-1}\) 和 \(a^{\lfloor\frac{b}{\sqrt t}\rfloor\sqrt t},a^{2\sqrt t},\dots,a^t\)，当询问 \(a^b\) 时可以 \(O(1)\) 得出答案为 \(a^{\sqrt t}\times a^{b\bmod \sqrt t}\)。本质上相当于将 \(b\) 进行 \(\sqrt t\) 进制分解。
当然，可以不按 \(\sqrt t\) 进行分块，例如，当按 \(\sqrt[3]t\) 分块时，可以得到 \(O(\sqrt[3]t)\sim O(1)\) 的光速幂，但常数会变大。
具体的说，当按 \(\sqrt[k]t\) 分块时，其时间复杂度为 \(O(k\sqrt[k]t)\sim O(k)\)。
设询问次数为 \(q\)，则总时间复杂度为 \(O(k\sqrt[k]t+qk)\)，可以根据 \(q,t\) 的大小酌情选择 \(k\)，当 \(\sqrt t\) 可以被接受时可以将 \(k\) 设为 \(2\)，当 \(t\) 很大时可以增大 \(k\) 以此减小预处理复杂度，但这样同时也会增大询问的复杂度。比如，当 \(t=10^{18},q=10^6\) 时 \(k\) 的最佳取值大约是 \(3.6\)。
这里给出 \(k=2\) 的实现。
void init(int x){
B=sqrt(V)+1;fac[0]=facs[0]=1;
for(int i=1;i<=B;i++) fac[i]=fac[i-1]*x%mod;
for(int i=1;i<=B;i++) facs[i]=facs[i-1]*fac[B]%mod;
}
int calc(int x){return fas[x/B]*fac[x%B]%mod;}