uoj#272. 【清華集訓2016】石家莊的工人階級隊伍比較堅強（矩陣+三維FWT）

傳送門

題解

抄代碼\(20\)分鐘，搞懂題解在幹嗎仨小時→_→

到今天才算真正搞明白\(FWT\)在幹嘛了

本題

首先轉移關係都是恆定的，設它爲一個矩陣\(B\)，那麼要求的就是\(f_n=f_0B^n\)

定義三進制不退位減法\(\ominus\)，三進制不退位加法\(\oplus\)，這兩個互爲逆運算（能夠類比一下二進制的異或）

根據\(B\)矩陣的定義以及剪刀石頭布的性質易知\(\forall k \lt 3^m,B_{i\oplus k,j\oplus k}=B_{i,j}\)

那麼易知\(B_{i,j}=B_{i\ominus j,0}\)。不難發現這個結論也能夠推廣到\(B^n\)，即有\(\forall k \lt 3^m,B^n_{i\oplus k,j\oplus k}=B^n_{i,j}\)，能夠用概括法證，也能夠感性理解一下

那麼對於\(f_n=f_0B^n\)，它的第\(i\)項就爲\[f_{n,i}=\sum_k f_{0,k}B^n_{k,i}=\sum_k f_{0,k}B^n_{0,i\ominus k}=\sum_{x\oplus y=i}f_{0,x}B^n_{0,y}\]
因而這玩意兒能夠寫成一個卷積的形式

咱們發現上面的式子只須要\(B\)的第一行就夠了，那麼由於\[B^n_{0,i}=\sum_k B^n_{0,k}B_{k,i}=\sum_k B^n_{0,k}B_{0,i\ominus k}=\sum_{x\oplus y=i} B^n_{0,x}B_{0,y}\]
那麼\(B\)的第一行作\(n\)次卷積，而後\(f\)再和它作一次卷積就是答案了

卷積

咱們如今須要找到一個三進制不退位加法的卷積變換。據LargestJN大佬說，若是下標運算是模意義下加法，這個卷積就叫作循環卷積

然而就算它叫雞卷也沒用咱們仍是不會求

首先一個卷積就是要知足\(T(a)T(b)=T(a\oplus b)\)（這裏\(\oplus\)爲任意運算），設\(T\)爲這個卷積的矩陣，\(x_{i,j}\)爲這個矩陣中的某個元素，那麼上面的形式就能夠表現爲\[(\sum_{j=0}^{n-1}x_{i,j}a_j)(\sum_{k=0}^{n-1}x_{i,k}a_k)=\sum_{j=0}^{n-1}x_{i,j}\sum_{k\oplus t=j}a_kb_t\]
因而考慮\(a_k\)和\(b_t\)的貢獻，得有\[x_{i,k}x_{i,t}a_kb_t=x_{i,k\oplus t}a_k,b_t\]
那麼就是對於\(T\)的每一行都得有\[x_{k}x_{t}=x_{k\oplus t}\]
因而\(T\)的每一行都是方程組\(x_{k}x_{t}=x_{k\oplus t}\)的一組解

然而只有卷積萬一沒有卷積的逆變換（管它叫雞卷好了）就gg了

因此咱們的\(T\)還得有逆矩陣，那麼根據線性代數的芝士，\(T\)的行列式不爲\(0\)，那麼方程組\(x_{k}x_{t}=x_{k\oplus t}\)至少要有\(n\)組不一樣的解

先來看看循環卷積的通常形式，長這樣\[A\times B=\sum_{i}\sum_{(j+k)mod \ n=i}A_jB_k\]
\(T\)知足\(x_ix_j=x_{(i+j)mod\ n}\)

而後咱們發現\(n\)次單位根\(\omega^i\)就是一組可行的解

不知道\(n\)次單位根的能夠回去看看\(FFT\)的原理

然而須要\(n\)組解，因而矩陣長這樣

逆矩陣爲\[\frac{1}{n} \begin{bmatrix} 1& 1 & 1& ... & 1\\ 1& w_n^{-1}& w_n^{-2}& ... & w_n^{-(n - 1)}\\ 1& w_n^{-2} & w_n^{-4}& ... & w_n^{-2(n- 1)}\\ ...& ...& ...& ...& ...\\ 1& w_n^{-(n - 1)}& w_n^{-2(n - 1)} & ... & w_n^{-(n - 1)(n - 1)} \end{bmatrix}\]

回到本題

由於是三進制不進位加法，因此就選三次單位根\(\omega\)。這裏能夠把全部的複數都表示爲\(a+b\omega\)的形式，那麼由於三次單位根有\(\omega^2+\omega+1=0\)，因此有\(\omega^2=-\omega-1\)，因此全部的係數都是\(0/1/-1\)，運算過程當中就不用取模了。兩個複數的乘法就是\[(a+b\omega)(c+d\omega)=ac+(ad+bc)\omega+bd(-\omega-1)=(ac-bd)+(ad+bc-bd)\omega\]
而後根據上面，\[T= \left[ \begin{matrix} 1 & 1 & 1 \\ 1 & \omega & \omega^2 \\ 1 & \omega^2 & \omega \end{matrix} \right]\]
以及咱們知道\[\left[ \begin{matrix} 1 & 1 & 1 \\ 1 & \omega & \omega^2 \\ 1 & \omega^2 & \omega \end{matrix} \right]\times \left[ \begin{matrix} 1 & 1 & 1 \\ 1 & \omega^2 & \omega \\ 1 & \omega & \omega^2 \end{matrix} \right]=3E\]
其中\(E\)是單位矩陣

