数论学习笔记

Post time: 2020-11-14 12:11:49

一、裴蜀定理

结论：

\(ax+by=c\) 有整数解，等价于 \(\gcd(a,b) | c\)

扩展一下：

\(a_1x_1+a_2x_2+...+a_nx_n=c\) 有整数解，等价于\(\gcd(a_1,a_2,...,a_n)|c\)

证明：略

二、欧几里得定理

\[\gcd(a,b)=\left\{ \begin{aligned} & a,b=0 \\ & \gcd(b,a\bmod b),b \neq 0 \end{aligned} \right. \]

证明：略

三、扩展欧几里得

首先，由欧几里得定理得：

\[ax+by=\gcd(a,b)=a'x'+b'y' \]

其中

\[\left\{ \begin{aligned} & a'=b \\ & b'=a\bmod b \end{aligned} \right. \]

那么

\[\begin{aligned} ax+by&=bx'+(a-\lfloor{a/b}\rfloor \times b)y'\\ &=bx'+ay'+b\times \lfloor{a/b}\rfloor\times y'\\ &=ay'+b(x'-\lfloor{a/b}\rfloor\times y')\\ \end{aligned} \]

得

\[\left\{ \begin{aligned} & x=y' \\ & y=x'-\lfloor{a/b}\rfloor\times y' \end{aligned} \right. \]

又知当 \(b=0\) 时，\(\gcd(a,b)=a\)，此时显然有一组解为

\[\left\{ \begin{aligned} & x=1 \\ & y=0 \\ \end{aligned} \right. \]

故我们可以 \(\gcd\) 中更新 \(x,y\) 的值。证毕。

对于一个形如 \(ax+by=c\) 的不定方程，若我们已经通过 \(Exgcd\) 求出了一组特解 \(x_0,y_0\)，可以得到此方程的整数解通解为：

\[\left\{ \begin{aligned} & x_t=x_0+b't \\ & y_t=y_0-a't \\ \end{aligned} \right. \]

其中

\[\left\{ \begin{aligned} & a'=\frac{a}{g} \\ & b'=\frac{b}{g} \\ \end{aligned} \right. ,g=\gcd(a,b) \]

模板题P5656

点击查看代码

ll exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){x=1,y=0;return a;}
	ll d=exgcd(b,a%b,y,x);
	y-=(a/b)*x;
	return d;
}

四、中国剩余定理

中国剩余定理用来解形似如下的方程组：

\[\left\{ \begin{aligned} & x\equiv b_1\pmod {a_1} \\ & x\equiv b_2\pmod {a_2} \\ & ...... \\ & x\equiv b_n\pmod {a_n} \\ \end{aligned} \right. \]

其中，对于 \(\forall i,j\in [1,n],\gcd(a_i,a_j)=1\) 。

构造一组解的方式是：

令 \(M=\prod_{i=1}^{n}a_i,M_i=\frac{M}{a_i},M_it_i\equiv 1\pmod {a_i}\)，即 \(t_i\) 为 \(M_i\) 在 \(\bmod a_i\) 意义下的乘法逆元，可由（三）中的 Exgcd 求出。

那么，\(x=\sum_{i=1}^{n}b_iM_it_i\) 就是原方程组的一个解。原方程组的通解为 \(x_0=x+kM,k\in Z\)。原方程组的最小非负整数解为 \(x_0 \bmod M\)。

证明：由定义得，对于 \(\forall i\in [1,n]\)，

对 \(\forall j\in[1,i-1] \bigcup [i+1,n],b_jM_jt_j\equiv 0\pmod {a_i}\);

对 \(i,b_iM_it_i\equiv b_i\pmod {a_i}\)。

故 \(\sum\limits_{j=1}^{n} b_jM_jt_j\equiv b_i\pmod {a_i}\)，即第 \(i\) 个方程组成立。证毕。

模板题P1495

点击查看代码

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const int N=10+13;
int n,a[N],b[N];
ll sum=1,res=0;
inline ll exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){x=1,y=0;return a;}
	ll g=exgcd(b,a%b,y,x);
	y-=(a/b)*x;
	return g;
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i){
		scanf("%d%d",&a[i],&b[i]);
		//if(b[i]==a[i]) ++a[i],b[i]=0;
	}
	for(int i=1;i<=n;++i) sum*=a[i];
	for(int i=1;i<=n;++i){
		ll x,y,tmp=sum/a[i];
		ll g=exgcd(tmp,a[i],x,y);
		x=(x%a[i]+a[i])%a[i];
		res+=b[i]*tmp*x;
	}
	printf("%lld\n",res%sum);
	return 0;
}

五、扩展中国剩余定理

扩展中国剩余定理用来解形似如下的方程组：

\[\left\{ \begin{aligned} & x\equiv b_1\pmod {a_1} \\ & x\equiv b_2\pmod {a_2} \\ & ...... \\ & x\equiv b_n\pmod {a_n} \\ \end{aligned} \right. \]

其中，\(a_1,a_2,...,a_n\) 不一定两两互质，\(n\leq 10^5\)。

此时我们会发现，由于不一定两两互质，导致上面中国剩余定理的构造在逆元上出现问题（可能没有逆元），从而导致算法崩溃。

于是我们看到 \(n\leq 10^5\)，想可不可以通过最坏 \(O(\log n)\) 复杂度合并两个同余方程。

我们先考虑这样一个问题，对于

\[\left\{ \begin{aligned} & x\equiv r_1\pmod {m_1} \\ & x\equiv r_2\pmod {m_2} \\ \end{aligned} \right. \]

这个同余方程组的合并。我们可以将其转化成这样一个形式，即

\[x=k_1m_1+r_1=k_2m_2+r_2 \]

移项，得

\[k_1m_1-k_2m_2=r_2-r_1 \]

其中，只有 \(k_1\) 和 \(k_2\) 是变量，转化成了一个不定方程。首先我们可以用裴蜀定理判断有没有根，即令 \(d=\gcd(m_1,m_2)\)，只需判断 \(d|(r_2-r_1)\) 的真假即可。

令 \(p_1=\frac{m_1}{d},p_2=\frac{m_2}{d}\)，则原式化为

\[k_1p_1-k_2p_2=\frac{r_2-r_1}{d} \]

考虑解这样一个方程

\[t_1p_1-t_2p_2=1 \]

可得

\[\left\{ \begin{aligned} & k_1=\frac{r_2-r_1}{d}t_1 \\ & k_2=\frac{r_2-r_1}{d}t_2 \\ \end{aligned} \right. \]

那么，原式中 \(x\) 就可以表示出来了，即

\[x=r_1+\frac{r2-r1}{d}t_1m_1 \]

我们设我们构造出来的通解为 \(x'\)，则

\[x\equiv x'\pmod {\mathrm{lcm}(m_1,m_2)} \]

由此可以得到，我们可以合并 \(n-1\) 次同余方程组，最终得到最后一个方程的最小正整数解即为答案。

一个小技巧：龟速乘——用快速幂的思想将乘法转为加法，使得带模乘法不越界。

模板题P4777

点击查看代码

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const int N=1e5+13;
int n;
ll a[N],b[N];
inline ll gcd(ll a,ll b){return !b?a:gcd(b,a%b);}
inline 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;
}
inline ll Mul(ll a,ll n,ll p){
	ll s=0;
	while(n){
		if(n&1) s=(s+a)%p;
		a=(a+a)%p;
		n>>=1;
	}
	return s;
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i) scanf("%lld%lld",&a[i],&b[i]);
	for(int i=2;i<=n;++i){
		ll g=gcd(a[i-1],a[i]),x,y;
		exgcd(a[i-1]/g,a[i]/g,x,y);
		a[i]*=a[i-1]/g;
		b[i]=(Mul(Mul((b[i]-b[i-1])/g,(x%a[i]+a[i])%a[i],a[i]),a[i-1],a[i])+b[i-1])%a[i];
		b[i]=(b[i]+a[i])%a[i];
	}
	printf("%lld\n",b[n]);
	return 0;
}

练习题P3868

点击查看代码

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const int N=10+13;
int n;
ll a[N],b[N];
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){x=1,y=0;return;}
	exgcd(b,a%b,y,x);
	y-=(a/b)*x;
}
inline ll Winter_mul(ll a,ll k,ll p){
	ll s=0;
	while(k){
		if(k&1) s=(s+a)%p;
		a=(a+a)%p;
		k>>=1;
	}
	return (s+p)%p;
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i) scanf("%lld",&b[i]);
	for(int i=1;i<=n;++i) scanf("%lld",&a[i]);
	for(int i=1;i<=n;++i) b[i]%=a[i];
	for(int i=2;i<=n;++i){
		ll g=gcd(a[i-1],a[i]),x,y;
		exgcd(a[i-1]/g,a[i]/g,x,y);
		a[i]*=a[i-1]/g;
		b[i]=b[i-1]+Winter_mul(Winter_mul((b[i]-b[i-1])/g,(x%a[i]+a[i])%a[i],a[i]),(a[i-1]%a[i]+a[i])%a[i],a[i]);
	}
	printf("%lld\n",b[n]);
	return 0;
}

六、扩展欧拉定理

前置： \(O(\sqrt p)\) 求 \(\varphi(p)\) 小公式

\[\varphi(p)=p\ \cdot (1-\frac{1}{p_1})\ \cdot (1-\frac{1}{p_2})\ \cdot ...\ \cdot (1-\frac{1}{p_n}) \]

其中，\(p=p_1^{k_1}\cdot p_2^{k_2}\cdot ...\cdot p_n^{k_n}\)。

欧拉定理：

由

\[a^{\varphi(m)}\equiv 1\pmod m \]

可得

\[a^b\equiv a^{b\bmod \varphi(m)}\pmod m,\gcd(a,m)=1 \]

扩展欧拉定理：

\[a^b\equiv\left\{ \begin{aligned} & a^b,b<\varphi(m) \\ & a^{(b\bmod\varphi(m))+\varphi(m)},b\geq\varphi(m) \end{aligned} \right.\pmod m,\ a,m\in Z \]

模板题P5091

点击查看代码

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
inline int phi(int x){
	int ans=x,n=sqrt(x);
	for(int i=2;i<=n;++i){
		if(x%i==0){
			ans=ans/i*(i-1);
			while(x%i==0) x/=i;
		}
	}
	if(x>1) ans=ans/x*(x-1);
	return ans;
}
inline int read(int m){
	int res=0,flag=0;char c=getchar();
	while(!isdigit(c)) c=getchar();
	while(isdigit(c)){
		res=res*10+c-'0';
		if(res>=m) flag=1,res%=m;
		c=getchar();
	}
	if(res>=m) flag=1,res%=m;
	return flag?res+m:res;
}
inline int quick(int a,int b,int p){
	int s=1;
	while(b){
		if(b&1) s=(1ll*s*a)%p;
		a=(1ll*a*a)%p;
		b>>=1;
	}
	return s;
}
int main(){
	int a,b,m;
	scanf("%d%d",&a,&m);
	b=read(phi(m));
	printf("%d\n",quick(a,b,m));
	return 0;
}

七、高斯消元

高斯消元可以在 \(O(n^3)\) 复杂度内解决 \(n\) 元一次方程组的问题。解决过程如下：

1. 把方程组化为一个 \(n\times (n+1)\) 的矩阵，分别表示 \(x_{1...n}\) 的系数和等号右边的结果。
2. 枚举每一列，即每一个元，找到这个元绝对值最大的系数的列，将这一列换到这一行。注意这里换行是不影响的。
3. 如果绝对值最大的系数仍然为零，那么说明这个元可以取任何数。
4. 进行加减消元，目的是最后令 \(a[i][i]=1\)，且 \(a[i+1...n][i]=0\)。
5. 最终得到一个只有对角线及其右上方有数的矩阵，直接从底向上代入即可。

注意这样一个性质：在消第 \(i\) 个元的时候，对于 \(\forall j\in [1,i-1]\)，都有 \(a[j][1...j-1]=0\).

模板题P3389

点击查看代码

#include<iostream>
#include<cstdio>
using namespace std;
const int N=100+13;
const double eps=1e-6;
inline double fabs(double x){return x<0?-x:x;}
int n;
double a[N][N],ans[N];
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n+1;++j) scanf("%lf",&a[i][j]);
	for(int i=1;i<=n;++i){//枚举每一列
		int r=i;
		for(int j=i+1;j<=n;++j)
			if(fabs(a[r][i])>=fabs(a[j][i])) r=j;
		if(fabs(a[r][i])<eps){puts("No Solution");return 0;}
		if(i!=r) swap(a[i],a[r]);//找到这一列系数最大的值，把它换到对角线上
		double div=a[i][i];//开始加减消元，注意在这里消元的时候已经保证了i左边的对角线左半边全为零。 
		for(int j=i;j<=n+1;++j) a[i][j]/=div;//把a[i][i]消为1 
		for(int j=i+1;j<=n;++j){
			div=a[j][i];
			for(int k=i;k<=n+1;++k) a[j][k]-=a[i][k]*div;//把每一行第i个系数化零 
		}
	}
	ans[n]=a[n][n+1];//对角线右半边还有数，可以从底往上消去 
	for(int i=n-1;i>=1;--i){
		ans[i]=a[i][n+1];
		for(int j=i+1;j<=n;++j) ans[i]-=a[i][j]*ans[j];//减掉这一行右边残留的其他元 
	}
	for(int i=1;i<=n;++i) printf("%.2lf\n",ans[i]);
	return 0;
}

八、线性基

这里的线性基指异或线性基。对于一个数组 \(s[1...n]\)，线性基是一个数组 \(p[1...m]\)，满足 \(p[i]\) 二进制第 \(i\) 位为 \(1\)，且这一位是值为 \(1\) 的最高位。对于 \(\forall i\in [1,n]\)，都 \(\exists j_1,j_2,...,j_k\in [1,m]\ (j_1<j_2<...<j_k)\)，使得 \(p[j_1] \otimes p[j_2]\otimes ...\otimes p_[j_k]=s[i]\)。

线性基的一种构造方法是，考虑对于 \(s\) 中每个数从高位向低位枚举第 \(i\) 位（\(1\) 位），如果此时 \(p_i=0\)，那么直接插入进线性基，否则将这个数异或上 \(p_i\)。感性理解一下即可。

如果需要求最大异或和，只需从高往低贪心即可，如果第 \(i\) 位可以是 \(1\)，那就通过取或不取使其为 \(1\)。

注意：异或的优先级很低，需要加括号！

模板题P3812

点击查看代码

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const int N=50+13;
int n;
ll p[N];
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i){
		ll x;scanf("%lld",&x);
		for(int j=51;j>=0;--j){
			if(x&(1ll<<j)){
				if(!p[j]){p[j]=x;break;}
				else x^=p[j];
			}
		}
	}
	ll ans=0;
	for(int i=51;i>=0;--i)
		if((ans^p[i])>ans) ans^=p[i];
	printf("%lld\n",ans);
	return 0;
}

