Newcoder 147 A.Circulant Matrix（FWT+递归）

Description

给出两个长度为$n$的序列$a_0,...,a_{n-1}$和$b_0,...,b_{n-1}$，以此定义矩阵$A_{i,j}=a_{i\odot j}$，其中$\odot$为异或运算，求解矩阵同余方程$Ax=b(modp)$

Input

第一行一整数$n$表示序列长度，之后输入两个长度为$n$的序列$a_i,b_i$

$(1\le n\le 222144,p=10^9+7,0\le a_i,b_i<p)$

Output

保证有唯一解，输出$x_0,...,x_{n-1}$，需保证$0\le x_i<p$

Sample Input

4
1 10 100 1000
1234 2143 3412 4321

Sample Output

4
3
2
1

Solution1

$A_{i,j}=a_{i\odot j}$，原矩阵方程等价于方程组$\sum _{j=1}^n{a}_{i\odot j}{x}_j=b_i$，也即$\sum _{j\odot k=i}{a}_j{x}_k=b_i$，这意味着$a\otimes x=b$，其中$\odot$表示异或运算，$\otimes$表示异或卷积，那么有$FWT(a)\cdot FWT(x)=FWT(b)$，其中$\cdot$表示点乘，故对$a,b$序列做$FWT$后得到$b_i\cdot a_i^{-1}$序列再$UFWT$即得到$x$序列

Code1

#include<cstdio>
using namespace std;
typedef long long ll;
#define maxn 262144+5
#define mod 1000000007
#define inv2 500000004
int mul(int x,int y)
{
	ll z=1ll*x*y;
	return z-z/mod*mod;
}
int add(int x,int y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x;
}
int Pow(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1)ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}
int n,a[maxn],b[maxn],x[maxn];
void FWT(int *a,int n)  
{  
    for(int d=1;d<n;d<<=1)  
        for(int i=0;i<n;i+=(d<<1))  
            for(int j=0;j<d;j++)  
	            {  
	                int x=a[i+j],y=a[i+j+d];  
	                a[i+j]=add(x,y),a[i+j+d]=add(x,mod-y);   
	            } 
}   
void UFWT(int *a,int n)  
{  
    for(int d=1;d<n;d<<=1)  
        for(int i=0;i<n;i+=(d<<1))  
            for(int j=0;j<d;j++)  
            {  
                int x=a[i+j],y=a[i+j+d];  
                a[i+j]=mul(inv2,add(x,y)),a[i+j+d]=mul(inv2,add(x,mod-y));   
            }  
} 
int main()
{
	scanf("%d",&n);
	for(int i=0;i<n;i++)scanf("%d",&a[i]);
	for(int i=0;i<n;i++)scanf("%d",&b[i]);
	FWT(b,n),FWT(a,n);
	for(int i=0;i<n;i++)x[i]=mul(b[i],Pow(a[i],mod-2));
	UFWT(x,n);
	for(int i=0;i<n;i++)printf("%d\n",x[i]);
	return 0;
}

Solution2

考虑$i\odot j$矩阵，由于$n$为$2$的整数次幂，简单分析可知$A$矩阵可以分块用两个$\frac{n}{2}\cdot \frac{n}{2}$的矩阵表示
$$A=\left[\begin{array}{c}{A}^1{A}^2\\ A^2{A}^1\end{array}\right]$$
同时将$x,b$分成前后两个等长的部分，$x=[x^1,x^2{]}^T,b=[b^1,b^2{]}^T$，那么原方程可以表示为
$$A^1{x}^1+A^2{x}^2=b^1,A^2{x}^1+A^1{x}^2=b^2$$
两式相加/做差得
$$(A^1+A^2)(x^1+x^2)=b^1+b^2,(A^1-A^2)(x^1-x^2)=(b^1-b^2)$$
令$a_i^1=a_i+a_{i+\frac{n}{2}},x_i^1=x_i+x_{i+\frac{n}{2}},b_i^1=b_i+b_{i+\frac{n}{2}},a_i^2=a_i-a_{i+\frac{n}{2}},x_i^2=x_i-x_{i+\frac{n}{2}},b_i^2=b_i-b_{i-\frac{n}{2}}$，

那么简单分析可知$a^1$序列可以用同样的方式生成矩阵$B^1=A^1+A^2$，$a^2$也可以生成$B^2=A^1-A^2$，故求解原来的$n\times n$规模的矩阵方程问题转化为两个$\frac{n}{2}\times \frac{n}{2}$规模的矩阵方程问题
$$B^1{x}^1=b^1,B^2{x}^2=b^2$$
递归下去，以一般方程作为递归终点即可，时间复杂度$O(nlogn)$，注意到这种解法其实就是模拟了$FWT$的过程

Code2

#include<cstdio>
using namespace std;
typedef long long ll;
#define maxn 262144+5
#define mod 1000000007
#define inv2 500000004
int n,a[19][maxn],b[19][maxn],x[19][maxn],y[19][maxn];
int mul(int x,int y)
{
	ll z=1ll*x*y;
	return z-z/mod*mod;
}
int add(int x,int y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x;
}
int Pow(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1)ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}
void Solve(int n,int m)
{
	if(m==0)
	{
		x[0][0]=mul(b[0][0],Pow(a[0][0],mod-2));
		return ;
	}
	for(int i=0;i<n/2;i++)
	{
		a[m-1][i]=add(a[m][i],a[m][i+n/2]);
		b[m-1][i]=add(b[m][i],b[m][i+n/2]);
	}
	Solve(n/2,m-1);
	for(int i=0;i<n/2;i++)y[m][i]=x[m-1][i];//x[i]+x[i+n/2]
	for(int i=0;i<n/2;i++)
	{
		a[m-1][i]=add(a[m][i],mod-a[m][i+n/2]);
		b[m-1][i]=add(b[m][i],mod-b[m][i+n/2]);
	}
	Solve(n/2,m-1);
	for(int i=0;i<n/2;i++)y[m][i+n/2]=x[m-1][i];//x[i]-x[i+n/2]
	for(int i=0;i<n/2;i++)
	{
		x[m][i]=mul(inv2,add(y[m][i],y[m][i+n/2]));
		x[m][i+n/2]=mul(inv2,add(y[m][i],mod-y[m][i+n/2]));
	}
}
int main()
{
	scanf("%d",&n);
	int m=0;
	while((1<<m)<n)m++;
	for(int i=0;i<n;i++)scanf("%d",&a[m][i]);
	for(int i=0;i<n;i++)scanf("%d",&b[m][i]);
	Solve(n,m);
	for(int i=0;i<n;i++)printf("%d\n",x[m][i]);
	return 0;
}