lp5431 【模板】乘法逆元2


挺套路的。
先计算出所有数的积的逆元,再计算除了这个数以外的数的积,然后乘一起,这样就完成了线性求逆元。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
using namespace std;

typedef long long ll;
const int N=5000005;
inline int rd(){
	int rt=0;char ch=getchar();
	while(!isdigit(ch)){ch=getchar();}
	while(isdigit(ch)){rt=rt*10+ch-'0';ch=getchar();}
	return rt;
}
int n,MOD,k;
int a[N],pr[N],sf[N];
inline int pw(int A,int X){
	int RT=1;
	while(X){
		if(X&1){
			RT=1ll*RT*A%MOD;
		}
		A=1ll*A*A%MOD;X>>=1;
	}
	return RT;
}

void init(){
	scanf("%d%d%d",&n,&MOD,&k);
	int pl=1;
	pr[0]=sf[n+1]=1;
	for(int i=1;i<=n;++i){
		a[i]=rd();
	}
	for(int i=1;i<=n;++i){
		pr[i]=1ll*pr[i-1]*a[i]%MOD;
	}
	pl=pw(pr[n],MOD-2);
	for(int i=n;i>=1;--i){
		sf[i]=1ll*sf[i+1]*a[i]%MOD;
	}
	int nk=1;
	int ans=0;
	for(int i=1;i<=n;++i){
		nk=1ll*nk*k%MOD;
		ans+=1ll*(1ll*pr[i-1]*sf[i+1]%MOD)*(1ll*pl*nk%MOD)%MOD;ans%=MOD;
	}
	printf("%d\n",ans);
	
}
int main(){
	init();
	return 0;
}

lp4980 【模板】Polya定理

Polya定理是一个关于置换群中组合计数的定理。
首先我们来了解Burnside引理。这个引理的证明较为复杂,在这里仅介绍其内容。
Burnside引理的内容是这样的:在置换意义下本质不同的染色方案数,等于单个置换不会改变的染色方案数的平均值。
那么,单个置换不会改变的染色方案数要怎么求呢?我们不妨把置换关系建成一张图,那么这张图必然由若干个环构成。也就是说,对于同一种置换,可能存在若干个环,而每个初始状态就是环上的一个节点。Polya定理描述的就是,某一种置换不会改变的方案数,等于颜色数的「这个置换对应的循环节个数」

阐述了上述两个定理,现在我们来看这一题。
我们观察到,这个环上存在 n 种置换。对于置换 i 而言,它的循环大小是:
$$
\frac{n}{\gcd(n,i)}
$$
这也就意味着,它的循环个数是:
$$
\gcd(n,i)
$$
那么我们要求的答案就是:
$$
\frac{\sum_{i=1}^nn^{(n,i)}}{n}
$$
然而这样计算答案的复杂度是 O(n) 的,我们无法接受。
我们不妨考虑枚举这个公因子。那么,容易证明的是,每一个公因子d对应的置换个数是\(\varphi(\frac{n}{d})\)

这是因为, d 是 \(\frac{n}{d}\) 个数的因数,而这些数中,如果它和\(\frac{n}{d}\)存在公因子的话,那么它和n的最大公因数必然大于d。

所以,我们求的就是:
$$
\frac{\sum_{d|n}\varphi(\frac{n}{d})n^{\frac{n}{d}}}{n}=\sum_{d|n}\varphi(\frac{n}{d})n^{\frac{n}{d}-1}
$$
我们只要枚举n的因数即可。

然而求这里的 \(\varphi\) 倒可能存在一些问题。事实上,由于我们要取的因数的值域是 \(10^9\) ,我们无法使用欧拉筛来预处理,而需要现场做。

这样的复杂度是不是对的呢?它的最坏情况是 \(O(Tn^{\frac{3}{4}})\) 的,但由于我并不会分析的一些性质,它不会达到这个值。因此可以通过此题。

另:括号忘记加,爆O泪两行。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
using namespace std;
typedef long long ll;
const ll MOD=1000000007;
const int N=1000005;
inline ll pw(ll A,ll X){
	ll RT=1;
	while(X){
		if(X&1){
			RT*=A;RT%=MOD;
		}
		A*=A;A%=MOD,X>>=1;
	}
	return RT;
}
int n;
inline int phi(int X){
	int RT=1;
	for(int i=2;i*i<=X;++i){
		if(X%i==0){
			X/=i;RT*=i-1;
			while(X%i==0){
				X/=i;RT*=i;
			}
		}
		if(X==1){
			break;
		} 
	}
	if(X>1){
		RT*=X-1;
	}
	return RT;
}
void init(){
	scanf("%d",&n);
	ll ans=0;
	for(int i=1;i*i<=n;++i){
		if(n%i==0){
			ans+=(1ll*phi(n/i)*pw(n,i-1))%MOD;
			if(i*i!=n){
				ans+=(1ll*phi(i)*pw(n,n/i-1))%MOD;
			}
			ans%=MOD;
//			printf("%d %lld\n",i,ans);
		}
	}
	printf("%lld\n",ans);
}
int main(){
	int T;
	scanf("%d",&T);
	while(T--){
		init();
	}
	return 0;
}

lp5395 【模板】第二类斯特林数·行

第二类斯特林数,指的是一组表示「将n个不同的元素划分为m个非空不相交集的方案数」的组合数。有时写作\(S(n,m)\)或\(\left\{\begin{matrix}n\\m\end{matrix}\right\}\)
从题意来看,我们可以容易地想到一个递推方案:每个新的元素,可以为它新开一个集合,或者放到已有的任何一个集合里面。所以我们得到一个递推式:
$$S(n,m)=S(n-1,m-1)+m*S(n-1,m)$$
然后问题是怎么求值。

首先我们有一个式子,我们称之为二项式反演的通项式。
$$ f(n)=\sum_{i=0}^nC_{n}^{i}g(n) \Leftrightarrow g(n)=\sum_{i=0}^n(-1)^{n-i}C_{i}{n}f(n)$$
举个例子:
$$\begin{matrix} f_{1}&=&g_{1}\\f_{2}&=&g_{1}&+&g_{2}\\f_{3}&=&g_{1}&+&2*g_{2}&+&g_{3}\\f_{4}&=&g_{1}&+&3*g_{2}&+&3*g_{3}&+&g_{4} \end{matrix}$$

由此可得:

$$\begin{matrix} g_{1}&=&f_{1}\\g_{2}&=&f_{2}&-&f_{1}\\g_{3}&=&f_{3}&-&2*g_{2}&-&g_{1}\\&=&f_{3}&-&2*f_{2}&+&f_{1}\\g_{4}&=&f_{4}&-&3*g_{3}&-&3*g_{2}&-&g_{1}\\&=&f_{4}&-&3*f_{3}&-&3*g_{2}&-&4*g_{1}\\&=&f_{4}&-&3*f_{3}&+&3*f_{2}&-&f_{1} \end{matrix}$$