九、拉格朗日插值

有一个 \(n\) 次多项式，给你 \((n+1)\) 个点，求出这个多项式。

爆力就是 \(O(n^3)\) 做高斯消元。

考虑类似 \(CRT\) 的构造，对于一个点在 \(k\) 处的取值，让构造出的这个多项式在 \(x_i\) 处取到 \(y_i\)，在其他位置都取 \(0\)。这样我们可以得到拉格朗日插值的基本式子：

\[f(x)=\sum_{i=0}^n y_i \prod_{j\neq i} \frac{x-x_j}{x_i-x_j} \]

如果直接做的话，带模复杂度 \(O(n^2 \log n)\)。如果先把分母都乘起来再求逆元，可以做到 \(O(n^2)\)。

模板题P4781

点击查看代码

#include<iostream>
#include<cstdio>
using namespace std;
const int N=2000+13,mod=998244353;
int n,x0,x[N],y[N],sum;
inline int Quick(int a,int k){
	int s=a%mod;--k;
	while(k){
		if(k&1) s=1ll*s*a%mod;
		a=1ll*a*a%mod;
		k>>=1;
	}
	return s%mod;
}
inline int Inv(int k){return Quick(k,mod-2);}
int main(){
	scanf("%d%d",&n,&x0);--n;
	for(int i=0;i<=n;++i) scanf("%d%d",&x[i],&y[i]);
	for(int i=0;i<=n;++i){
		int res1=1,res2=1;
		for(int j=0;j<=n;++j){
			if(i==j) continue;
			res1=1ll*res1*(x0-x[j]+mod)%mod;
			res2=1ll*res2*(x[i]-x[j]+mod)%mod;
		}
		sum=(sum+1ll*y[i]*res1%mod*Inv(res2)%mod)%mod;
	}
	printf("%d\n",sum);
	return 0;
}

