[uoj272]石家庄的工人阶级队伍比较坚强
假设$x,y\in \{0,1,2\}$,则$x$能赢$y$(根据题中定义)当且仅当$x-y\equiv 1(mod\ 3)$
定义$\ominus$为两数3进制下不退位的减法,$S_{x}$表示$x$在3进制下1的个数,则$u=W(x,y)=S_{x\ominus y}$,类似地,有$v=W(y,x)=S_{y\ominus x}$
记$B_{x,y}=b_{S_{x\ominus y},S_{y\ominus x}}$,令$k=x\ominus y$,根据$\ominus$的性质,有$k\ominus 0=x\ominus y$以及$0\ominus k=y\ominus x$,代入后就可以得到$B_{x,y}=b_{S_{k\ominus 0},S_{0\ominus k}}=B_{k,0}$,再记$g_{k}=B_{k,0}$,可以预处理出来
原题也就是$f_{i,x}=\sum_{y=0}^{3^{m}-1}g_{x\ominus y}f_{i-1,y}$,再令$\oplus$为3进制下按位加法,有$(x\ominus y)\oplus y=x$,也就是一个$\oplus$的卷积,记这个卷积为$\otimes$运算,有$f_{i}=g\otimes f_{i-1}$
根据卷积的结合律,有$f_{t}=g^{t}\otimes f_{0}$(其中$g^{t}$也是上面的$\otimes$),也就是要快速计算$\otimes$运算
我们直接来考虑原问题的一个拓展,即令$\oplus$为$B$进制下不进位加法,$\otimes$为$\oplus$卷积($B$是一个较小的数),如何对两个序列快速计算$\otimes$这个运算
为了方便起见,不妨假设序列长度$n=B^{k}$,不足可以补0
对于两个序列$a_{i}$和$b_{i}$,考虑构造一个$n\times n$的矩阵$A_{i,j}$(将序列看作一个$1\times n$的矩阵),令$a'=aA$,$b'=bA$,再令$c'=a'*b'$(这里的*指对应位置相乘),最后$c=c'A^{-1}$
$\forall n=B^{k}(k\ge 1)$,若矩阵$A$满足以下两个条件,则最终有$c=a\otimes b$:
1.其左上角$B\times B$的矩阵在$n=B$时满足$c=a\otimes b$
2.$A_{i,j}=A_{\lfloor\frac{i}{B}\rfloor,\lfloor\frac{j}{B}\rfloor}A_{i\ mod\ B,j\ mod\ B}$,且$A^{-1}$具有同样性质
证明如下——
注意到当$n=B$时,左上角$B\times B$的矩阵即为$A$本身,那么由条件1可知其正确
归纳$n=B^{k}(k\ge 1)$时正确,考虑$n=B^{k+1}$时正确性:
对于每一个$0\le i<B^{k+1}$,都可以被表示为$i_{0}B+i_{1}$(其中$0\le i_{0}<B^{k}$,$0\le i_{1}<B$)
对于$a_{i}$和$b_{i}$,代入式子即有
$$
c_{i}=\sum_{j=0}^{B^{k+1}-1}c'_{j}A^{-1}_{j,i}=\sum_{j_{1}=0}^{B-1}\sum_{j_{0}=0}^{B^{k}-1}c'_{j_{0}B+j_{1}}A^{-1}_{j_{0}B+j_{1},i_{0}B+i_{1}}=\sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}\sum_{j_{0}=0}^{B^{k}-1}c'_{j_{0}B+j_{1}}A^{-1}_{j_{0},i_{0}}
$$
考虑$c'_{i}$,同样可以代入得到
$$
c'_{i}=a'_{i}b'_{i}=(\sum_{j_{1}=0}^{B-1}A_{j_{1},i_{1}}\sum_{j_{0}=0}^{B^{k}-1}a_{j_{0}B+j_{1}}A_{j_{0},i_{0}})(\sum_{j_{1}=0}^{B-1}A_{j_{1},i_{1}}\sum_{j_{0}=0}^{B^{k}-1}b_{j_{0}B+j_{1}}A_{j_{0},i_{0}})
$$
(关于$a'_{i}$和$b'_{i}$用$a_{i}$和$b_{i}$表达,与$c_{i}$用$c'_{i}$表达相同,将$A^{-1}$改为$A$即可)
对$a_{i}$和$b_{i}$分别构造出$B$个序列,即$d_{j_{0},i}=a_{iB+j_{0}}$以及$e_{j_{0},i}=b_{iB+j_{0}}$,考虑$d_{x}$和$e_{y}$分别作为$a_{i}$和$b_{i}$通过前面的方式计算,可以算出$c'$,记作$f'_{x,y}$,即有
$$
f'_{x,y,i}=(\sum_{j=0}^{B^{k}-1}d_{x,j}A_{j,i})(\sum_{j=0}^{B^{k}-1}e_{y,j}A_{j,i})=(\sum_{j=0}^{B^{k}-1}a_{jB+x}A_{j,i})(\sum_{j=0}^{B^{k}-1}b_{jB+y}A_{j,i})
$$
可以发现这个式子与上面的$c'_{i}$有关,具体来说即
$$
c'_{i}=\sum_{x=0}^{B-1}A_{x,i_{1}}\sum_{y=0}^{B-1}A_{y,i_{1}}f'_{x,y,i_{0}}
$$
整体将$c'_{i}$的式子代入$c_{i}$中,也就得到了
$$
c_{i}=\sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}\sum_{j_{0}=0}^{B^{k}-1}A^{-1}_{j_{0},i_{0}}\sum_{x=0}^{B-1}A_{x,j_{1}}\sum_{y=0}^{B-1}A_{y,j_{1}}f'_{x,y,j_{0}}=\sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}\sum_{x=0}^{B-1}A_{x,j_{1}}\sum_{y=0}^{B-1}A_{y,j_{1}}\sum_{j_{0}=0}^{B^{k}-1}A^{-1}_{j_{0},i_{0}}f'_{x,y,j_{0}}
$$
显然最后一个式子也就是$f'_{x,y}A^{-1}$,记作$f_{x,y}$,根据归纳有$f_{x,y}=d_{x}\otimes e_{y}$,联系其定义即
$$
f_{x,y,i}=\sum_{0\le p_{0},q_{0}<B^{k},p_{0}\oplus q_{0}=i}d_{x,p_{0}}e_{y,q_{0}}
$$
再次代入$c_{i}$中,即
$$
c_{i}=\sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}\sum_{x=0}^{B-1}A_{x,j_{1}}\sum_{y=0}^{B-1}A_{y,j_{1}}\sum_{0\le p_{0},q_{0}<B^{k},p_{0}\oplus q_{0}=i_{0}}d_{x,p_{0}}e_{y,q_{0}}=\sum_{0\le p_{0},q_{0}<B^{k},p_{0}\oplus q_{0}=i_{0}}\sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}\sum_{x=0}^{B-1}A_{x,j_{1}}d_{x,p_{0}}\sum_{y=0}^{B-1}A_{y,j_{1}}e_{y,q_{0}}
$$
对于$p_{0}$和$q_{0}$,构造$u_{i}=d_{i,p_{0}}$以及$v_{i}=e_{i,q_{0}}$,将两者作为$a_{i}$和$b_{i}$进行计算,关于$x$的枚举也就是算出$u'_{j_{1}}$,第二个关于$y$的枚举也就是$v'_{j_{1}}$,两者对应位置相乘,也就是$c'_{j_{1}}$,再乘上前面的$A_{j_{1},i_{1}}^{-1}$,也就是$c_{i_{1}}$
(这里的$c'_{j_{1}}$和$c_{i_{1}}$指代入$u_{i}$和$v_{i}$后用前面的方式计算得到,与之前的$c'_{i}$与$c_{i}$无关)
进一步的,根据归纳有$c_{i_{1}}=(u\otimes v)_{i_{1}}=\sum_{0\le p_{1},q_{1}<B,p_{1}\oplus q_{1}=i_{1}}u_{p_{1}}v_{q_{1}}$,代入即
$$
c_{i}=\sum_{0\le p_{0},q_{0}<B^{k},p_{0}\oplus q_{0}=i_{0}}\sum_{0\le p_{1},q_{1}<B,p_{1}\oplus q_{1}=i_{1}}a_{p_{0}B+p_{1}}b_{q_{0}B+q_{1}}=\sum_{0\le p,q<B^{k+1},p\oplus q=i}a_{p}b_{q}=(a\otimes b)_{i}
$$
(最后就利用到了$\oplus$是按位运算的)
确定矩阵$A$的条件后,我们还要具体的构造出一个矩阵$A$,通过矩阵左上角$B\times B$的部分就可以根据第二个条件唯一确定$A$了,也就是如何构造$n=B$时的矩阵
不难发现这就是一个幂次对$B$取模的加法卷积,联系FFT,记$\omega=\cos\frac{2\pi}{B}+\sin\frac{2\pi}{B}i$,取$A_{i,j}=\omega^{ij}$,逆矩阵将则令$A^{-1}_{i,j}=\frac{\omega^{-(ij)}}{B}$,下面来证明一下:
当我们确定了$A_{i,j}=\omega^{i}$,相当于将$a_{i}$看作一个$B$次多项式,求出其在$\omega^{i}$的值
根据余数定理,当我们知道一个多项式在$x=\omega^{i}$上的值,事实上也就是该多项式模$x-\omega^{i}$后的结果,而求出了所有$\omega^{i}$,那么也就知道了原多项式对$\prod_{i=0}^{B-1}(x-\omega^{i})$取模的结果
再将点值相乘, 也就确定了最终多项式对$\prod_{i=0}^{B-1}(x-\omega^{i})$取模后的结果
注意到$\omega$满足$\omega^{B}=1$且$\forall 1\le i<B,\omega^{i}\ne 1$,则$\prod_{i=0}^{B-1}(x-\omega^{i})=x^{B}-1$,而$x^{n}\equiv x^{n\ mod\ B}(mod\ x^{B}-1)$,因此恰好也就求出了下标对$B$取模的卷积结果
关于$\prod_{i=0}^{B-1}(x-\omega^{i})=x^{B}-1$,由于其可以看作$B$对点值$\{x=\omega^{i},y=0\}$,且每一个$x$各不相同,由此恰好可以确定一个$B$次多项式,且$x^{B}-1$恰好满足此性质,因此即相等
然而,本题需要对$p$取模,如果用实数表示的误差无法解决,代入$B=3$,可得$\omega=\frac{\sqrt{3}}{2}i-\frac{1}{2}$
考虑题中关于$p$的限制:不存在正整数$x,y$满足$\frac{1}{x}+\frac{1}{y}=\frac{3}{p}$
若$p\equiv 0(mod\ 3)$,取$x=y=\frac{2p}{3}$即满足条件,因此$\frac{1}{3}$在模$p$意义下有意义,但并不能保证3在模$p$意义下有二次剩余,换言之$\sqrt{3}$不一定有意义($\frac{1}{2}$有意义也可以证,但不关心)
为了使系数为整数,考虑将$\omega$作为一个整体,显然每一个虚数与$a+b\omega$一一对应,同时根据$\omega^{3}=1$以及$\omega^{2}=-\omega-1$(这两个条件也就是根据$\prod_{i=0}^{B-1}(x-\omega^{i})=x^{B}-1$得到),当产生高次项也都可以降为一次来表示
最终答案必然是整数,因此$b=0$,即为$a$
另外,计算矩阵乘法可以利用条件2来做,以由$a$计算$a'$为例子,即
$$
a'_{i}=\sum_{j=0}^{B^{k}-1}a_{j}A_{j,i}=\sum_{j_{1}=0}^{B-1}A_{j_{1},i_{1}}\sum_{j_{0}=0}^{B^{k-1}-1}a_{j_{0}B+j_{1}}A_{j_{0},i_{0}}
$$
对于每一个$j_{1}$,后面也就是一个$B^{k-1}$规模的子问题(矩阵乘法),递归后将其求出后,再代入该式即可,这样单次矩乘时间复杂度为$o(kB^{k+1})$
由于要快速幂,总复杂度为$o((m+\log_{2}t)3^{m})$可以通过
补充一下与本题无关的一些事情:
$B=2$时这个问题就是FWT,且做法与其完全相同
当$B$足够大时即为FFT,但时间复杂度会退化为$o(n^{2})$
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 600005 4 #define M 13 5 int n,m,t,mod,S[N],base[N][M],b[M][M]; 6 struct Complex{ 7 int a,b; 8 Complex(){ 9 a=b=0; 10 } 11 Complex(int x){ 12 a=x,b=0; 13 } 14 Complex(int x,int y){ 15 a=x,b=y; 16 } 17 Complex operator + (const Complex &k)const{ 18 return Complex((a+k.a)%mod,(b+k.b)%mod); 19 } 20 Complex operator * (const Complex &k)const{ 21 int s=mod-1LL*b*k.b%mod; 22 return Complex((1LL*a*k.a+s)%mod,(1LL*a*k.b+1LL*b*k.a+s)%mod); 23 } 24 }inv3,A[3][3],invA[3][3],a[N],ans[N]; 25 Complex pow(Complex n,int m){ 26 Complex s=n,ans=Complex(1); 27 while (m){ 28 if (m&1)ans=ans*s; 29 s=s*s; 30 m>>=1; 31 } 32 return ans; 33 } 34 void DFT(Complex *a){ 35 Complex aa[3]; 36 for(int i=0,s=1;i<m;i++,s*=3) 37 for(int j=0;j<n;j++) 38 if (!base[j][i]){ 39 for(int k=0;k<3;k++)aa[k]=Complex(); 40 for(int k=0;k<3;k++) 41 for(int l=0;l<3;l++)aa[k]=aa[k]+a[j+l*s]*A[l][k]; 42 for(int k=0;k<3;k++)a[j+k*s]=aa[k]; 43 } 44 } 45 void IDFT(Complex *a){ 46 Complex aa[3]; 47 for(int i=0,s=1;i<m;i++,s*=3) 48 for(int j=0;j<n;j++) 49 if (!base[j][i]){ 50 for(int k=0;k<3;k++)aa[k]=Complex(); 51 for(int k=0;k<3;k++) 52 for(int l=0;l<3;l++)aa[k]=aa[k]+a[j+l*s]*invA[l][k]; 53 for(int k=0;k<3;k++)a[j+k*s]=aa[k]; 54 } 55 } 56 int main(){ 57 scanf("%d%d%d",&m,&t,&mod); 58 n=1; 59 for(int i=1;i<=m;i++)n*=3; 60 for(int i=0;i<n;i++){ 61 S[i]=S[i/3]+(i%3==1); 62 base[i][0]=i%3; 63 for(int j=1;j<m;j++)base[i][j]=base[i/3][j-1]; 64 } 65 int k=mod,phi=mod; 66 for(int i=2;i*i<=k;i++) 67 if (k%i==0){ 68 phi=phi/i*(i-1); 69 while (k%i==0)k/=i; 70 } 71 if (k>1)phi=phi/k*(k-1); 72 inv3=pow(Complex(3),phi-1); 73 for(int i=0;i<3;i++) 74 for(int j=0;j<3;j++){ 75 A[i][j]=pow(Complex(0,1),i*j); 76 invA[i][j]=pow(Complex(0,1),6-i*j)*inv3; 77 } 78 for(int i=0;i<n;i++){ 79 scanf("%d",&k); 80 a[i]=Complex(k); 81 } 82 for(int i=0;i<=m;i++) 83 for(int j=0;j<=m-i;j++)scanf("%d",&b[i][j]); 84 for(int i=0;i<n;i++){ 85 k=0; 86 for(int j=m-1;j>=0;j--) 87 if (!base[i][j])k=k*3; 88 else k=k*3+3-base[i][j]; 89 ans[i]=Complex(b[S[i]][S[k]]); 90 } 91 DFT(a); 92 DFT(ans); 93 for(int i=0;i<n;i++)ans[i]=a[i]*pow(ans[i],t); 94 IDFT(ans); 95 for(int i=0;i<n;i++)printf("%d\n",ans[i].a); 96 }