然后我们考虑\(m^n\)的组合意义,也就是,将\(n\)个不同的球,放到\(m\)个不同的盒子(可以有空盒)的方案数。
我们不妨枚举有装东西的盒子的个数。令它为\(i\),那么选取这些盒子的方案数即是\(C_{m}^{i}\)。
选取了盒子之后,问题就转化为了将\(n\)个不同物品装到\(m\)个不同盒子里,使得每个盒子非空。这就相当于,将\(n\)个不同物品装到\(m\)个相同盒子里,使得每个盒子非空的方式——也就是\(S(n,m)\),乘上盒子的排列顺序——也就是\(m!\)。
形式化的说,就是:
$$m^n=\sum_{i=0}^{m}S(n,i)i!C_{m}^{i}$$
我们发现这个式子的形式符合二项式反演的通项式。
因而我们将\(m^n\)作为\(f\),将\(S(n,i)i!\)作为\(g\)。
那么反演得到:
$$S(n,m)=\frac{1}{m!}\sum_{i=0}^{m}(-1)^{m-i}i^nC_{m}^{i}$$
考虑\(C_{m}^{i}=\frac{m!}{i!(m-i)!}\),我们可以将上式化简得到:
$$S(n,m)=\sum_{i=0}^{m}\frac{(-1)^{m-i}i^n}{i!(m-i)!}$$
于是我们有:
$$\left\{\begin{matrix}n\\m\end{matrix}\right\}=\sum_{i=0}^{m}\frac{(-1)^{m-i}i^n}{i!(m-i)!}$$
考虑函数卷积的本质,对于卷积\(f=gu\),其本质是: $$f(n)=\sum_{i=0}^{n}g(i)u(n-i)$$ 因而,我们发现,上式本质上就是: $$\left\{\begin{matrix}n\\m\end{matrix}\right\}=\sum_{i=0}^{m}\frac{(-1)^{m-i}}{(m-i)!}*\frac{i^n}{i!}$$
上一个NTT卷积一下就好了。

#include<iostream>
#include<cstdio>

#define Swap(A,B) (A^=B^=A^=B)
const long long P=167772161;
const long long g0=3,gi=55924054;
int L=1,R[1<<21|1];
long long invL;
inline int pw(int A,int X){
	long long BS=A,RT=1;
	while(X){
		if(X&1){
			RT=RT*BS%P;
		}
		BS=BS*BS%P;
		X>>=1;
	}
	return RT;
}
inline void prpr(int LEN){
	int B=0;
	while(L<=LEN){
		L<<=1;
		++B;
	}
	invL=P-(P-1)/L;
	for(int i=0;i<L;++i){
		R[i]=R[i>>1]>>1|(i&1)<<(B-1);
	}
}
inline void FNTT(int *A,int typ){
	for(int i=0;i<L;++i){
		if(R[i]<i){
			Swap(A[R[i]],A[i]);
		}
	}
	int gn,g,X,Y,M;
	for(int i=2;i<=L;i<<=1){
		M=i>>1;
		gn=pw(~typ?g0:gi,(P-1)/i);
		for(int j=0;j<L;j+=i){
			g=1;
			for(int k=0;k<M;++k,g=1ll*g*gn%P){
				X=A[j+k],Y=1ll*g*A[M+j+k]%P;
				A[j+k]=(X+Y)%P,A[M+j+k]=(X-Y)%P;
			}
		}
	}
}
int n,a[1<<21|1],b[1<<21|1],inv[1<<21|1];
void init(){
	scanf("%d",&n);
	prpr(n+1<<1);
	inv[0]=inv[1]=1;
	for(int i=2;i<=n;++i){
		inv[i]=1ll*(P-P/i)*inv[P%i]%P;
	}
	for(int i=1;i<=n;++i){
		inv[i]=1ll*inv[i-1]*inv[i]%P;
	}
	for(int i=0;i<=n;++i){
		a[i]=1ll*pw(-1,i)*inv[i]%P;
		b[i]=1ll*pw(i,n)*inv[i]%P;
	}
	FNTT(a,1);
	FNTT(b,1);
	for(int i=0;i<L;++i){
		a[i]=1ll*a[i]*b[i]%P;
	}
	FNTT(a,-1);
	for(int i=0;i<=n;++i){
		printf("%d ",1ll*(a[i]+P)*invL%P);
	}
}

int main(){
	init();
	return 0;
}

lp3455 POI2007 ZAP-Queries

曾经有一道题,叫做YY的GCD,它求的是这样一个值:
$$\begin{equation}\begin{split}\sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{p\in P}[gcd(x,y)==p] \end{split}\end{equation}$$
然后观察这道题的求和式:
$$\begin{equation}\begin{split}\sum_{x=1}^{n}\sum_{y=1}^{m}[gcd(x,y)==d] \end{split}\end{equation}$$
长得很像对不对?
事实上就是前者的简化版。
剩下的操作就很套路了。
$$\begin{equation}\begin{split}&\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(x,y)==1]\\=&\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{k|gcd(x,y)}\mu(k)\\=&\sum_{x=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{dk}\rfloor}\mu(k)\\=&\sum_{k=1}^{min(n,m)}\mu(k)\lfloor\frac{n}{dk}\rfloor\lfloor\frac{m}{dk}\rfloor \end{split}\end{equation}$$
然后预处理出莫比乌斯函数跑数论分块即可。

#include<iostream>
#include<cstdio>

const int N = 100000+5;

int p[N/10],mu[N],n,m,d;
bool ip[N];
void prpr(){
	p[0]=0;mu[1]=1;ip[1]=1;
	for(int i=2;i<=100000;++i){
		if(!ip[i]){
			p[++p[0]]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=p[0]&&p[j]*i<=100000;++j){
			ip[i*p[j]]=1;
			if(!(i%p[j])){
				mu[i*p[j]]=0;
				break;
			}else{
				mu[i*p[j]]=-mu[i];
			}
		}
	}
	for(int i=2;i<=100000;++i){
		mu[i]+=mu[i-1]; 
	}
}
long long ans;
void init(){
	scanf("%d%d%d",&n,&m,&d);
	n>m?n^=m^=n^=m:0;
	ans=0;
	int k=0;
	for(int i=1;i<=n;i=k+1){
		k=std::min(n/(n/i),m/(m/i));
		ans+=1ll*(n/(i*d))*(m/(i*d))*(mu[k]-mu[i-1]);
	}
	printf("%lld\n",ans);
}

int main(){
	prpr();
	int T;
	scanf("%d",&T);
	while(T--){
		init();
	}
	return 0;
}

lp4245 【模板】任意模数NTT