那麼咱們就能夠把後面那個當作逆矩陣，帶進去算

由於\(DFT\)和\(IDFT\)本質都是分治作向量對矩陣的乘法，它\(IDFT\)的時候每一次都多乘了個\(3\)，總共有\(\log_3 n\)層，那麼總共多乘了\(3^{\log_3 n}=n\)次，那麼只要最後把全部元素都除以一個\(n\)就吼啦！

奇怪的模數問題……

最後是一個比較奇怪的模數問題，由於咱們最後要除以\(3^m\)，要求\(3^m\)有逆元，也就是\(3\nmid p\)，而後由於它這個奇怪的模數\(p\)，假設\(k=\frac{p}{3}\)是個正整數，那麼有\[\frac{1}{k + 1} + \frac{1}{k(k + 1)} = \frac{1}{k} = \frac{3}{p}\]
矛盾，因而\(3\nmid p\)

而後……就真的沒而後了……

//minamoto
#include<cstdio>
#include<cstring>
#include<map>
#define R register
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
char sr[1<<21],z[20];int C=-1,Z=0;
inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
void print(R int x){
    if(C>1<<20)Ot();if(x<0)sr[++C]='-',x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++C]=z[Z],--Z);sr[++C]='\n';
}
const int N=6e5+5;
int lim,m,P,t,b[25][25],a[N],inv2,invn;
inline int add(R int x){return x>=P?x-P:x;}
inline int dec(R int x){return x<0?x+P:x;}
inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
void exgcd(int &x,int &y,int a,int b){
    if(!b)return (void)(x=1,y=0);
    exgcd(y,x,b,a%b),y-=a/b*x;
}
inline int inv(int a){
    int x,y;exgcd(x,y,a,P);
    return (x%P+P)%P;
}
struct complex{
    int x,y;
    complex(int X=0,int Y=0):x(X),y(Y){}
    inline complex operator +(const complex &b){return complex(add(x+b.x),add(y+b.y));}
    inline complex operator -(const complex &b){return complex(dec(x-b.x),dec(y-b.y));}
    inline complex operator *(const complex &b){
        return complex(dec(mul(x,b.x)-mul(y,b.y)),dec(add(mul(x,b.y)+mul(y,b.x))-mul(y,b.y)));
    }
    inline bool operator <(const complex &b)const{
        return x==b.x?y<b.y:x<b.x;
    }
}f[N],g[N];map<complex,complex>mp;
complex ksm(complex x,R int y){
    complex res(1,0);
    for(;y;y>>=1,x=x*x)if(y&1)res=res*x;
    return res;
}
complex calc1(complex b){return complex(dec(-b.y),dec(b.x-b.y));}
complex calc2(complex b){return complex(dec(b.y-b.x),dec(-b.x));}
void DFT(complex *A){
    for(R int mid=1;mid<lim;mid*=3)
        for(R int j=0;j<lim;j+=mid*3)
            for(R int k=0;k<mid;++k){
                complex x=A[j+k],y=A[j+k+mid],z=A[j+k+(mid<<1)];
                A[j+k]=x+y+z;
                A[j+k+mid]=x+calc1(y)+calc2(z);
                A[j+k+(mid<<1)]=x+calc2(y)+calc1(z);
            }
}
void IDFT(complex *A){
    for(R int mid=1;mid<lim;mid*=3)
        for(R int j=0;j<lim;j+=mid*3)
            for(R int k=0;k<mid;++k){
                complex x=A[j+k],y=A[j+k+mid],z=A[j+k+(mid<<1)];
                A[j+k]=x+y+z;
                A[j+k+mid]=x+calc2(y)+calc1(z);
                A[j+k+(mid<<1)]=x+calc1(y)+calc2(z);
            }
    for(R int i=0;i<lim;++i)A[i].x=mul(A[i].x,invn);
}
int main(){
//  freopen("testdata.in","r",stdin);
    m=read(),t=read(),P=read();
    lim=1;fp(i,1,m)lim*=3;
    if(P==1){
        fp(i,0,lim-1)print(0);
        return Ot(),0;
    }
    fp(i,0,lim-1)a[i]=read();
    fp(i,0,m)fp(j,0,m-i)b[i][j]=read();
    fp(i,0,lim-1){
        int tmp=i,cntw=0,cntl=0;
        while(tmp){
            int k=tmp%3;
            k==1?++cntw:k==2?++cntl:0;
            tmp/=3;
        }g[i]=b[cntw][cntl],f[i]=a[i];
    }
    inv2=inv(2),invn=inv(lim);
    DFT(f),DFT(g);
    fp(i,0,lim-1)f[i]=f[i]*(mp.count(g[i])?mp[g[i]]:(mp[g[i]]=ksm(g[i],t)));
    IDFT(f);
    fp(i,0,lim-1)print(f[i].x);
    return Ot(),0;
}