Loading

$k$ 进制 FWT 小记

上文：FWT 快速沃尔什变换

概念

\(k\) 进制的 DWT 本质上是二进制操作的一个扩展操作。

原来的位运算也推广到了 \(k\) 进制上。

\(\text{or}\) 卷积

  • 定义：两个数 \(x_{1...k},y_{1...k}\) 的 \(k\) 进制或运算为：\(z_i=\max(x_i,y_i)\)。

换言之，在每个维度上取 \(\max\)。

上文中，我们提到了 FWT 的本质，变换矩阵需要满足 \(w(i,j)w(i,k)=w(i,j\oplus k)\)。

对于或卷积，考虑某一位，构造一个矩阵满足 \(w(i,j)w(i,k)=w(i,\max(j,k))\)。

钦定 \(j<k\)，那么 \(w(i,j)w(i,k)=w(i,k)\)。

当 \(w(i,k)=0\) 时，\(w(i,j)\) 随意。

当 \(w(i,k)\not =0\) 时，\(w(i,j)=1\)。

钦定整个矩阵只有 \(0\) 和 \(1\)，那么 \(w(i,j)=0,w(i,k)=1\) 是不合法的。换言之，每一行中所有 \(1\) 在所有 \(0\) 前面。

随便构造一个有逆的，这里采用这种：

\[\begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 \\ 1 & 1 & 1 & 1 \end{bmatrix} \]

其实就是 \(k\) 维前缀和。

\(\text{and}\) 卷积也很类似，就是高位后缀和。

\(\text{xor}\) 卷积

  • 定义：\(z_i=(x_i+y_i)\bmod k\)

即两个数的 \(k\) 进制不进位加法。

此时有 \(w(x,y)w(x,z)=w(x,(y+z)\bmod k)\)。

我们发现把单位根带入就行了：\(w(x,y)=\omega_k^y\)

带入上式 \(\omega_k^y\cdot \omega_k^z = \omega_k^{(y+z)\bmod k}\)

但是每行所有数都相同，没有逆，考虑带入范德蒙德矩阵：

\[\begin{bmatrix} 1 & 1 & 1 & ... & 1 \\ 1 & \omega_k^1 & \omega_k^2 & ... & \omega_k^{k-1} \\ 1 & \omega_k^2 & \omega_k^4 & ... & \omega_k^{2(k-1)} \\ ... & ... & ... & ... & ...\\ 1 & \omega_k^{k-1} & \omega_k^{2(k-1)} & ... & \omega_k^{(k-1)(k-1)} \\ \end{bmatrix} \]

发现仍然成立，此时 \(w(x,y)=\omega_k^{xy}\)。而且他的逆是已知的：

\[ \frac 1k \begin{bmatrix} 1 & 1 & 1 & ... & 1 \\ 1 & \omega_k^{-1} & \omega_k^{-2} & ... & \omega_k^{-(k-1)} \\ 1 & \omega_k^{-2} & \omega_k^{-4} & ... & \omega_k^{-2(k-1)} \\ ... & ... & ... & ... & ...\\ 1 & \omega_k^{-(k-1)} & \omega_k^{-2(k-1)} & ... & \omega_k^{-(k-1)(k-1)} \\ \end{bmatrix} \]

但是又有一个问题：\(k\) 不一定是 \(2\) 的幂，模意义下的单位根不一定存在。

考虑扩域，用 \(x\) 代替 \(\omega_k\)，那么每个数可以表示成 \(x\) 的 \(k-1\) 次多项式。我们把多项式表示为 \(x^k-1\) 模意义下，可以沿用单位根的“指数取模”的性质（即 \(x^k\equiv 1\pmod {x^k-1}\)）

但是这样是错误的，这个模数可能导致存在零因子，即可能存在一个多项式的逆元无法表示。

本质原因是：\(x^k-1\) 可以因式分解，等于 \((x-1)(1+x+x^2+x^3+...+x^{k-1})\)，所以 \(x-1\) 不存在逆元。而在范德蒙德矩阵和他的逆矩阵的严格互逆证明需要用到类似于单位根反演的东西，但是在单位根反演的等比数列求和式子中，\(x-1\) 作为除数不一定存在逆元。相关式子详见：单位根反演小记

所以我们需要一个不能因式分解的模数，这里引入分圆多项式。

分圆多项式

定义分圆多项式为 \(\Phi_k(x)=\prod\limits_{1\le i\le n,\space i\bot k} (x-\omega_k^i)\)

分圆多项式的非复数域计算：

\[\Phi_k(x)=\sum_{d|k} (x^d-1)^{\mu(\frac kd)} \]

基本性质：

  • \(x^k-1=\prod_{i=1}^n (x-\omega_k^i)\)

  • 类似于欧拉反演，我们有 \(x^k-1=\prod_{d|k} \Phi_d(x)\)

这里我们不加证明地给出两条重要性质：

  • \(\Phi_k(x)\) 不能因式分解

  • \(\Phi_k(x)|(x^k-1)\)