我们首先考虑一个很暴力的玩法,直接找一个很大很大很大的模数,然后用int128上一个朴素NTT,再将结果对题目给的模数取模,似乎就可以了。
然而仔细考虑一下发现这个做法实在太暴力了(才不是因为找不到合适的模数呢。)我们考虑另一种方案。
现在假设我们已经用好几种模数求出最终的各种值了,那么原值是多少呢?很显然可以用中国剩余定理来求出。
我们计算了一下,发现只要用三个质数,值域就足以用中国剩余定理求出最终值了。
然而我们知道,模运算没有交换律。故而我们不能用常规的中国剩余定理合并,应该需要使用一种类似于扩展中国剩余定理的递推合并方法来合并它们的解。
这就完成了任意模数的NTT。
另外要注意数组大小…

#include<iostream>
#include<cstdio>

inline void Swap(int &A,int &B){
	A^=B^=A^=B;
}

inline long long mlt(long long A,long long X,long long P){
	long long BS=A,RT=0;
	while(X){
		if(X&1){
			RT+=BS;
			RT%=P;
		}
		BS+=BS;
		BS%=P;
		X>>=1;
	}
	return RT;
}

inline long long pw(long long A,long long X,long long P){
	long long BS=A,RT=1;
	while(X){
		if(X&1){
			RT*=BS;
			RT%=P;
		}
		BS*=BS;
		BS%=P;
		X>>=1;
	}
	return RT;
}

const long long MOD[3]={469762049,998244353,1004535809};
const long long MM=468937312667959297;
const long long g0=3,gi[3]={156587350,332748118,334845270};

int L=1,inv[3],R[1<<21|1];
inline void prpr(int LEN){
	int B=0;
	while(L<=LEN){
		L<<=1;
		++B;
	}
	for(int i=0;i<3;++i){
		inv[i]=MOD[i]-(MOD[i]-1)/L;
	}
	for(int i=0;i<L;++i){
		R[i]=R[i>>1]>>1|(i&1)<<(B-1);
	}
}

inline void NTT(int *A,int typ,int P){
	for(int i=0;i<L;++i){
		if(R[i]<i){
			Swap(A[i],A[R[i]]);
		}
	}
	int bs,nw,X,Y,M;
	for(int i=2;i<=L;i<<=1){
		M=i>>1;
		bs=pw(~typ?g0:gi[P],(MOD[P]-1)/i,MOD[P]);
		for(int j=0;j<L;j+=i){
			nw=1;
			for(int k=0;k<M;++k,nw=1ll*nw*bs%MOD[P]){
				X=A[j+k],Y=1ll*nw*A[j+k+M]%MOD[P];
				A[j+k]=(X+Y)%MOD[P];
				A[j+k+M]=(X-Y+MOD[P])%MOD[P];
			}
		}
	}
}
int C[1<<21|1],D[1<<21|1],ans[3][1<<21|1];
inline void FNTT(int *A,int *B,int P){
	for(int i=0;i<L;++i){
		C[i]=A[i],D[i]=B[i];
	}
	NTT(C,1,P);
	NTT(D,1,P);
	for(int i=0;i<L;++i){
		C[i]=1ll*C[i]*D[i]%MOD[P];
	}
	NTT(C,-1,P);
	for(int i=0;i<L;++i){
		ans[P][i]=(int)(1ll*(C[i]+MOD[P])*inv[P]%MOD[P]);
	}
}
int a[1<<21|1],b[1<<21|1],n,m,p;
long long t[3];
void init(){
	scanf("%d%d%d",&n,&m,&p);
	for(int i=0;i<=n;++i){
		scanf("%d",a+i);
	}
	for(int i=0;i<=m;++i){
		scanf("%d",b+i);
	}
	prpr(n+m);
	FNTT(a,b,0);FNTT(a,b,1);FNTT(a,b,2);
	t[0]=pw(MOD[1]%MOD[0],MOD[0]-2,MOD[0]);t[1]=pw(MOD[0]%MOD[1],MOD[1]-2,MOD[1]);t[2]=pw(MM%MOD[2],MOD[2]-2,MOD[2]);//分别求出要用到的三个逆元。 
	long long T1,T2;
	for(int i=0;i<=n+m;++i){
		T1=(mlt(1ll*ans[0][i]*MOD[1]%MM,t[0],MM)+mlt(1ll*ans[1][i]*MOD[0]%MM,t[1],MM))%MM;//直接用CRT合并一二式。 
		T2=((ans[2][i]-T1)%MOD[2]+MOD[2])%MOD[2]*t[2]%MOD[2];//用EXCRT将第三式合并之。 
		printf("%d ",(int)(((T2%p)*(MM%p)%p+T1%p)%p));//求得最终解。 
	}
}

int main(){
	init();
	return 0;
}

lp4777 【模板】扩展中国剩余定理(EXCRT)

