P6097 题解
Doqe
·
·
题解
前置文章。其包含了这里的一些概念。
原始算法
在这题里我们希望：
c^0=W(a^0,b^0),c^1=W(a^0,b^1)+W(a^1,b^0)
好消息是总项数只有 3，直接做也是 O(3^n) 的，坏消息是难以用前置文章中的通用思路构造出小于 3 的项数。
这个形式形如一个真正的卷积，我们搓出一个形式元 x，有：
W(a^0+a^1x,b^0+b^1x)=W(a^0,b^0)+(W(a^0,b^1)+W(a^1,b^0))x+W(a^1,b^1)x^2
也就是 c^0+c^1x+W(a^1,b^1)x^2。
那就好办了。
$$d^0=W(a^0,a^1),d^1=W(a^0+a^1x,b^0+b^1x)$$
$$c^0=d^0,c'^1=d^1-d^0$$
这里的 $c'^1$ 是多项式，真正答案在于最低有效位，可以推出 $i$ 的最低有效位的指数是 $\operatorname{popcount}(i)$。
分析复杂度，在 “FMT” 和 “IFMT” 过程之中，所有操作是加法和乘 $x$，均为线性；其余操作是执行 $2^n$ 次的卷积。总复杂度 $O(2^nn^2)$。
（PS：如果你让 $x=\epsilon$，足够小的话，最后的结果形如 $Ax^k+\sum_{i\ge k+1}Bx^{i}$，剩余部分完全能够被区分出来。只不过这要求 $x=O(nV)$，也就是需要上高精度，导致和一般的子集卷积复杂度相同。）
（PS：这部分和一般不使用 $x$ 的算法，在张量分析下可以使用普通张量秩和边界秩分别解释，即 $\sum a_i\circ b_i\circ c_i=X$ 还是 $\sum a_i\circ b_i\circ c_i=Xx^k+Yo(x^k)$（如果不知道这表示什么东西的话可以忽略），其天然对应占位多项式。）
此时的代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const unsigned P=1e9+9,N=1<<20|3;
unsigned a[N][21],b[N][21];
int n;
inline void ADD(unsigned&A,unsigned B){A=(A+B>=P?A+B-P:A+B);}
void _FMT(unsigned a[][21])
{
for(int i=0;i<n;++i)
for(int j=0;j<(1<<n);++j)if(!(j>>i&1))
{
unsigned *A=a[j|(1<<i)],*B=a[j];
memmove(A+1,A,80);A[0]=B[0];
for(int i=1;i<=n;++i)ADD(A[i],B[i]);
}
}
void _IFMT(unsigned a[][21])
{
for(int i=0;i<n;++i)
for(int j=0;j<(1<<n);++j)if(!(j>>i&1))
{
unsigned *A=a[j|(1<<i)],*B=a[j];
for(int i=0;i<=n;++i)ADD(A[i],P-B[i]);
}
}
unsigned long long tmp[21];
int main()
{
cin.tie(0)->sync_with_stdio(0);
cin>>n;
for(int i=0;i<(1<<n);++i)cin>>a[i][0];
for(int i=0;i<(1<<n);++i)cin>>b[i][0];
_FMT(a),_FMT(b);
for(int i=0;i<(1<<n);++i)
{
memset(tmp,0,sizeof tmp);
for(int j=0;j<=n;++j)
for(int k=0;j+k<=n;++k)
(tmp[j+k]+=1ull*a[i][j]*b[i][k]))%=P;
for(int j=0;j<=n;++j)
a[i][j]=tmp[j]%P;
}
_IFMT(a);
for(int i=0;i<(1<<n);++i)cout<<a[i][__builtin_popcount(i)]<<" ";cout<<endl;
}
```
### 改进算法：位置对齐
如果初始输入直接放入 $x^{\operatorname{popcount}(i)}$ 的系数，“FMT”中的 `memmove` 就可以省去。更改部分的代码：
```cpp
void _FMT(unsigned a[][21])
{
for(int i=0;i<n;++i)
for(int j=0;j<(1<<n);++j)if(!(j>>i&1))
{
unsigned *A=a[j|(1<<i)],*B=a[j];
for(int i=0;i<=n;++i)ADD(A[i],B[i]);
}
}
//do other thing
int main()
{
cin.tie(0)->sync_with_stdio(0);
cin>>n;
for(int i=0;i<(1<<n);++i)cin>>a[i][__builtin_popcount(i)];
for(int i=0;i<(1<<n);++i)cin>>b[i][__builtin_popcount(i)];
//do other thing
}
```
### 进一步：线性空间
如果我们强制线性空间，那么多项式 $f(x)$ 是无法存储的。
但是我们可以存储 $f(B)$，其中 $B$ 是一个初始固定的值。接下来可以尝试使用拉格朗日插值求系数。
但是不需要求系数这么麻烦。我们知道对于 $S$，记 $c=\operatorname{popcount}(S)$，则 $f_S(x)=ans_Sx^c+\sum_{j>c}z_{S,j}x^j$。
我们可以对 $g_S(x)=\frac{f_S(x)}{x^c}$ 求 $g_S(0)$，即使 $\frac{f_S(0)}{0^c}$ 看起来似乎有问题，但如果中间看成多项式操作就完全能够解释。（这里是因为足够好的连续函数除法会导致一些点上的间断，但是结果仍然可以通过极限弥补成“连续”的，即 $g_S(0)=\lim_{B\to 0}g_S(B)$）
于是我们可以直接预处理拉格朗日插值的系数 $c_j$，然后对于所有 $S$ 有 $g_S(0)=\sum_j c_jg_S(j)$，选取 $j\ge O(\log n)$，这样 $j$ 超过了多项式最大度数（似乎是 $\log n+O(1)$？），最后枚举 $1\le B\le j$，对于所有 $S$ 求 $g_S(B)$ 并做出对应贡献。
时间复杂度不变，空间从 $O(2^nn)$ 转到 $O(2^n)$。
```cpp
void _FMT(unsigned long long a[])
{
for(int i=0;i<n;++i)
for(int j=0;j<(1<<n);++j)if(!(j>>i&1))
ADD(a[j|(1<<i)],a[j]);
}
void _IFMT(unsigned long long a[])
{
for(int i=0;i<n;++i)
for(int j=0;j<(1<<n);++j)if(!(j>>i&1))
ADD(a[j|(1<<i)],P-a[j]);
}
const int M=50;
long long K[M+3];
long long pw[M+3],ipw[M+3];
long long inv[M+3];
int main()
{
cin.tie(0)->sync_with_stdio(0);
inv[1]=1;
for(int i=2;i<=M;++i)inv[i]=(P-1ll*P/i)*inv[P%i]%P;
cin>>n;
for(int i=0;i<(1<<n);++i)cin>>a[i];
for(int i=0;i<(1<<n);++i)cin>>b[i];
for(int i=1;i<=M;++i)
{
long long A=1,B=1;
for(int j=1;j<=M;++j)if(i!=j)
A=A*(P-j)%P;
for(int j=1;j<=M;++j)if(i!=j)
if(i<j) B=B*(P-inv[j-i])%P;
else B=B*(inv[i-j])%P;
K[i]=A*B%P;
}
for(int _=1;_<=M;++_)
{
pw[0]=ipw[0]=1;
for(int j=1;j<=M;++j)pw[j]=pw[j-1]*_%P,ipw[j]=ipw[j-1]*inv[_]%P;
for(int i=0;i<(1<<n);++i)
A[i]=a[i]*pw[__builtin_popcount(i)]%P,
B[i]=b[i]*pw[__builtin_popcount(i)]%P;
_FMT(A),_FMT(B);
for(int i=0;i<(1<<n);++i)A[i]=1ll*A[i]*B[i]%P;
_IFMT(A);
for(int i=0;i<(1<<n);++i)(ans[i]+=K[_]*(A[i]*ipw[__builtin_popcount(i)]%P))%=P;
}
for(int i=0;i<(1<<n);++i)cout<<ans[i]<<" ";cout<<endl;
}
```