十、BSGS（大步小步算法）

BSGS用来解决这样一类问题：

已知 \(a,b,p\)，求最小的 \(x\) 使得 \(a^x\equiv b\pmod p\)，或声明无解。

设 \(x=i\times t-j\)，其中 \(t=\lceil\sqrt p\rceil,1\leq i \leq t,0\leq j\leq t-1\)，则

\[\begin{aligned} a^{i\times t-j}&\equiv b\pmod p\\ (a^i)^t&\equiv b\times a^j\pmod p \end{aligned} \]

将 \(b\times a_j(0\leq j\leq t-1)\) 插入一个哈希表内，然后枚举 \(i=1\sim t\) 判断即可。

时间复杂度是 \(O(\sqrt p)\) 的，哈希表可以使用 unordered_map 实现 \(O(1)\) 的插入和查询。注意 hash 在 c++11 中是关键字，慎用。

模板题P3846

点击查看代码

#include<iostream>
#include<cstdio>
#include<cmath>
#include<unordered_map>
using namespace std;
typedef long long ll;
unordered_map<int,int> ha;
inline int qpow(int a,int k,int p){int s=1;for(;k;k>>=1,a=(ll)a*a%p)if(k&1)s=(ll)s*a%p;return s;}
int main(){
	ha.clear();int p,a,b;
	scanf("%d%d%d",&p,&a,&b);
	int t=sqrt(p)+1;int tmp=1;
	for(int i=0;i<t;++i) ha[(ll)b*tmp%p]=i+1,tmp=(ll)tmp*a%p;
	a=tmp;
	for(int i=1;i<=t;++i){
		tmp=qpow(a,i,p);
		if(ha[tmp]) return printf("%d\n",i*t-ha[tmp]+1),0;
	}
	puts("no solution");
	return 0;
}