$$ \left\{ \begin{matrix}x\equiv a_{1} \pmod{p_1} \\x\equiv a_{2}\pmod{p_2}\\\cdots\\x\equiv a_{n}\pmod{p_i}\end{matrix} \right.$$
对于这个方程,当所有模数两两互质的时候,可以使用中国剩余定理来求解。
然而,倘若并不满足这个条件,普通的中国剩余定理就未必能够奏效了。
我们尝试将中国剩余定理的成立条件进行推广,从而得到一种被称为「扩展中国剩余定理」的算法,以解决模数并不保证互质的情况下的线性同余方程组求解。
对于每个方程的解,它构成了一个解集。而我们需要找到的是这个解集中最小的且满足下一个方程的正整数解。
容易证明,一些线性同余方程组的解集一定是关于所有模数的最小公约数同余的一个数。
于是,我们令\(x_k\)为方程\(1\to k\)的最小正整数解:
$$\begin{equation}\begin{split}P,st:P&=&lcm({p_{i},i\in [1,k]})\\x_{k+1}&=&x_{k}+tP\\x_{k+1}&\equiv&a_{k+1}&\pmod{p_{k+1}}\\x_{k}+tP&\equiv&a_{k+1}&\pmod{p_{k+1}}\\tP&\equiv&a_{k+1}-x_{k}&\pmod{p_{k+1}}\end{split}\end{equation}$$
最后的那个式子可以转化形式然后上扩欧。
这就完成了一种递推求解不保证模数互质的线性同余方程组的方案。

#include<iostream>
#include<cstdio>

long long MOD=1;


inline long long mlt(long long A,long long X,long long P){
	long long BS=A,RT=0;
	while(X){
		if(X&1){
			RT+=BS;
			RT%=P;
		}
		BS+=BS;
		BS%=P;
		X>>=1;
	}
	return RT;
}
inline long long gcd(long long A,long long B){
	return B?gcd(B,A%B):A;
}
inline long long exgcd(long long A,long long B,long long &X,long long &Y,long long D=0){
	return B?(D=exgcd(B,A%B,Y,X),Y-=A/B*X,D):(X=1,Y=0,A);
}
long long n,a[100005],b[100005];
void init(){
	scanf("%lld",&n);
	for(int i=1;i<=n;++i){
		scanf("%lld%lld",&b[i],&a[i]);
	}
	long long X,Y,T,C,ans=a[1];
	MOD=b[1];
	for(int i=2;i<=n;++i){
		C=(a[i]-ans%b[i]+b[i])%b[i];
		X=exgcd(MOD,b[i],T,Y);
		if(C%X!=0){
			puts("F**KING HARD!");
			return;
		}
		T=mlt(T,C/X,b[i]/X);
		ans+=T*MOD;
		MOD/=X;
		MOD*=b[i];
		ans=(ans%MOD+MOD)%MOD;
		
	}
	printf("%lld\n",(ans+MOD)%MOD);
}

int main(){
	init();
	return 0;
}

lp3868 TJOI2009 猜数字

这是一道中国剩余定理(又名孙子定理、CRT)的模板题。
它最早来自于「孙子算经」,形式化的说,它是一种用于求解如下的线性同余方程组的解的算法。
$$ \left\{ \begin{matrix}x\equiv a_{1} \pmod{p_1} \\x\equiv a_{2}\pmod{p_2}\\\cdots\\x\equiv a_{n}\pmod{p_i}\end{matrix} \right.$$
我们令:
$$\begin{equation}\begin{split}P&=&\prod_{i=1}^{n}p_{i}\\t_{i},st: \frac{P}{p_{i}}t_{i}&\equiv &1\pmod{p_i}\\\exists x,st: x&\equiv &\sum_{i=1}^n \frac{P}{p_{i}}a_{i}t_{i}\\x&\in&\{x_0+kP,k\in Z\pmod{P}\}\end{split}\end{equation}$$
然而这是怎么证明的呢?我们有如下推导:
$$\begin{equation}\begin{split}\forall i,\forall j,st:i\neq j, a_{i}t_{i}\frac{P}{p_{i}}&\equiv&0&\pmod{p_i}\\a_{i}t_{i}\frac{P}{p_i}&\equiv&a_{i} &\pmod{p_i}\\\forall j,\sum_{i=1}^{n}a_{i}t_{i}\frac{P}{p_{i}}&\equiv&a_{j}t_{j}\frac{P}{p_{j}} \equiv a_{j}&\pmod{p_{j}}\\\exists x_0&\equiv& \sum_{i=1}^{n}a_{i}t_{i}\frac{P}{p_{i}}&\pmod{p_{j}}\\\end{split}\end{equation}$$
故而我们就可以求得上述方程的最小正整数解。
另外,看到这题的数据范围就可以猜到又是龟速乘登场的时刻了。
这一题还有个坑点就是读入的数可能是负数…

#include<iostream>
#include<cstdio>

long long MOD=1;


inline long long mlt(long long A,long long X){
	long long BS=A,RT=0;
	while(X){
		if(X&1){
			RT+=BS;
			RT%=MOD;
		}
		BS+=BS;
		BS%=MOD;
		X>>=1;
	}
	return RT;
}
inline long long exgcd(long long A,long long B,long long &X,long long &Y,long long D=0){
	return B?(D=exgcd(B,A%B,Y,X),Y-=A/B*X,D):(X=1,Y=0,A);
}
inline long long pw(long long A,long long X){
	long long BS=A,RT=1;
	while(X){
		if(X&1){
			RT=mlt(RT,BS)%MOD;
		}
		BS=mlt(BS,BS)%MOD;
		X>>=1;
	}
	return RT;
}
long long n,a[100005],b[100005];
void init(){
	scanf("%lld",&n);
	for(int i=1;i<=n;++i){
		scanf("%lld",&a[i]);
	}
	for(int i=1;i<=n;++i){
		scanf("%lld",&b[i]);
		a[i]=(a[i]%b[i]+b[i])%b[i];//注意负数!! 
		MOD*=b[i];
	}
	long long X,Y,T,ans=0;
	for(int i=1;i<=n;++i){
		T=MOD/b[i];
        exgcd(T,b[i],X,Y);
        X=(X%b[i]+b[i])%b[i];;
        ans=(ans+mlt(mlt(T,X),a[i]))%MOD;
	}
	printf("%lld\n",(ans+MOD)%MOD);
}

int main(){
	init();
	return 0;
}

lp2044 NOI2012 随机数生成器

我们不妨将原式递归展开,可以得到形如此的式子:
$$f_{n}=a^nf_{0}+\frac{(a^{n}-1)c}{a-1}$$
对于这个式子的值,我们看似可以上一个快速幂直接求得。
然而仔细观察上式的情况,发现如果出现了奇妙的模数时,\(a-1\)是不一定有逆元的!
这就很令人难受了。
我们找到原来的问题,我们不得不求下述式子的值:
$$f_{n}=a^nf_{0}+\sum_{i=0}^{n-1}a^ic=a^nf_{0}+(\sum_{i=0}^{n-1}a^i)c$$
求值的最大障碍是这个式子:
$$\sum_{i=0}^{n-1}a^i$$
将这个式子分奇偶讨论后递归求解。
我们令:
$$u_{x}=\sum_{i=0}^{x-1}a^i$$
若是偶数,则有:
$$u_x=\sum_{i=0}^{\frac{x-1}{2}}a^i+a^{\frac{x+1}{2}} (\sum_{i=0}^{\frac{x-1}{2}}a^i)$$
即,
$$u_x=(a^{\frac{x+1}{2}}+1)u_{\frac{x-1}{2}}$$
若是奇数,则只需将第一项单独计算。
这样就可以在对数复杂度内对式子求解了。
另外,考虑到模数的数据范围,我们需要使用一种被称为「龟速乘」,也就是「快速幂」在乘法意义上的等同的操作,使得不会出现乘法,也就不会爆范围。

#include<iostream>
#include<cstdio>

long long MOD,a,c,x0,n,g;
//必须注意,不能重载一个全部是标准型的运算符。 
inline long long mlt(long long A,long long X){
	long long BS=A,RT=0,CNT=X;
	while(CNT){
		if(CNT&1){
			RT+=BS;
			RT%=MOD;
		}
		CNT>>=1;
		BS+=BS;
		BS%=MOD;
	}
	return RT;
}
inline long long pw(long long A,long long X){
	long long BS=A,RT=1,CNT=X;
	while(CNT){
		if(CNT&1){
			RT=mlt(RT,BS)%MOD;
		}
		CNT>>=1;
		BS=mlt(BS,BS)%MOD;
	}
	return RT;
}
inline long long calc(long long X){
	if(X==2){
		return (a+1)%MOD;
	}
	if(X==1){
		return 1;
	}
	if(X&1){
		return (pw(a,X-1)+calc(X-1))%MOD;
	}
	return mlt((pw(a,(X+1)>>1)+1),calc(X>>1));
}
void init(){
	scanf("%lld%lld%lld%lld%lld%lld",&MOD,&a,&c,&x0,&n,&g);
	printf("%lld",((mlt(pw(a,n),x0)+mlt(calc(n),c))%MOD)%g);
}

int main(){
	init();
	return 0;
}

lp2257 YY的GCD

不妨设\(n<m\)
首先将原提问形式化,可以得到原式为:
$$\sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{p\in P}[gcd(x,y)==p] $$
将所有的\(x,y\)都除以\(gcd(x,y)\),则原式可转化为:
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}[gcd(x,y)==1]$$
然后,根据莫比乌斯函数的基本性质,我们有:
$$\sum_{d|n}\mu(d)=[n==1]$$
替换入原式可以得到:
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}\sum_{d|gcd(x,y)}\mu(d)$$
将\(d|gcd(x,y)\)一式化简可得:
$$d|gcd(x,y)=d|x\&\&d|y$$
代入原式:
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}\sum_{d|x\&\&d|y}\mu(d)$$
将所有的\(x,y\)都除以\(d\),则原式可转化为:
$$\sum_{x=1}^{\lfloor \frac{n}{pd} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{pd} \rfloor}\sum_{p\in P}\mu(d)$$
由于\(x,y\)对后半部分式子无影响,且\(pd\le n\)原式可变式为:
$$\sum_{p\in P}\sum_{d=1}^{\lfloor \frac{n}{p} \rfloor}\mu(d)\lfloor \frac{n}{pd} \rfloor \lfloor \frac{m}{pd} \rfloor$$
为了减少枚举的数量,不妨设\(q=pd\),那么我们又可以将原式转化为:
$$\sum_{q=1}^{n}\lfloor \frac{n}{q} \rfloor \lfloor \frac{m}{q} \rfloor \sum_{p\in P\&\&pd==q}\mu(d)$$
又,根据狄利克雷卷积的基本式
$$g=f u \Rightarrow g(x)=\sum_{a b=x} f(a) u(b)$$
我们发现,上式靠右的如下部分可以作出变换:
$$\sum_{p\in P\&\&pd==q}\mu(d)=g(q)=\sum_{pd=q}[p\in P]\mu(d)$$
若令:
$$f(p)=[p\in P]$$
则我们可以发现:
$$ g(q)=f(p) \mu(d)$$
故而这一部分的值可以用狄利克雷卷积来\(O(nln_n)\)预处理。
一看数据范围,完了:\(10^7\),\(O(nln_n)\)的复杂度显然是不够优秀的。我们需要将它优化到\(O(n)\)