所以我们可以直接使用分圆多项式作为模数，这样正确定就得到了保障。

但是每次做一次运算后都要取一次模，时间复杂度不优秀。

注意到第二条重要性质 \(\Phi_k(x)|(x^k-1)\)，说明我们其实可以把多项式先对 \(x^k-1\) 取模，最后再对 \(\Phi_k(x)\) 取模。

注意此时的正确性仍然是有保证的。

整个过程中枚举每一位 \(O(m)\)，枚举所有下标 \(O(k^{m-1})\)，线性变换的计算（含多项式平移）\(O(k^3)\)，总时间复杂度 \(O(mk^{m+2})\)。

例题

题意：有 \(n\) 个数 \(a_{1...n}\)，维护一个值 \(s\)。每次选择一个数，不进位地累加进 \(s\) 中，可以重复选数。求对于 \(i\in [0,n-1]\)，最终 \(s=i\) 的方案数。

\(1\le n\le 10^5,\space 0\le a_i<10^5\)

题目意思就是求 \((1+\sum\limits_{i=1}^nx^{a_i})^n\)，是 FWT 的板子，这里是 \(10\) 进制。

这里直接给出代码。

点击查看代码 
#include<bits/stdc++.h>
#define ll long long
#define ull unsigned ll
#define pir pair<ll,ll>
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
using namespace std;
const ll maxn=1e5+10, k=10;
ll n,a[maxn];
struct node{
	ull dat[10]={};
	ull operator[] (ll x) const {return dat[x];}
	ull &operator[] (ll x) {return dat[x];}
	void operator+= (const node &t) {for(ll i=0;i<k;i++) dat[i]+=t[i];} 
	void operator*= (const ull x) {for(ll i=0;i<k;i++) dat[i]*=x;} 
	node operator* (const node &y){
		static ull a[20];
		for(ll i=0;i<(k<<1);i++) a[i]=0;
		for(ll i=0;i<k;i++)
			for(ll j=0;j<k;j++)
				a[i+j]+=dat[i]*y[j];
		node res;
		for(ll i=0;i<k;i++) res[i]=a[i]+a[i+k];
		return res;
	}
	node operator+ (const node &y){
		node res;
		for(ll i=0;i<k;i++) res[i]=dat[i]+y[i];
		return res;
	}
	node mov(ll x){
		x=(x%k+k)%k; node res;
		for(ll i=0;i<k;i++) res[(i+x)%k]=dat[i];
		return res;
	}
}b[maxn],ans[maxn],Mod;
ull pw[maxn];
void fwt(node t[],ll op){
	for(ll i=0;i<5;i++)
		for(ll j=0;j<pw[5];j+=pw[i+1])
			for(ll x=0;x<pw[i];x++){
				node ret[10];
				for(ll p=0;p<k;p++)
					for(ll q=0;q<k;q++)
						ret[q]+=t[j+pw[i]*p+x].mov(p*q*op);
				for(ll p=0;p<k;p++) t[j+pw[i]*p+x]=ret[p];
			}
}
ull inv;
ull getans(node &t){
	for(ll i=k-1;i>=4;i--){
		ll w=t[i];
		for(ll j=0;j<=4;j++) t[i-j]-=w*Mod[4-j];
	}
	ull x=t[0];
	x>>=5; return x*inv%(1ull<<58);
}
int main(){
	pw[0]=1, pw[1]=10, pw[2]=100, pw[3]=1000, pw[4]=10000, pw[5]=100000;
	scanf("%lld",&n);
	for(ll i=1;i<=n;i++){
		scanf("%lld",a+i);
		++b[a[i]][0];
	}
	fwt(b,1);
	for(ll i=0;i<pw[5];i++){
		ll r=n; ans[i][0]=1;
		while(r){
			if(r&1) ans[i]=ans[i]*b[i];
			b[i]=b[i]*b[i]; r>>=1;
		}
	}
	fwt(ans,-1);
	Mod[0]=Mod[2]=Mod[4]=1, Mod[1]=Mod[3]=-1;
	inv=14757395258967641293ull;
	inv=inv*inv*inv*inv*inv;
	for(ll i=0;i<n;i++){
		printf("%lld\n",getans(ans[i]));
	}
	return 0;
}

题意：求 \(m\) 位 \(k\) 进制下 \(\prod\limits_{i=1}^n (1+x^{a_i})\)。

\(1\le n\le 10^6,\space 1\le m\le 7,\space 5\le k\le 6,\space m+k\le 12\)

此时 \(n\) 很大，我们不能把每个数都 FWT 一下。

考虑直接计算 \(\prod\limits_{i=1}^n fwt(1+x^{a_i})\) 即所有集合幂级数 FWT 后的结果。

考虑对于他的第 \(x\) 位，\(a_{1...n}\) 中每个数的贡献都是形如 \((1+\omega_k^i)\) 的形式，其中 \(i\in [0,k-1]\)。

那么我们只需要对于 \(i\in[0,k-1]\)，分别计算有多少次贡献。