十一、min-max 容斥

设全集为 \(U\)，如果对于全集的子集 \(S\)，其所有子集 \(T(T\subseteq S)\)，\(\min(T)\) 都很好求出，但是 \(\max(S)\) 不好求，可以使用 min-max 容斥快速求出。

上式子：

\[\max(S)=\sum_{T\subseteq S}(-1)^{|T|+1}\min(T) \]

将 \(\max\) 和 \(\min\) 互换一下也成立。

十二、Lucas 定理

结论：

\[\binom{n}{m}\equiv \binom{n\bmod p}{m\bmod p}\times\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\pmod{p} \]

证明：

引理 1

对于质数 \(p\) 和任意的整数 \(j\in [1,p)\)，都有

\[\binom{p}{j}\equiv 0\pmod p \]

证明：

\[\binom{p}{j}=\frac{p!}{j!(p-j)!}=\frac{(p-1)!}{j!(p-j)!}\cdot p\equiv 0\pmod p \]

引理 2

对于任意的质数 \(p\)，都有

\[(1+x)^p\equiv 1+x^p\pmod p \]

证明：

\[(1+x)^p=1+x^p+\sum_{i=0}^p \binom{p}{i}x^i\equiv 1+x^p\pmod p \]

设 \(n=k_1p+b_1,m=k_2p+b_2\)，则