我们考虑对这个函数\(g\)进行线性筛。
我们发现,上式可以进行如下变换:
$$g(q)=\sum_{pd=q\&\&p\in P}\mu(d)=\sum_{p\in P\&\&p|q}\mu(\frac{q}{p})$$
本质上就是枚举其质因子之和并求其除以该质因子的莫比乌斯函数值之和。
我们分三类讨论。
对于\(g(x)\)
如果\(x\in P\),显然\(g(x)=\mu(1)=1\)
如果\(x=a^vb,a\in P,v>1\),则显然枚举其他的质数是没有意义的。\(g(x)=\mu(b)\)
如果\(x=ab,a\in P\),则,\(g(x)\)的值,等价于在计算\(g(b)\)的和的时候每一次计算都多乘上一个\(a\)的贡献,也就是乘上负一。然后将这个值再加上枚举\(a\)的情况,也就是\(\mu(b)\),即\(g(x)=-g(b)+\mu(b)\)
然后大力上一个线性筛即可。

对于左半部分,上一个数论分块即可。

#include<iostream>
#include<cstdio>

const int N = 10000000+5;

int p[N/10],g[N],mu[N],n,m;
bool ip[N];
void prpr(){
	p[0]=0;mu[1]=1;ip[1]=1;
	for(int i=2;i<=10000000;++i){
		if(!ip[i]){
			p[++p[0]]=i;
			mu[i]=-1;
			g[i]=1;
		}
		for(int j=1;j<=p[0]&&p[j]*i<=10000000;++j){
			ip[i*p[j]]=1;
			if(!(i%p[j])){
				g[i*p[j]]=mu[i];
				mu[i*p[j]]=0;
				break;
			}else{
				mu[i*p[j]]=-mu[i];
				g[i*p[j]]=-g[i]+mu[i];
			}
		}
	}
	for(int i=2;i<=10000000;++i){
		g[i]+=g[i-1]; 
	}
}
long long ans;
void init(){
	scanf("%d%d",&n,&m);
	n>m?n^=m^=n^=m:0;
	ans=0;
	int k=0;
	for(int i=1;i<=n;i=k+1){
		k=std::min(n/(n/i),m/(m/i));
		ans+=1ll*(n/i)*(m/i)*(g[k]-g[i-1]);
	}
	printf("%lld\n",ans);
}

int main(){
	prpr();
	int T;
	scanf("%d",&T);
	while(T--){
		init();
	}
	return 0;
}

CF1091 Good Bye 2018

不知不觉就到了2018的最后一天。这个博客也有两个多月了。
我的姿势水平固然有了一定长进,但和巨神的距离却是越拉越远了。
这场CF也是,手速场我却没能有足够的手速。尤其是D题那种超简单的题目却死磕了半天才磕出来,F也没做出来。
总之,祝大家新年快乐_(:з」∠)_


CF1091A
依题意即可。
千万别掉进分类讨论的坑里。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
using namespace std;
int a,b,c;

void init(){
	scanf("%d%d%d",&a,&b,&c);
	--b,c-=2;
	printf("%d",min(min(a,b),c)*3+3);
}
int main(){
	init();
	return 0;
}

CF1091B
我们将每个宝藏都加上每一个向量箭头,并标记目的地。
答案是任意一个被标记n次的。
可以用map或者hashmap维护。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
#include<map>
using namespace std;
int a[1005],b[1005],c[1005],d[1005];
typedef pair<int,int> pii;
typedef map<pii,int> mppii;
mppii mp;
int n;
void init(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i){
		scanf("%d%d",a+i,b+i);
	}
	for(int i=1;i<=n;++i){
		scanf("%d%d",c+i,d+i);
	}
	for(int i=1;i<=n;++i){
		for(int j=1;j<=n;++j){
			pii nw(a[i]+c[j],b[i]+d[j]);
			++mp[nw];
		}
	}
	mppii::iterator it=mp.begin();
	while(it!=mp.end()){
		if(it->second==n){
			printf("%d %d\n",it->first.first,it->first.second);
			break;
		}
		++it;
	} 
}
int main(){
	init();
	return 0;
}

