转置原理多点求值 和 电子科大校赛2022 简单算术
转置原理多点求值
rushcheyo 《转置原理及其应用》 https://jkloverdcoi.github.io/2020/08/04/转置原理及其应用/slide.pdf
https://blog.csdn.net/weixin_43960287/article/details/122691212OIer学习的功利性很强情有可原，但是我要的是严谨的证明。
转置原理
若矩阵\(A\)可逆，则\(A\)可表示为初等矩阵的积。
令\(A=P_1P_2\cdots P_n\)，则\(A^T=P_n^T\cdots P_2^TP_1^T\)。而三种初等矩阵的转置分别为
-
\(E(i(k))^T=E(i(k))\)不变。
-
\(E(i,j)^T=E(i,j)\)不变。
-
\(E(i,j(k))^T=E(j,i(k))\)相反。
多项式卷积的转置
我们可以把\(n\)次多项式的卷积\(H(x)=F(x)\times G(x)\bmod x^{n+1}\)看成由多项式系数构成的列向量\(F\)经过等价变换\(M_G\)得到的新的由多项式系数构成的列向量\(H\)。即
注意\(M_G\)只写了\(n+1\)行。
模拟下三角矩阵通过行变换转化成单位矩阵的过程手动拆分一下\(M_G\)，
考察一下\(M_G^T\)的形式
所以即使\(M_G\)不可逆（\(g_0=0\)），我们的\(M_G\)仍然是满足使用转置定理的条件的。而\(M_G^T\)本质上就是将\(M_G\)行列交换了，也符合我们的认知。
考察一下\(M_G^T\)对\(F\)的作用。
\(H'(x)=\sum_{i=0}^n\sum_{j=0}^if_ig_jx^{i-j}\)。我们定义这种运算叫做多项式卷积的转置（减法卷积），记作\(H'(x)=F(x)\times^TG(x)\)。显然卷积的转置不满足交换律。
注意这里为了使得转置前后矩阵乘法有意义，\(M_G\)舍去了\(n+1\)行之后的内容。我们规定\(n\)次多项式和\(m\)次多项式的卷积是\(n+m\)次的，而他们卷积的转置是\(n-m\)次的。
多项式综合操作的转置
现在我们对多项式\(F(x)\)进行了一系列加减乘除操作，那么如果把这些操作综合起来转置再对\(F(x)\)作用，结果是什么呢？
-
首先加减法可以看成是列向量的加减，因为乘除操作对加法具有分配律，所以完全可以将所有加减号全部拆开。也就是说加减法对\(F(x)\)根本没有影响，不在考虑范围内。多项式的加减法本身时间复杂度并不高可以直接做，所以这个结论也符合我们的认知。
-
其次是取模操作，\(\bmod x^{n+1}\)相当于把\(n+1\)次及其以后的项系数乘\(0\)。这相当于第二类初等矩阵的作用，因此转置前后没变。
-
最后是卷积（除法是与多项式的逆卷积），卷积的转置是新操作，上面分析过了。
这些操作综合起来转置的话，先要确定被操作的多项式是哪个，再将操作次序颠倒。
例如对\(H(x)\times(G(x)\times F(x)+A(x))\)中\(F(x)\)接受的操作转置，那么就是\(F(x)\times^TG(x)\times^TH(x)\)。
另外分析操作的转置的时候既可以将\(F(x)\)看成多项式又可以将\(F\)看成列向量。总之哪个方便用哪个分析。
多项式多点求值问题
有一个\(n\)次多项式\(F(x)=\sum_{i=0}^nf_ix^i\)，给出\(n+1\)个不同的数\(a_0,a_1,\dots,a_n\)，求\(b_0=F(a_0),b_1=F(a_1),\dots,b_n=F(a_n)\)。
不难发现规定求\(n+1\)个点值得出的结论适用于求任意个点值。
考虑矩阵
则我们要求的就是\(B=AF\)。虽然\(AF\)不好求，但是我们发现\(B'=A^TF\)好求。
令形式幂级数\(B'(x)=\sum_{i=0}^nb'_ix^i\)，则
\(B'(x)\)可通过分治FFT求出。令\(P_{i,i}(x)=f_i, Q_{i,i}(x)=1-a_ix\)，则
最后\(B'(x)=\frac{P_{0,n}(x)}{Q_{0,n}(x)}\)。
既然求\(A^TF\)时通过一系列多项式操作得到了\(B'(x)\)，那么我们求\(AF\)的时候只需要把这些操作颠倒顺序再转置就可得到\(B(x)\)。现在的关键在于\(A^TF\)究竟对\(F\)进行了那些操作。
-
首先\(P_{0,n}(x)\)与\(F(x)\)有关，而\(Q_{0,n}^{-1}(x)\)与\(F(x)\)无关。最后一个操作是\(P_{0,n}(x)\times Q_{0,n}^{-1}(x)\)，所以转置后的第一个操作是\(B_{0,n}(x)=F(x)\times^TQ_{0,n}^{-1}(x)\)。
-
然后在递归的过程中，被操作的\(F\)的某个元\(f_i\)只与\(P_{l,m}(x)\)和\(P_{m+1,r}(x)\)的一个有关。因此分左右两侧递归，左侧传参\(B_{l,m}(x)=B_{l,r}(x)\times^TQ_{m+1,r}(x)\)，右侧传参\(B_{m+1,r}(x)=B_{l,r}(x)\times^TQ_{l,m}(x)\)。
-
最后递归到叶子的时候，\(F\)的转置操作已经全部完成。也就是说这时候叶结点的参值就是我们要求的点值。
Luogu5050 多项式多点求值
我们发现\(B_{0,n}(x)=F(x)\times^TQ_{0,n}^{-1}(x)\)这个式子中\(B_{0,n}(x)\)的次数并不是\(n\)。
仔细分析发现\(B'(x)=P_{0,n}(x)\times Q_{0,n}^{-1}(x)\bmod x^{n+1}\)，也就是说最后一步是取模。（注意分治FFT合并的过程中并没有取模）
这个取模的操作等于是将\(x^{n+1}\)项至\(x^{2n}\)项的系数乘以\(0\)。转置后操作不变，我们就应该先将\(F(x)\)用\(0\)补齐到\(2n\)次项。
总时间复杂度\(O(n\log^2 n)\)，因为不涉及多项式取模常数小很多。
而这个做法的优势在于做减法卷积时可以利用循环卷积NTT长度减少一半。
穷尽手段减少NTT次数从而进一步卡常：https://www.cnblogs.com/chasedeath/p/13073178.html
据说这个做法能1s做1e6，但是如果用
vector封装那么连多项式乘法都不能做1e6。😓
// L=ceil(log2(2*n))
constexpr int L=17, N=1<<L;
using Poly=vector<int>;
int omg[2][N+1], rev[N+1];
int fac[N+1], inv[N+1], ifac[N+1];
void initNTT(){
omg[0][0]=1, omg[0][1]=fastPow(3, (MOD-1)/N);
omg[1][0]=1, omg[1][1]=fastPow(omg[0][1], MOD-2);
rev[0]=0, rev[1]=1<<(L-1);
fac[0]=fac[1]=1;
inv[0]=inv[1]=1;
ifac[0]=ifac[1]=1;
for(int i=2; i<=N; ++i){
omg[0][i]=mul(omg[0][i-1], omg[0][1]);
omg[1][i]=mul(omg[1][i-1], omg[1][1]);
rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
fac[i]=mul(fac[i-1], i);
inv[i]=mul(MOD-MOD/i, inv[MOD%i]);
ifac[i]=mul(ifac[i-1], inv[i]);
}
}
void NTT(Poly &a, int dir){
int lim=a.size(), len=log2(lim);
for(int i=0; i<lim; ++i){
int r=rev[i]>>(L-len);
if(i<r) swap(a[i], a[r]);
}
for(int i=1; i<lim; i<<=1)
for(int j=0; j<lim; j+=i<<1)for(int k=0; k<i; ++k){
int t=mul(omg[dir][N/(i<<1)*k], a[j+i+k]);
a[j+i+k]=add(a[j+k], MOD-t), a[j+k]=add(a[j+k], t);
}
if(dir)for(int i=0; i<lim; ++i)
a[i]=mul(a[i], inv[lim]);
}
Poly operator~(Poly a){
int n=a.size();
Poly b={fastPow(a[0], MOD-2)};
a.resize(1<<(int)ceil(log2(n)));
for(int lim=2; lim<2*n; lim<<=1){
Poly c(a.begin(), a.begin()+lim);
c.resize(lim<<1), NTT(c, 0);
b.resize(lim<<1), NTT(b, 0);
for(int i=0; i<lim<<1; ++i) b[i]=mul(2+MOD-mul(c[i], b[i]), b[i]);
NTT(b, 1), b.resize(lim);
}
return b.resize(n), b;
}
Poly log(Poly a){
int n=a.size();
Poly b=~a;
for(int i=0; i<n-1; ++i) a[i]=mul(a[i+1], i+1);
a.resize(n-1);
int lim=1<<(int)ceil(log2(2*n-2));
a.resize(lim), NTT(a, 0);
b.resize(lim), NTT(b, 0);
for(int i=0; i<lim; ++i) a[i]=mul(a[i], b[i]);
NTT(a, 1), a.resize(n);
for(int i=n-1; i>=1; --i) a[i]=mul(a[i-1], inv[i]);
return a[0]=0, a;
}
// a[0]=0
Poly exp(Poly a){
int n=a.size();
Poly b={1};
a.resize(1<<(int)ceil(log2(n)));
for(int lim=2; lim<2*n; lim<<=1){
b.resize(lim); Poly c=log(b);
c[0]=add(1+a[0], MOD-c[0]);
for(int i=1; i<lim; ++i) c[i]=add(a[i], MOD-c[i]);
c.resize(lim<<1), NTT(c, 0);
b.resize(lim<<1), NTT(b, 0);
for(int i=0; i<lim<<1; ++i) b[i]=mul(c[i], b[i]);
NTT(b, 1), b.resize(lim);
}
return b.resize(n), b;
}
Poly operator*(Poly a, Poly b){
int n=a.size()+b.size()-1, lim=1<<(int)ceil(log2(n));
a.resize(lim), NTT(a, 0);
b.resize(lim), NTT(b, 0);
for(int i=0; i<lim; ++i) a[i]=mul(a[i], b[i]);
NTT(a, 1), a.resize(n);
return a;
}
Poly deMul(Poly a, Poly b){
int n=a.size(), m=b.size();
if(n<m) return {};
reverse(b.begin(), b.end());
int lim=1<<(int)ceil(log2(n)); // n+m-1 -> n
a.resize(lim), b.resize(lim);
NTT(a, 0), NTT(b, 0);
for(int i=0; i<lim; ++i) a[i]=mul(a[i], b[i]);
NTT(a, 1);
for(int i=0; i<=n-m; ++i) a[i]=a[i+m-1];
a.resize(n-m+1);
return a;
}
Poly q[N];
int ans[N];
#define lc (x<<1)
#define rc (x<<1|1)
#define mid ((l+r)>>1)
void goUp(int x, int l, int r, const Poly&a){
if(l==r){
q[x]={1, (MOD-a[l])%MOD};
return;
}
goUp(lc, l, mid, a), goUp(rc, mid+1, r, a);
q[x]=q[lc]*q[rc];
}
void goDown(int x, int l, int r, const Poly&b){
if(l==r) {ans[l]=b[0]; return;}
goDown(lc, l, mid, deMul(b, q[rc])),
goDown(rc, mid+1, r, deMul(b, q[lc]));
}
void multiPoint(Poly f, Poly a){
int n=max(f.size(), a.size());
f.resize(n<<1), a.resize(n);
goUp(1, 0, n-1, a);
goDown(1, 0, n-1, deMul(f, ~q[1])); // length=3*n -> 2*n
}
int main(){
initNTT();
int n=read<int>(), m=read<int>();
Poly f(n+1), a(m);
for(int i=0; i<=n; ++i) read(f[i]);
for(int i=0; i<m; ++i) read(a[i]);
multiPoint(f, a);
for(int i=0; i<m; ++i) printf("%d\n", ans[i]);
return 0;
}
电子科大校赛2022N 简单算术
给出\(n\)个正整数\(a_1,a_2,\dots,a_n\)，求
\(n\leq 10^5\)。
题解
放眼全局，形如\((a_i+a_i)\)的项只出现一次，而\((a_i+a_j)\)（\(i\neq j\)）出现了两次。因为原式形式不统一，不便于进行推导，所以对原式进行拆分
设\(F(l,r)=\prod_{i=l}^r\prod_{j=l+1}^r(a_i+a_j)\)，则问题转化成了计算\(F(1,n)\)。
这种任意两项都组合一个贡献的形式和排列逆序对非常相似
考虑分治。令\(m=\lfloor\frac{l+r}{2}\rfloor\)
考虑每个\(a_i\)的所有贡献，发现它们十分类似。
若设\(G_{m+1, r}(x)=\prod_{j=m+1}^r(x+a_j)\)，则要求的就是\(\prod_{i=l}^mG_{m+1, r}(a_i)\)。
这类似于多项式变量添系数求积的问题。可惜我们要求的不是多项式而是具体值，不能用这个技巧。把这个连乘换成求和我将绝杀，可惜换不得。
直接调用多项式多点求值，时间复杂度\(T(n)=2T(n/2)+O(n\log^2 n)=O(n\log^3 n)\)。鉴于转置原理多点求值的常数非常小，所以是可以过的。
总而言之，这题是一个板子题，没有什么有趣的地方。而ACM可以带板子，完全没有代码难度，就看选手平时积累了。
电子科大300个队中在考场上AC的有2个队，颇具实力。
Poly q[N];
int val[N];
#define lc (x<<1)
#define rc (x<<1|1)
#define mid ((l+r)>>1)
void goUp(int x, int l, int r, const Poly&a){
if(l==r){
q[x]={1, (MOD-a[l])%MOD};
return;
}
goUp(lc, l, mid, a), goUp(rc, mid+1, r, a);
q[x]=q[lc]*q[rc];
}
void goDown(int x, int l, int r, const Poly&b){
if(l==r) {val[l]=b[0]; return;}
goDown(lc, l, mid, deMul(b, q[rc])),
goDown(rc, mid+1, r, deMul(b, q[lc]));
}
void multiPoint(Poly f, Poly a){
int n=max(f.size(), a.size());
f.resize(n<<1), a.resize(n);
goUp(1, 0, n-1, a);
goDown(1, 0, n-1, deMul(f, ~q[1])); // length=3*n -> 2*n
}
Poly a, f[N];
int solve(int x, int l, int r){
if(l==r){
f[x]={a[l], 1};
return 1;
}
int ans=mul(solve(lc, l, mid), solve(rc, mid+1, r));
f[x]=f[lc]*f[rc];
multiPoint(f[rc], Poly(a.begin()+l, a.begin()+mid+1));
for(int i=0; i<=mid-l; ++i) ans=mul(ans, val[i]);
return ans;
}
int main(){
freopen("N.in", "r", stdin);
freopen("N.out", "w", stdout);
initNTT();
int n=read<int>();
a.resize(n);
for(int i=0; i<n; ++i) read(a[i]);
int ans=1;
for(int i=0; i<n; ++i) ans=mul(ans, 2*a[i]);
ans=mul(ans, fastPow(solve(1, 0, n-1),2));
printf("%d\n", ans);
cerr<<(double)clock()/CLOCKS_PER_SEC<<endl;
return 0;
}
另外我发现\(n=10^5\)且数据随机时，答案大概率是\(0\)。这或许可以衍生出一些乱搞做法……