然后我们会发现，其实这个计算过程和 FWT 很像。具体的，我们把计数数组 \(cnt\) 拿来 FWT，然后 \(fwt(cnt)_x\) 上的多项式各项系数便恰好是每个 \((1+\omega_k^i)\) 的贡献次数，这个可以自行脑补。

接下来我们要计算 \((1+\omega_k^i)\) 的幂，这个可以光速幂。

时间复杂度 \(O(n+mk^{m+2})\)。

点击查看代码 
#include<bits/stdc++.h>
#define ll long long
#define ull unsigned ll
#define pir pair<ll,ll>
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
using namespace std;
const ll maxn=1e6+10, mod=998244353;
ll n,k,m,a[maxn];
void ad(ll &x,const ll &y) {x=(x+y>=mod? x+y-mod:x+y);}
ll power(ll a,ll b=mod-2){
	ll s=1;
	while(b){
		if(b&1) s=s*a%mod;
		a=a*a%mod; b>>=1; 
	} return s;
}
struct node{
	ll dat[6]={};
	ll operator[] (ll x) const {return dat[x];}
	ll &operator[] (ll x) {return dat[x];}
	void operator+= (const node &t) {for(ll i=0;i<k;i++) ad(dat[i],t[i]);}
	void operator*= (const ll x) {for(ll i=0;i<k;i++) dat[i]=dat[i]*x%mod;}
	node operator* (const node &t){
		static ll a[11];
		for(ll i=0;i<(k<<1)-1;i++) a[i]=0;
		for(ll i=0;i<k;i++)
			for(ll j=0;j<k;j++)
				ad(a[i+j],dat[i]*t[j]%mod);
		node res;
		for(ll i=0;i<k;i++) res[i]=(a[i]+a[i+k])%mod;
		return res;
	}
	node mov(ll x){
		x=(x%k+k)%k;
		node res;
		for(ll i=0;i<k;i++) res[(i+x)%k]=dat[i];
		return res;
	}
}cnt[maxn],pw1[7][1010],pw2[7][1010],ans[maxn],Mod;
ll pw[maxn], sq;
void fwt(node t[],ll op){
	for(ll i=0;i<m;i++)
		for(ll j=0;j<pw[m];j+=pw[i+1])
			for(ll x=0;x<pw[i];x++){
				node ret[6];
				for(ll p=0;p<k;p++)
					for(ll q=0;q<k;q++)
						ret[q]+=t[j+p*pw[i]+x].mov(p*q*op);
				for(ll p=0;p<k;p++) t[j+p*pw[i]+x]=ret[p];
			}
}
node qry(ll x,ll n){
	return pw1[x][n/sq]*pw2[x][n%sq];
} ll deg, inv;
ll getans(node &t){
	for(ll i=k-1;i>=deg;i--){
		ll w=t[i];
		for(ll j=0;j<=deg;j++) ad(t[i-j],mod-w*Mod[deg-j]%mod);
	}
	ll x=t[0];
	return x*inv%mod;
}
int main(){
	scanf("%lld%lld%lld",&n,&k,&m);
	pw[0]=1;
	for(ll i=1;i<=m;i++) pw[i]=pw[i-1]*k;
	sq=sqrt(n)+1;
	for(ll i=0;i<k;i++){
		node g; ++g[0], ++g[i];
		pw2[i][0][0]=pw1[i][0][0]=1;
		for(ll j=1;j<=sq;j++) pw2[i][j]=pw2[i][j-1]*g;
		for(ll j=1;j<=sq;j++) pw1[i][j]=pw1[i][j-1]*pw2[i][sq];
	}
	for(ll i=1;i<=n;i++){
		ll x;
		scanf("%lld",&x);
		for(ll j=0,v=1;j<m;j++,v*=10) a[i]+=x/v%10*pw[j];
		++cnt[a[i]][0];
	}
	fwt(cnt,1);
	for(ll i=0;i<pw[m];i++){
		ans[i][0]=1;
		for(ll j=0;j<k;j++){
			node tmp=qry(j,cnt[i][j]);
			ans[i]=ans[i]*tmp;
		}
	}
	fwt(ans,-1);
	if(k==5) Mod[0]=Mod[1]=Mod[2]=Mod[3]=Mod[4]=1, deg=4;
	else Mod[0]=Mod[2]=1, Mod[1]=mod-1, deg=2;
	inv=power(pw[m]);
	for(ll i=0;i<pw[m];i++){
		printf("%lld\n",getans(ans[i]));
	}
	return 0;
}

附：\(k=1\sim 6\) 的分圆多项式表

\(k\) \(\Phi_k(x)\)
\(1\) \(x-1\)
\(2\) \(x+1\)
\(3\) \(x^2+x+1\)
\(4\) \(x^2+1\)
\(5\) \(x^4+x^3+x^2+x+1\)
\(6\) \(x^2-x+1\)
posted @ 2024-03-19 11:49  Lgx_Q  阅读(24)  评论(0)  编辑  收藏  举报