CF1091C
观察题意以后我们发现,对于\(k\),当且仅当\(n\%k == 0\)的时候,\(k\)会和其他的本质不同。
依据剩余系的性质容易证明。
那么我们就得到了一个朴素的求答案公式:
$$\forall x,st:n\% x==0,ans_{x}=\sum_{i=0}^{n/gcd(n,x)-1}(1+(ix-1)\%n+1) $$
然后考虑到\(n\%x==0\),所以
$$ans_{x}=\sum_{i=0}^{n/gcd(n,x)-1}(1+(ix-1)\%n)==\sum_{i=0}^{n/gcd(n,x)-1}(1+ix)$$
套上等差数列求和公式即可。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
using namespace std;
inline long long gcd(long long A,long long B){
	return B?gcd(B,A%B):A;
}
long long n;
long long ans[100005],tp=0;

inline void calc(int X){
	long long tms=n/gcd(n,X);
	ans[++tp]=tms+((((tms*(tms-1))>>1)*X));
}

void init(){
	scanf("%I64d",&n);
	for(int i=1;i*i<=n;++i){
		if(!(n%i)){
			calc(i);
			if(i*i!=n){
				calc(n/i);
			}
		}
	}
	sort(ans+1,ans+1+tp);
	tp=unique(ans+1,ans+1+tp)-1-ans;
	for(int i=1;i<=tp;++i){
		printf("%I64d ",ans[i]);
	}
}
int main(){
	init();
	return 0;
}

CF1091D
由于长度为\(n\),显然序列的值必然是\(\frac{(n*(n+1))}{2}-SUM_i+SUM_{i+1}\), 其中\(SUM_{i}\)表示第\(i\)个排列的某个前缀和。
又排列的性质我们会发现,对于任意两个相邻的全排列,它的任意长度前缀和要么增加且前缀序列变化,要么前缀和保持不变且前缀序列保持不变。
故而,如果一个序列要满足要求,它只有两种情况——它本身就是一个排列,或者它横跨的两个排列必然有一部分前缀相同。
对于第一种情况我们可以看成它们的前缀有\(0\)个相同的。
问题转化为了,求所有排列中,两种相邻排列前缀相同的长度的和。
通过打表找规律,或者对全排列性质的精确推演,我们得到了这样一个求和式:
$$ ans=\sum_{pre=0}^{n-2}(n!-(\Pi_{i=n}^{n-pre+1}i))$$
但是暴力计算的复杂度最坏是\(n^2\)的,仍然无法通过此题。
我们考虑优化这个求和。优化的点在于\(\Pi_{i=n}^{n-pre+1}i\),我们(毫不)惊讶地发现,它等价于\(\frac{n!}{n-pre}\),所以我们可以预处理阶乘和阶乘的逆元。
这样就做完了。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
using namespace std;
const long long MOD=998244353;

long long fac[1000005],inv[1000005];
long long n;

void init(){
	scanf("%I64d",&n);
	fac[0]=fac[1]=inv[0]=inv[1]=1;
	for(int i=2;i<=n;++i){
		fac[i]=fac[i-1]*i%MOD;
		inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
	}
	for(int i=2;i<=n;++i){
		inv[i]=inv[i-1]*inv[i]%MOD;
	}
	long long ans=fac[n];
	for(int i=1;i<=n-2;++i){
		ans+=(fac[n]-fac[n]*inv[n-i]%MOD+MOD)%MOD;
		ans%=MOD;
	}
	printf("%I64d",ans);
}
int main(){
	init();
	return 0;
}

lp5104 红包发红包

在\(\lbrack 0,1 \rbrack \)中任取一个数,期望是\(\frac{1}{2}\)
同理,在\(\lbrack 0,w \rbrack \)中任取一个数,期望是\(\frac{w}{2}\)
故而,取\(k\)个数的期望是\(\frac{w}{2^k}\)
费马小定理加快速幂即可。

#include<iostream>
#include<cstdio>
const long long MOD = 1000000007;

inline long long pw(int A,long long X){
	long long RT=1,BS=A;
	while(X){
		if(X&1){
			RT*=BS;
			RT%=MOD;
		}
		BS*=BS;
		BS%=MOD;
		X>>=1;
	}
	return RT;
}

inline long long inv(int X){
	return pw(X,MOD-2);
}

long long w,n,k;

void init(){
	scanf("%lld%lld%lld",&w,&n,&k);
	printf("%lld\n",(w*(inv(pw(2,k))))%MOD);
}

int main(){
	init();
	return 0;
}

lp5077 Tweetuzki 爱等差数列

令首项为\(x\),长度为\(l\)由等差数列求和式可得:
$$s=\frac{l(l-1)}{2}+l*x$$
故,题意等价于求:
使得式子:
$$2s=l(l-1+2x)$$
成立的最小\(x\)
即求使得式子成立的最大\(l\)
又,易知:
$$l<\sqrt{2s},s\le10^{12}$$
故,
$$l<10^{6}$$
显然这个复杂度是可以接受的。

#include<iostream>
#include<cstdio>
#include<cmath>
long long s;
inline void init(){
	scanf("%lld",&s);
	long long ans=0,len;
	for(long long i=1;i*i<=2*s;++i){
		if((2*s%i==0)&&(2*s/i-i+1>0)&&((s-(i*(i-1)/2))%i==0)){
			len=i;ans=((2*s)/i-i+1)>>1;
		}
	}
	printf("%lld %lld",ans,ans+len-1);
}

int main(){
	init();
	return 0;
} 

lp2260 清华集训2012 模积和

题目大意:求式子:
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j) (i \ne j)$$
的值。

那么容斥以后等价于求:
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j)-\sum_{i=1}^{min(n,m)}(n\ mod\ i)*(m\ mod\ i)$$

首先求:
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j)$$

首先,我们有定理:
$$ n\ mod\ p=n-p*\lfloor \frac{n}{p} \rfloor$$
然后配合上\(\Sigma\)的运算法则,我们可以把原式转化为:
$$n^{2}*m^{2}-m^{2}\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor-n^{2}\sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor-\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor \sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor$$
我们设:
$$ P=\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor,Q=\sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor$$
所以我们可以用数论分块来分别计算\(P,Q\)的值。

而,对于第二部分:
$$R=\sum_{i=1}^{min(n,m)}(n\ mod\ i)*(m\ mod\ i)$$
我们不妨令\(n<m\)
那么将原式展开,可得:
$$R=\sum_{i=1}^{n}n*m-n*i\lfloor \frac{m}{i} \rfloor-m*i\lfloor \frac{n}{i} \rfloor+i^{2}\lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor$$
这个式子显然也是可以数论分块的。但是我们需要计算\(i^{2}\)对于式子的贡献。
用一些简单的证明方法可以得到:
$$\sum_{i=1}^{n}i^{2}=\frac{n(2n+1)(n+1)}{6}$$
那么它的区间和即是:
$$\sum_{i=1}^{r}i^{2}-\sum_{i=1}^{l-1}i^{2}=\frac{(2r^2+2rl+2l^2-r-3l+2)(r-l+1)}{6}$$
于是我们就可以统计整个结果了。

