题意
(n) 局石头剪刀布，设每局的贡献为赢的次数与输的次数之 (gcd) ，求期望贡献乘以 (3^{2n}) ，定义若 (xy=0) 则，(gcd(x,y)=x+y)
思路
不难得出
[ans=3^nsum_{i=0}^nsum_{j=0}^{n-i}{nchoose i}{n-ichoose j}gcd(i,j)
]
对于正整数 (n) ，有如下表达式
[n=sum_{d|n}varphi(d)
]
那么
[ans=3^nsum_{i=1}^nsum_{j=1}^{n-i}{nchoose i}{n-ichoose j}sum_{d|gcd(i,j)}varphi(d)+3^ncdot2sum_{i=1}^n{nchoose i}i
]
注意欧拉函数是不取到零的，若要严谨的写，零应该特判。
在外层枚举 (d) ，并只枚举为 (d) 倍数的 (i,j) ，
[ans=3^nsum_{d=1}^nvarphi(d)sum_{i=1}^{lfloor{nover d}
floor}sum_{j=1}^{lfloor{nover d}
floor-i}{nchoose id}{n-idchoose jd}+3^ncdot2sum_{i=1}^n{nchoose i}i
]
另外，利用表达式
[{nchoose m}={nover m}{n-1choose m-1}
]
还可以化简后式（虽然对复杂度没有影响。。）
[ans=3^nn!sum_{d=1}^nvarphi(d)sum_{i=1}^{lfloor{nover d}
floor}sum_{j=1}^{{lfloor{nover d}
floor}-i}{1over (id)!}{1over (jd)!}cdot{1over(n-id-jd)!}+6^nn
]
变成了卷积的形式，后面于 (i+j) 有关的项就当作最后乘出多项式的常数。
复杂度的证明如下
[egin{array}{}
T(n)&=displaystylesum_{i=1}^n {nover i}log{nover i}\
&geq displaystylesum_{i=1}^n {nover i}log{n}\
&= log ndisplaystylesum_{i=1}^n {nover i}\
&geq nlog nln n
end{array}
]
那么此题可解。
代码
#include<bits/stdc++.h>
#define FOR(i,x,y) for(int i=(x),i##END=(y);i<=i##END;++i)
#define DOR(i,x,y) for(int i=(x),i##END=(y);i>=i##END;--i)
using namespace std;
template<typename T,typename _T>inline bool chk_min(T &x,const _T y){return y<x?x=y,1:0;}
template<typename T,typename _T>inline bool chk_max(T &x,const _T y){return x<y?x=y,1:0;}
typedef long long ll;
const double PI=acos(-1.0);
const int N=1<<17|5;
namespace _Maths
{
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
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 Pow(ll a,ll p,ll P)
{
ll res=1;
for(;p>0;p>>=1,(a*=a)%=P)if(p&1)(res*=a)%=P;
return res;
}
ll inv(ll a,ll P){ll x,y;exgcd(a,P,x,y);return (x%P+P)%P;}
};
using namespace _Maths;
int P;
struct Complex
{
double x,y;
Complex operator +(const Complex &_)const{return (Complex){x+_.x,y+_.y};}
Complex operator -(const Complex &_)const{return (Complex){x-_.x,y-_.y};}
Complex operator *(const Complex &_)const{return (Complex){x*_.x-y*_.y,x*_.y+y*_.x};}
Complex operator /(const int &_)const{return (Complex){x/_,y/_};}
};
namespace _Polynomial
{
const int K=(1<<15)-1,L=15;
Complex A[N<<1],B[N<<1],C[N<<1],D[N<<1];
Complex w[N<<1];int r[N<<1];
void DFT(Complex *a,int op,int n)
{
FOR(i,0,n-1)if(i<r[i])swap(a[i],a[r[i]]);
for(int i=2;i<=n;i<<=1)
for(int j=0;j<n;j+=i)
for(int k=0;k<i/2;k++)
{
Complex u=a[j+k],t=w[op==1?n/i*k:(n-n/i*k)&(n-1)]*a[j+k+i/2];
a[j+k]=u+t,a[j+k+i/2]=u-t;
}
if(op==-1)FOR(i,0,n-1)a[i]=a[i]/n;
}
void multiply(const int *a,const int *b,int *c,int n1,int n2)
{
int n=1;
while(n<n1+n2-1)n<<=1;
FOR(i,0,n1-1)A[i]=(Complex){a[i]&K,a[i]>>L};
FOR(i,0,n2-1)B[i]=(Complex){b[i]&K,b[i]>>L};
FOR(i,n1,n-1)A[i]=(Complex){0,0};
FOR(i,n2,n-1)B[i]=(Complex){0,0};
FOR(i,0,n-1)r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
FOR(i,0,n-1)w[i]=(Complex){cos(2*PI*i/n),sin(2*PI*i/n)};
DFT(A,1,n),DFT(B,1,n);
FOR(i,0,n-1)
{
int j=(n-i)&(n-1);
C[i]=(Complex){0.5*(A[i].x+A[j].x),0.5*(A[i].y-A[j].y)}*B[i];
D[i]=(Complex){0.5*(A[i].y+A[j].y),0.5*(A[j].x-A[i].x)}*B[i];
}
DFT(C,-1,n),DFT(D,-1,n);
FOR(i,0,n1+n2-2)
{
ll s=C[i].x+0.5,t=C[i].y+0.5,u=D[i].x+0.5,v=D[i].y+0.5;
c[i]=(s+((t+u)%P<<L)%P+(v%P<<L<<L)%P)%P;
}
}
};
int A[N],B[N],C[N<<1];
int phi[N],fac[N],ifac[N];
int n;
void init()
{
fac[0]=fac[1]=1;FOR(i,2,n)fac[i]=(ll)fac[i-1]*i%P;
ifac[0]=ifac[1]=1;FOR(i,2,n)ifac[i]=(ll)(P-P/i)*ifac[P%i]%P;
FOR(i,2,n)ifac[i]=(ll)ifac[i-1]*ifac[i]%P;
FOR(i,1,n)phi[i]=i;
FOR(i,2,n)if(phi[i]==i)
for(int j=i;j<=n;j+=i)
phi[j]=phi[j]/i*(i-1);
}
int main()
{
int T;
scanf("%d",&T);
while(T--)
{
scanf("%d%d",&n,&P);
init();
ll ans=0;
FOR(d,1,n)
{
ll sum=0;
FOR(i,1,n/d)A[(i)-1]=ifac[i*d];
FOR(i,1,n/d)B[(i)-1]=ifac[i*d];
_Polynomial::multiply(A,B,C,n/d,n/d);
FOR(i,1,n/d)(sum+=(ll)C[(i)-2]*ifac[n-i*d]%P)%=P;
(ans+=sum*phi[d]%P)%=P;
}
ans=(ans*Pow(3,n,P)%P*fac[n]%P+Pow(6,n,P)*n%P)%P;
printf("%lld
",(ans+P)%P);
}
return 0;
}