HDU-6088 Rikka with Rock-paper-scissors (容斥+MTT)

HDU-6088(容斥+MTT)

考虑计算两个人赢得次数（设为\(a,b\)）都是\(d\)的倍数的方案数

注意一下\(a+b>0\)

对于\(a,b\)，它的方案数为\(C(n,a)C(n-a,b)\)，即\(\frac{n!}{a!b!(n-a-b)!}\)

多以对于每个\(d\)计算一遍所有可行的\(a,b\)能组成的总方案数

这个可以\(MTT\)(因为模数不确定)

\(MTT\)的实现，可以参考，在文章下面的扩展介绍里

由于是\(d\)的倍数会包含所有\(d|gcd(a,b)\)的情况，所以还需要容斥一下

复杂度\(n\ln n \log n\)，以及极其大的常数

#include<bits/stdc++.h>
using namespace std;

#define double long double

#define reg register
typedef long long ll;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)

template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); } 
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); } 

char IO;
int rd(){
	int s=0,f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const double PI=acos((double)-1);

const int N=(1<<18)+4;
ll S;

bool be;

int n,P;
ll qpow(ll x,ll k,ll P) {
	x%=P;
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}
ll qmul(ll x,ll y,ll P){
	y=(y%P+P)%P;
	ll res=0;
	for(;y;y>>=1,x+=x,(x>=P&&(x-=P))) if(y&1) res+=x,(res>=P&&(res-=P));
	return res;
}

int rev[N];
struct Cp{
	double x,y;
	Cp operator + (const Cp t) const { return (Cp){x+t.x,y+t.y}; }
	Cp operator - (const Cp t) const { return (Cp){x-t.x,y-t.y}; }
	Cp operator * (const Cp t) const { return (Cp){x*t.x-y*t.y,x*t.y+y*t.x}; }
}A[N],B[N],w[N];

void FFT(int n,Cp *a,int f){
	rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
	for(reg int i=1;i<n;i<<=1) {
		int len=n/i/2;
		for(reg int l=0;l<n;l+=2*i) {
			for(reg int j=l;j<l+i;j++) {
				Cp t=a[j+i]*w[len*(j-l)];
				a[j+i]=a[j]-t;
				a[j]=a[j]+t;
			}
		}
	}
	if(f==-1) rep(i,0,n-1) a[i].x/=n,a[i].y/=n;
}


ll C[N],Inv[N],Fac[N];

ll Down(double x){
	return x<0?(ll)(x-0.5):(ll)(x+0.5);
}//四舍五入取整
bool ed;
int main(){
	rep(kase,1,rd()) {
		n=rd(),P=rd(),S=sqrt(P)+1;
		Inv[0]=Inv[1]=Fac[0]=Fac[1]=1;
		rep(i,2,n) {
			Inv[i]=(P-P/i)*Inv[P%i]%P;
			Fac[i]=Fac[i-1]*i%P;
		}
		rep(i,1,n) Inv[i]=Inv[i-1]*Inv[i]%P;
		int lst=-1;
		rep(d,1,n) {
			int lim,R=1,c=-1;
             //MTT
			for(lim=0;lim*d<=n;lim++) {
				ll t=Inv[lim*d];
				A[lim]=(Cp){1.0*(t/S),1.0*(t%S)};
				B[lim]=(Cp){1.0*(t%S),0.0};
			}
			lim--;
			while(R<=lim*2) R<<=1,c++;
			rep(j,lim+1,R) A[j].x=A[j].y=B[j].x=B[j].y=0;

			if(lst!=R) {
				rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
				w[0]=(Cp){1,0},w[1]=(Cp){cos(2*PI/R),sin(2*PI/R)};
				rep(i,1,R-1) if(i%5000==3) w[i]=(Cp){cos(2*PI/R*i),sin(2*PI/R*i)};
				else w[i]=w[i-1]*w[1];
				lst=R;
			} else rep(i,0,R-1) w[i].y=-w[i].y;

			FFT(R,A,1),FFT(R,B,1);
			rep(i,0,R-1) A[i]=A[i]*A[i],B[i]=B[i]*B[i];

			rep(i,0,R-1) w[i].y=-w[i].y;
			FFT(R,A,-1),FFT(R,B,-1);
			C[d]=0;
			rep(i,1,lim) {
				ll a=Down(A[i].x),b=Down(B[i].x),ab=Down(A[i].y),t=0;
				a+=b;
				t=a%P*S%P*S%P+ab%P*S%P+b%P;
				t=t*Fac[n]%P*Inv[n-d*i]%P;
				C[d]=(C[d]+t)%P;
			}
		}
		drep(i,n,1)	for(reg int j=i+i;j<=n;j+=i) C[i]=(C[i]-C[j])%P;//容斥
		ll ans=0;
		rep(i,1,n) ans=(ans+C[i]*i)%P;
		rep(i,1,n) ans=ans*3%P;
		ans=(ans%P+P)%P;
		printf("%lld\n",ans);
	}
}