现在我们来考虑数论分块。

数论分块是一种处理形如
$$ \sum_{i=1}^{n}\lfloor \frac{n}{i} \rfloor$$
的公式的求值的方法。

而对于函数:

$$f(x)=\lfloor \frac{a}{x} \rfloor$$

可以发现,对于连续的一段\(x\)的区间,它的值是相同的。
具体来说是,令值为\(k\)则,令区间的左端点为\(l\),则整个区间的值为:
$$k=\lfloor \frac{n}{l} \rfloor$$
那么右端点\(r\)很显然是满足方程:
$$\lfloor \frac{n}{x} \rfloor=k$$
的最大值。
依据整除的性质得到的引理一,我们发现:
$$r=\lfloor \frac{n}{k} \rfloor$$
从而,区间为:
$$[l,\lfloor \frac{n}{\lfloor \frac{n}{l} \rfloor} \rfloor$$

附:对于引理一的证明:

引理一:求证方程:
$$\lfloor \frac{n}{x} \rfloor=k$$
的解的最大值\(t\)为
$$t=\lfloor \frac{n}{k} \rfloor$$
我们令方程的一个合法解\(x\)满足:
$$x*k+r=n\ (0\le r < x)\ (1)$$
并令:
$$t=x+d$$
那么,有:
$$\lfloor \frac{n}{x+d} \rfloor=k$$
则:
$$(x+d)*k+r’=n\ (0\le r < x+d)\ (2)$$
联立方程\((1),(2)\)可得:
$$d*k+r’=r$$
由整除的定义可得:
$$d=\lfloor \frac{r}{k} \rfloor$$
很显然,当
$$x=\lfloor \frac{n}{k} \rfloor$$
时,
$$r=n\ mod\ k,r<k,\lfloor \frac{r}{k} \rfloor=0$$
此时\(x\)取到最大值。
证毕。

另:尽量用模块化编程,而不是去化简式子。这样可以用复杂度换取正确率。

#include<iostream>
#include<cstdio>
const long long MOD = 19940417;
const long long inv2 = 9970209;
const long long inv6 = 3323403;
long long n,m;
inline long long sm1(long long A){
    return (A*(A+1)%MOD)*inv2%MOD;
}
inline long long sm2(long long A){
    return ((A*(A+1)%MOD)*(2*A+1)%MOD)*inv6%MOD;
}
void init(){
    scanf("%lld%lld",&n,&m);
    n>m?n^=m^=n^=m:0; 
    long long P,Q,R;
    int l,r;
    r=0,P=0;
    while(r<n){
        l=r+1;
        r=(n/l)?(n/(n/l)):n;
        P+=(((r-l+1)*(n/l)%MOD)*(l+r)%MOD*inv2)%MOD;
        P%=MOD;
//        printf("P:%lld\n",P);
    }
    P%=MOD;
    r=0,Q=0;
    while(r<m){
        l=r+1;
        r=(m/l)?(m/(m/l)):m;
        Q+=(((r-l+1)*(m/l)%MOD)*(l+r)%MOD*inv2)%MOD;
        Q%=MOD;
//        printf("Q:%lld\n",Q);
    }
//    printf("%lld %lld\n",P,Q);
    Q%=MOD;
    R=(n*m%MOD)*n%MOD;
    long long ans=(n*m%MOD)*(n*m%MOD)%MOD-(n*n%MOD)*Q%MOD-(m*m%MOD)*P%MOD+P*Q%MOD;
    ans%=MOD;
    ans+=MOD;
    ans%=MOD;
    r=0;
    while(r<n){
        l=r+1;
        r=std::min((n/l)?(n/(n/l)):n,(m/l)?(m/(m/l)):m);
        CNT1=m/l;
        CNT2=n/l;
        R-=((((sm1(r)-sm1(l-1))%MOD)*(m/l)%MOD)*n%MOD+(((sm1(r)-sm1(l-1))%MOD)*(n/l)%MOD)*m%MOD);
        R+=(((sm2(r)-sm2(l-1))%MOD)*((m/l)*(n/l)%MOD))%MOD;
        R%=MOD;
    }
    ans-=R;
    ans%=MOD;
    ans+=MOD;
    ans%=MOD;
    printf("%lld",ans);
    
}
int main(){
    init();
    return 0;
}

 

NOIP2018 货币系统

一道看起来像数论题的完全背包问题。
幸好我前一天没有仔细看数论,使得我想的是写暴力…要不然我可能真的会推一个小时的\(EXGCD\)或\(CRT\)或者高斯消元之类的东西。
总之是小凯的疑惑迷惑了我,以为是一道数论题,所以到最后我都以为自己写的是暴力。
结果发现其实是正解思路的,已经写到了正解的倒数第二步,随便一优化就通过了。
但是始终以为是打暴力。
果然还是太菜了。
题意简述:你有\(n\)个正整数,要求,确定一个最小的\(m\),使得用\(m\)种正整数数能够表示的数集与题目给出数能够表示的数集的完全一致。
首先我们证明一个定理:答案的\(m\)个数一定是原有数的子集。
下面我们分两类考虑这个问题。
如果新的数可以被原有数表示,那么显然它是没有意义的。
如果新的数不能被原有数表示,那么显然它是不合法的。
那么我们只需要考虑删除哪些数即可。
我们发现,一个数需要删除,当且仅当它可以被其他数表达出来。
反证法:如果它不可被其他数表示的话,那么删除它以后它就无法表示出来了。
并且,我们可以发现,这个问题满足最优子结构:如果一个数可以被当前的系统表达出来,那么删除当前系统中任意一些可以被表达出来的数,都不会导致这个数不能被表达出来。
故而,先删哪个其实都不影响答案的正确性。
因此我们可以从小往大删(可能可以节省时间)
暴力地考虑,我们需要枚举值域内每一个已经可以被表达出来的数,然后将这个数加上当前数的任意倍(值域范围内)
这么做最坏可能是\(O(T*n*a^2)\)的。
考虑优化。用欧拉筛的思想,对于每一个数,我们只需要加上当前数的一倍即可。
因为,如果加上当前数的两倍的话,就会在之后再一次访问到这个数,这时候就产生了访问次数上的浪费。
这样做是\(O(T*n*a)\)。
那就做完了。

#include<cstdio>
#include<iostream>
#include<algorithm> 
#include<cstring>

int n,a[105];
bool usf[25005];
void init(){
	std::memset(usf,0,sizeof(usf));
	scanf("%d",&n);
	for(int i=1;i<=n;++i){
		scanf("%d",a+i);
	}
	if(n==1){
		puts("1");
		return;
	}
	std::sort(a+1,a+1+n);
	int cnt=n;
	usf[0]=1;
	for(int i=1;i*a[1]<=25000;++i){
		usf[i*a[1]]=1;
	}
	for(int i=2;i<=n;++i){
		if(usf[a[i]]){
			a[i]=0;
			--cnt;
		}else{
			for(int j=0;j+a[i]<=25000;++j){
				if(!usf[j]){
					continue;
				}else{
					usf[j+a[i]]=1;
				}
			}
		}
	}
	printf("%d\n",cnt);
}
int main(){
	int T;
	scanf("%d",&T);
	while(T--){
		init();
	}
	return 0;
}

 

线性求逆元

$$我们有一个质数p$$
$$\forall i\in Z,def\ q_{i}=[\frac{p}{i}],r_{i}=p\%i\ st:\ p = iq_{i}+r_{i} \ (1)$$
$$def\ i^{-1}*i≡1\ (mod\ p)$$
$$(1)*i^{-1}$$
$$=>p*i^{-1}≡q_{i}+r_{i}*i^{-1}\ (mod\ p)$$
$$∵(p*i^{-1})\%p≡0\ (mod\ p)$$
$$∴r_{i}*i^{-1}≡-q_{i}\ (mod\ p)$$
$$i^{-1}≡-q_{i}*r_{i}^{-1}\ (mod\ p)$$
$$i^{-1}≡(p-[\frac{p}{i}])*(p\%i)^{-1}\ (mod\ p)$$
代码非常简短:

#include<iostream>
#include<cstdio>
int MOD;
int inv[1000005],n;
int main(){
	scanf("%d%d",&n,&MOD);
	inv[1]=1;
	for(int i=2;i<=n;++i){
		inv[i]=(MOD-(MOD/i))*inv[MOD%i]%MOD;
	}
	for(int i=1;i<=n;++i){
		printf("%d ",inv[i]); 
	}
}

 

lp2822 NOIP2016 组合数问题

基础数论题。
首先我们知道,任何数的逆元模k,都不可能等于零。
故而我们不必考虑k是否是质数。
然后我们考虑先预处理出哪些\(n,m\)满足\(C_{n}^{m}≡0(mod\ k)\)
接着前缀和即可。
由于\(n,m\)较小,可以考虑用递推。

#include<iostream>
#include<cstdio> 
#include<cstring>
using namespace std;
int sm[2005][2005],C[2005][2005];
int n,m,k;
void prpr(){
    memset(sm,0,sizeof(sm));
    C[0][0]=C[1][1]=1;
    for(int i=1;i<=2000;++i){
        C[i][0]=1,C[i][i]=1;
    }
    for(int i=2;i<=2000;++i){
        for(int j=1;j<=i;++j){
            C[i][j]=(C[i-1][j-1]+C[i-1][j])%k;
        }
    }
    for(int i=2;i<=2000;++i){
        for(int j=1;j<=i;++j){
            sm[i][j]=sm[i-1][j]+sm[i][j-1]-sm[i-1][j-1];
            if(C[i][j]==0){
                ++sm[i][j];
            }
        }
        sm[i][i+1]=sm[i][i];
    }
}
void init(){
    scanf("%d%d",&n,&m);
    printf("%d\n",(m>n)?(sm[n][n]):(sm[n][m]));
}
int main(){
    int T;
    scanf("%d%d",&T,&k);
    prpr();
    while(T--){
        init();
    }
    return 0;
}

 

Lucas定理

首先给出关于Lucas定理的简要证明:
定义:
$$ a=a_{k}*p^{k}+a_{k-1}*p^{k-1}+…+a_{1}*p+a_{0}=\sum_{i=1}^{k}a_{i}^p{i} $$
$$ b=b_{k}*p^{k}+b_{k-1}*b^{k-1}+…+b_{1}*p+b_{0}=\sum_{i=1}^{k}b_{i}^p{i} $$
求证:
$$ C_{a}^{b}=C_{a_{k}}^{b_{k}}*C_{a_{k-1}}^{b_{k-1}}…C_{a_{0}}^{b_{0}} $$

首先我们证明引理一:
$$ (1+x)^p≡1+x^p\ (mod\ p) $$
根据组合数基本性质,我们有:
$$\forall j\in [1,p-1],C_{p}^{j}=\frac{p}{j}*C_{p-1}^{j-1}≡0(mod\ p) $$

$$∴(1+x)^p≡\sum_{i=1}^{p}C_{p}^{i}*x^i≡1+x^p(mod\ p) $$
于是我们得到结论:
$$(1+x)^a≡\prod_{i=1}^{k}(1+x)^{p^k*a_{k}}≡\prod_{i=1}^{k}(1+x^{p^{k}})^{a_{k}}(mod\ p)\ (1)$$
根据进制的基本性质和幂的基本性质,我们有:
$$b=\sum_{i=1}^{k}b_{i}^p{i},x^b=\prod_{i=1}^{k}x^{p^k*b_{i}}$$
并且我们知道,用上述方法表示\(p\)进制数,是完全等价的。即,两者的集合构成双射。

故而我们比较\((1)\)式展开后左右各项,可以得到:
$$\forall b\in [1,a],C_{a}^{b}≡\prod_{i=1}^{k}C_{a_{i}}^{b_{i}}(mod\ p)$$

证毕。

故而对于实际处理问题,只需要逆向秦九韶即可。

 

#include<iostream>
#include<cstdio>
using namespace std;
long long n,m,p,inv[100005],fac[100005];
#define Max(_A,_B) ((_A)>(_B)?(_A):(_B))
#define Min(_A,_B) ((_A)<(_B)?(_A):(_B))

//y取x 
inline long long C0(long long x,long long y){
    return ((x>y)?0:(x?fac[y]*inv[x]*inv[y-x]%p:1))%p;
}
//模p意义下的y取x 
inline long long C(long long x,long long y){
    return ((x>y)?0:((y>=p)?C(x/p,y/p)*C0(x%p,y%p):C0(x%p,y%p)))%p;
}
void init(){
    scanf("%lld%lld%lld",&n,&m,&p);
    fac[0]=fac[1]=inv[0]=inv[1]=1;
    for(int i=2;i<=p;++i){
        fac[i]=fac[i-1]*i%p;
        inv[i]=(p-p/i)*inv[p%i]%p;
    }
    for(int i=2;i<=p;++i){
        inv[i]*=inv[i-1];
        inv[i]%=p;
    }
    printf("%lld\n",C(n,n+m)%p);
} 
int main(){
    int T;
    scanf("%d",&T);
    while(T--){
        init();
    }
    return 0;
}