Link
就是(K)进制FWT。
注意到题目保证(
otexists frac1x+frac1y=frac3p),那么(2^{-1}mod p,3^{-1}mod p)都是存在的。
一个小结论是转移矩阵DFT之后,本质不同的值最多只有(frac{m(m+1)}2)项,那么我们可以用一个map记忆化来优化复杂度。
最后时间复杂度为(O(3^mm+m^2log t))。
#include<map>
#include<cstdio>
#include<cctype>
const int N=531441;
int P,n,m,t,v[13][13];
struct complex{int x,y;complex(int a=0,int b=0):x(a),y(b){}}a[N],b[N];
int read(){int x=0,c=getchar();while(isspace(c))c=getchar();while(isdigit(c))(x*=10)+=c&15,c=getchar();return x;}
int inc(int a,int b){return a+=b-P,a+(a>>31&P);}
int dec(int a,int b){return a-=b,a+(a>>31&P);}
int mul(int a,int b){return 1ll*a*b%P;}
void exgcd(int&x,int&y,int a,int b){!b? (x=1,y=0):(exgcd(y,x,b,a%b),y-=a/b*x);}
int inv(int a){int x,y;exgcd(x,y,a,P);return (x+P)%P;}
int operator<(complex a,complex b){return a.x<b.x||(a.x==b.x&&a.y<b.y);}
complex operator+(complex a,complex b){return {inc(a.x,b.x),inc(a.y,b.y)};}
complex operator*(complex a,complex b){return {dec(mul(a.x,b.x),mul(a.y,b.y)),dec(inc(mul(a.x,b.y),mul(a.y,b.x)),mul(a.y,b.y))};}
complex operator^(complex a,int k){complex r(1);for(;k;k>>=1,a=a*a)if(k&1)r=a*r;return r;}
complex F(complex a,int f){return ~f? complex{dec(0,a.y),dec(a.x,a.y)}:complex{dec(a.y,a.x),dec(0,a.x)};}
void FWT(complex*a,int f)
{
for(int i=1;i<n;i*=3)
for(int j=0;j<n;j+=i*3)
for(int k=0;k<i;++k)
{
complex x=a[k+j],y=a[i+j+k],z=a[i+i+j+k];
a[j+k]=x+y+z,a[i+j+k]=x+F(y,f)+F(z,-f),a[i+i+j+k]=x+F(y,-f)+F(z,f);
}
if(!~f) for(int i=0,x=inv(n);i<n;++i) a[i].x=mul(a[i].x,x);
}
void dfs(int x,int id,int w,int l)
{
if(x==m) return b[id].x=v[w][l],void();
id*=3,dfs(x+1,id,w,l),dfs(x+1,id+1,w+1,l),dfs(x+1,id+2,w,l+1);
}
int main()
{
m=read(),t=read(),P=read(),n=1;std::map<complex,complex>mp;
for(int i=0;i<m;++i) n*=3;
for(int i=0;i<n;++i) a[i].x=read();
for(int i=0;i<=m;++i) for(int j=0;i+j<=m;++j) v[i][j]=read();
dfs(0,0,0,0),FWT(a,1),FWT(b,1);
for(int i=0;i<n;++i) a[i]=a[i]*(mp.count(b[i])? mp[b[i]]:mp[b[i]]=b[i]^t);
FWT(a,-1);
for(int i=0;i<n;++i) printf("%d
",a[i].x);
}