\[\begin{aligned} (1+x)^n&=(1+x)^{k_1p}\cdot(1+x)^{b_1}\\ &\equiv (1+x^p)^{k_1}\cdot (1+x)^{b_1}\\ &\equiv \big(\sum_{i=0}^{k_1}\binom{k_1}{i}x^{ip}\big)\big(\sum_{j=0}^{b_1}\binom{b_1}{j}x^{j}\big)\pmod p \end{aligned} \]

取其中的 \(x^m\) 项，则有

\[\binom{n}{m}\cdot x^m\equiv \binom{k_1}{k_2}\cdot x^{k_2p}\cdot\binom{b_1}{b_2}\cdot x^{b_2}\pmod p \]

由于 \(p\) 是质数，两边同除以 \(x^m\)，可得

\[\binom{n}{m}\equiv \binom{k_1}{k_2}\cdot\binom{b_1}{b_2}\equiv \binom{n\bmod p}{m\bmod p}\times\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\pmod{p} \]

模板P3807

点击查看代码

#include<iostream>
#include<cstdio>
typedef long long ll;
const int N=1e5+13;
int mul[N],invmul[N];
inline int qpow(int a,int k,int p){int s=1;for(;k;k>>=1,a=(ll)a*a%p)if(k&1)s=(ll)s*a%p;return s;}
inline void init(int mod){
	mul[0]=invmul[0]=1;int n=mod-1;
	for(int i=1;i<=n;++i) mul[i]=(ll)mul[i-1]*i%mod;
	invmul[n]=qpow(mul[n],mod-2,mod);
	for(int i=n-1;i;--i) invmul[i]=(ll)invmul[i+1]*(i+1)%mod;
}
inline int C(int n,int m,int p){if(n<m)return 0;return (ll)mul[n]*invmul[m]%p*invmul[n-m]%p;}
int Lucas(int n,int m,int p){
	if(!m) return 1;
	return (ll)Lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int main(){init(100000);int T;scanf("%d",&T);while(T--){
	int n,m,p;scanf("%d%d%d",&n,&m,&p);
	init(p);
	printf("%d\n",Lucas(n+m,m,p));
}
	return 0;
}

十三、扩展 Lucas

要计算

\[\binom{n}{m}\bmod P \]

其中不保证 \(P\) 是质数。

首先想到可以将 \(P\) 质因数分解，即

\[P=\prod p_i^{k_i} \]

那么如果能把每个 \(\binom{n}{m}\bmod p_i^{k_i}\) 都求出来，就可以使用 CRT 合并了。

现在把问题转化为求

\[\binom{n}{m}\bmod P^k(P\in prime) \]

这个东西是

\[\frac{n!}{m!(n-m)!}\bmod P^k \]

但是因为 \(m!,(n-m)!\) 不一定和 \(P^k\) 互质，所以不一定有逆元。考虑把式子化成

\[\frac{\frac{n!}{P^x}}{\frac{m!}{P^y}\cdot\frac{(n-m)!}{P^z}}P^{x-y-z}\bmod P^k \]

其中 \(\frac{n!}{P^x},\frac{m!}{P^y},\frac{(n-m)!}{P^z}\) 都与 \(P^k\) 互质。

现在把问题转化为求

\[\begin{aligned} \frac{n!}{P^x}\bmod P^k&=(P\cdot 2P\cdot 3P\cdot\ldots\cdot \lfloor\frac{n}{P}\rfloor P)(1\cdot 2\cdot 3\cdot\ldots\cdot n)\\ &=P^{\lfloor\frac{n}{P}\rfloor}(\lfloor\frac{n}{P}\rfloor)!\prod_{i=1,i\not\equiv0\pmod P}^n i\\ &\equiv P^{\lfloor\frac{n}{P}\rfloor}(\lfloor\frac{n}{P}\rfloor)!\big(\prod_{i=1,i\not\equiv 0\pmod P}^{P^k}i\big)^{\lfloor\frac{n}{P^k}\rfloor}\prod_{i=\lfloor\frac{n}{P^k}\rfloor P^k,i\not\equiv 0\pmod P}^n i\pmod {P^k} \end{aligned} \]

设 \(f(n)=\frac{n!}{P^x}\bmod P^k\)，则有递推式

\[f(n)=f(\lfloor\frac{n}{P}\rfloor)\big(\prod_{i=1,i\not\equiv 0\pmod P}^{P^k}i\big)^{\lfloor\frac{n}{P^k}\rfloor}\prod_{i=\lfloor\frac{n}{P^k}\rfloor P^k,i\not\equiv 0\pmod P}^n \]

这一步计算 \(f(n)\) 的复杂度就是 \(O(P\log_P n)\)。

考虑怎么计算 \(P^{x-y-z}\)。设 \(g(n)=x\)，其中 \(x\) 就是上边 \(\frac{n!}{P^x}\) 的那个。通过上面那个 \(f\) 的递推式可以得到

\[g(n)=\lfloor\frac{n}{P}\rfloor+g(\lfloor\frac{n}{P}\rfloor) \]

这一步的复杂度是 \(O(\log_P n)\)。

所以总的式子就是

\[\begin{aligned} \binom{n}{m}&\equiv\frac{n!}{m!(n-m)!}\\ &\equiv\frac{\frac{n!}{P^x}}{\frac{m!}{P^y}\cdot\frac{(n-m)!}{P^z}}P^{x-y-z}\\ &\equiv \frac{f(n)}{f(m)f(n-m)}P^{g(n)-g(m)-g(n-m)} \pmod {P^k}(P\in prime) \end{aligned} \]

然后拿 CRT 合并就行了。总的复杂度就是 \(O(P\log_P n+\log^2 P)\)，可以通过。

模板P4720

点击查看代码

#include<iostream>
#include<cstdio>
typedef long long ll;
const int N=1e6+13,M=30+13;
int prm[N],pcnt,a[M],b[M];
bool vis[N];
inline void init(int lim){
	for(int i=2;i<=lim;++i){
		if(!vis[i]) prm[++pcnt]=i;
		for(int j=1;j<=pcnt&&i*prm[j]<=lim;++j){
			vis[i*prm[j]]=1;
			if(i%prm[j]==0) break;
		}
	}
}
inline int qpow(int a,ll k,int p){int s=1;for(;k;k>>=1,a=(ll)a*a%p)if(k&1)s=(ll)s*a%p;return s;}
inline void exgcd(int a,int b,int &x,int &y){!b?(x=1,y=0):(exgcd(b,a%b,y,x),y-=(a/b)*x);}
inline int inv(int a,int p){int x,y;exgcd(a,p,x,y);return (x%p+p)%p;}
int f(ll n,int p,int pk){
	if(!n) return 1;
	int res1=1,res2=1;
	for(int i=1;i<=pk;++i)
		if(i%p) res1=(ll)res1*i%pk;
	res1=qpow(res1,n/pk,pk);
	for(ll i=pk*(n/pk);i<=n;++i)
		if(i%p) res2=(ll)res2*(i%pk)%pk;
	return (ll)f(n/p,p,pk)*res1%pk*res2%pk;
}
ll g(ll n,int p){return n<p?0:(n/p)+g(n/p,p);}
inline int solve(ll n,ll m,int p,int pk){
	return (ll)f(n,p,pk)*inv(f(m,p,pk),pk)%pk*inv(f(n-m,p,pk),pk)%pk*qpow(p,g(n,p)-g(m,p)-g(n-m,p),pk)%pk;
}
using namespace std;
inline int exLucas(ll n,ll m,int p){
	int tmp=p,cnt=0;
	for(int i=1;i<=pcnt&&prm[i]*prm[i]<=p;++i){
		if(p%prm[i]) continue;
		int pk=1;
		while(!(tmp%prm[i])) pk*=prm[i],tmp/=prm[i];
		a[++cnt]=pk,b[cnt]=solve(n,m,prm[i],pk);
	}
	if(tmp!=1) a[++cnt]=tmp,b[cnt]=solve(n,m,tmp,tmp);
	int ans=0;
	for(int i=1;i<=cnt;++i){
		int M=p/a[i],ni=inv(M,a[i]);
		ans+=(ll)M*ni%p*b[i]%p,ans%=p;
	}
	return ans;
}
int main(){
	init(1000000);
	ll n,m;int p;
	scanf("%lld%lld%d",&n,&m,&p);
	printf("%d\n",exLucas(n,m,p));
	return 0;
}

十四、狄利克雷前缀和

有一个序列 \(a_{1\ldots n}\)，要求一个序列 \(b_{1\ldots n}\) 满足 \(b_k=\sum\limits_{i|k}a_i\)。

朴素的做法就是直接枚举每个数向其倍数贡献。狄利克雷前缀和说的是，考虑把每个数分解质因子之后，每个质因子都相当于是求了一个高维前缀和，所以可以直接枚举质因子求前缀和。

点击查看代码

for(int i=1;i<=pcnt;++i)
	for(int j=1;j<=n&&j*prm[i]<=n;++j)
		a[j*prm[i]]+=a[j];

十五、莫比乌斯反演

前置知识：数论分块，狄利克雷卷积。不一一赘述。

1. 莫比乌斯函数

定义

\[\mu(n)=\left\{ \begin{aligned} 1\ \ \ &n=1\\ 0\ \ \ &n 含有平方因子\\ (-1)^k\ \ \ &k 为 n 本质不同质因子个数 \end{aligned} \right. \]

性质

\[\sum_{d|n}\mu(d)=\left\{ \begin{aligned} 1\ \ \ n=1\\ 0\ \ \ n\not=1 \end{aligned} \right. \]

证明：

令 \(n=\prod_{i=1}^k p_i^{c_i},n'=\prod_{i=1}^k p_i\)

则根据二项式定理可得 \(\sum_{d|n}\mu(d)=\sum_{d|n'}\mu(d)=\sum_{i=0}^k C_k^i\cdot(-1)^i=(1+(-1))^k\)

那么 \(\sum_{d|n}\mu(d)=1\) 当且仅当 \(k=0\)，也就是 \(n=1\) 时取到。

反演常用结论：\([\gcd(i,j)=1]=\sum_{d|\gcd(i,j)}\mu(d)\)

扩展：

\[\varphi * \mathrm{1}=\mathrm{id} \]

证明：

令 \(n=\prod_{i=1}^k p_i^{c_i}\)。由于 \(\varphi\) 是积性函数，所以只需要证明 \(n'=p^c\) 时 \(\varphi * \mathrm{1}=\sum_{d|n'}\varphi(d)=\mathrm{id}\) 即可。

那么

\[\begin{aligned} \varphi*\mathrm{1}(n')&=\sum_{d|n'}\varphi(d)\\ &=\sum_{i=0}^c \varphi(p^i)\\ &=1+p^0\cdot (p-1)+p^1\cdot (p-1)+\ldots+p^{c-1}\cdot (p-1)\\ &=1+\frac{p^c-1}{p-1}\cdot (p-1)\\ &=p^c\\ &=\mathrm{id}(n') \end{aligned} \]

并且，该式子两边同时卷 \(\mu\) 可以得到 \(\varphi(n)=\sum_{d|n}d\cdot\mu(\frac{n}{d})\)，即 \(\varphi=\mathrm{id}* \mu\)。

2. 莫比乌斯变换（反演）

若 \(f(n)=\sum_{d|n}g(d)\)，那么有 \(g(n)=\sum_{d|n}\mu(d)f(\frac{n}{d})\)。

这个可以用 \(\mu * \mathrm{1}=\epsilon\) 来解释。把第一个式子两边卷上 \(\mu\) 就可以得到第二个式子了。注意任何函数卷上 \(\epsilon\) 都是自己。