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

最新推荐文章于 2024-11-13 23:40:34 发布

v5zsq

最新推荐文章于 2024-11-13 23:40:34 发布

阅读量284

点赞数

CC 4.0 BY-SA版权

分类专栏： Newcoder 递归 FWT

本文链接：https://blog.youkuaiyun.com/V5ZSQ/article/details/82774030

Newcoder 同时被 3 个专栏收录

186 篇文章

订阅专栏

递归

18 篇文章

订阅专栏

FWT

7 篇文章

订阅专栏

本文介绍了一种解决特定类型的矩阵同余方程的方法，该方程涉及到异或运算。通过利用快速沃尔什-哈达玛变换(FWT)简化求解过程，实现了高效的算法设计。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Description

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

Input

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

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

Output

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

Sample Input

4
1 10 100 1000
1234 2143 3412 4321

Sample Output

4
3
2
1

Solution1

$Ai,j=ai⊙jA_{i,j}=a_{i\odot j}$ ，原矩阵方程等价于方程组 $∑j=1nai⊙jxj=bi\sum\limits_{j=1}^n a_{i\odot j}x_j=b_i$ ，也即 $∑j⊙k=iajxk=bi\sum\limits_{j\odot k=i}a_jx_k=b_i$ ，这意味着 $a⊗x=ba\otimes x=b$ ，其中 $⊙\odot$ 表示异或运算， $⊗\otimes$ 表示异或卷积，那么有 $FWT(a)⋅FWT(x)=FWT(b)FWT(a)\cdot FWT(x)=FWT(b)$ ，其中 $⋅\cdot$ 表示点乘，故对 $a, b$ 序列做 $F W T$ 后得到 $bi⋅ai−1b_i\cdot a_i^{-1}$ 序列再 $U F W T$ 即得到 $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⊙ji\odot j$ 矩阵，由于 $n$ 为 $2$ 的整数次幂，简单分析可知 $A$ 矩阵可以分块用两个 $n2⋅n2\frac{n}{2}\cdot \frac{n}{2}$ 的矩阵表示
$\left[\begin{matrix} A^1\ A^2\\ A^2\ A^1 \end{matrix}\right]$
同时将 $x, b$ 分成前后两个等长的部分， $x=[x^1,x^2]^T,b=[b^1,b^2]^T$ ，那么原方程可以表示为
$A^1x^1+A^2x^2=b^1,A^2x^1+A^1x^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)$
令 $ai1=ai+ai+n2,xi1=xi+xi+n2,bi1=bi+bi+n2,ai2=ai−ai+n2,xi2=xi−xi+n2,bi2=bi−bi−n2a^1_{i}=a_i+a_{i+\frac{n}{2}},x^1_i=x_i+x_{i+\frac{n}{2}},b^1_i=b_i+b_{i+\frac{n}{2}},a^2_{i}=a_i-a_{i+\frac{n}{2}},x^2_i=x_i-x_{i+\frac{n}{2}},b^2_i=b_i-b_{i-\frac{n}{2}}$ ，

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

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;
}