转置原理多点求值 和 电子科大校赛2022N 简单算术

转置原理多点求值

rushcheyo 《转置原理及其应用》 https://jkloverdcoi.github.io/2020/08/04/转置原理及其应用/slide.pdf
https://blog.csdn.net/weixin_43960287/article/details/122691212

OIer学习的功利性很强情有可原，但是我要的是严谨的证明。

转置原理

若矩阵\(A\)可逆，则\(A\)可表示为初等矩阵的积。

令\(A=P_1P_2\cdots P_n\)，则\(A^T=P_n^T\cdots P_2^TP_1^T\)。而三种初等矩阵的转置分别为

1. \(E(i(k))^T=E(i(k))\)不变。

2. \(E(i,j)^T=E(i,j)\)不变。

3. \(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\)。即

\[F=\begin{bmatrix} f_0\\ f_1\\ \vdots\\ f_n \end{bmatrix}, M_G=\begin{bmatrix} g_0 & 0 & 0 & \cdots & 0\\ g_1 & g_0 & 0 & \cdots & 0\\ g_2 & g_1 & g_0 & \cdots & 0\\ \vdots & \vdots & \vdots & & \vdots\\ g_n & g_{n-1} & g_{n-2} & & g_0 \end{bmatrix}, H=M_GF \]

注意\(M_G\)只写了\(n+1\)行。

模拟下三角矩阵通过行变换转化成单位矩阵的过程手动拆分一下\(M_G\)，

\[\begin{split} \because\left(\prod_{i=1}^{n+1}E\left(i(\frac{1}{g_0})\right)\right)\left(\prod_{i=1}^{n+1}\prod_{j=i+1}^{n+1}E\left(j,i(-\frac{g_{j-i}}{g_0})\right)\right)M_G=E_{n+1}\\ \therefore M_G=\left(\prod_{i=1}^{n+1}\prod_{j=i+1}^{n+1}E\left(j,i(\frac{g_{j-i}}{g_0})\right)\right)\left(\prod_{i=1}^{n+1}E(i(g_0))\right)E_{n+1} \end{split} \]

考察一下\(M_G^T\)的形式

\[\begin{split} M_G^T=\left(\prod_{i=1}^{n+1}E(i(g_0))\right)\left(\prod_{i=1}^{n+1}\prod_{j=i+1}^{n+1}E\left(i,j(\frac{g_{j-i}}{g_0})\right)\right)E_{n+1}\\ =\begin{bmatrix} g_0 & g_1 & g_2 & \cdots & g_n\\ 0 & g_0 & g_1 & \cdots & g_{n-1}\\ 0 & 0 & g_0 & \cdots & g_{n-2}\\ \vdots & \vdots & \vdots & & \vdots\\ 0 & 0 & 0 & \cdots & g_0 \end{bmatrix} \end{split} \]

所以即使\(M_G\)不可逆（\(g_0=0\)），我们的\(M_G\)仍然是满足使用转置定理的条件的。而\(M_G^T\)本质上就是将\(M_G\)行列交换了，也符合我们的认知。

考察一下\(M_G^T\)对\(F\)的作用。

\[M_G^TF=\begin{bmatrix} \sum_{i=0}^nf_ig_i\\ \sum_{i=1}^nf_ig_{i-1}\\ \sum_{i=2}^nf_ig_{i-2}\\ \vdots\\ f_ng_0 \end{bmatrix}=H' \]

\(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\)个点值得出的结论适用于求任意个点值。

考虑矩阵

\[A=\begin{bmatrix} 1 & a_0 & \cdots & a_0^n\\ 1 & a_1 & \cdots & a_1^n\\ \vdots & \vdots & & \vdots\\ 1 & a_n & \cdots & a_n^n \end{bmatrix}, F=\begin{bmatrix} f_0\\ f_1\\ \vdots\\ f_n \end{bmatrix} \]

则我们要求的就是\(B=AF\)。虽然\(AF\)不好求，但是我们发现\(B'=A^TF\)好求。

\[\begin{split} A^TF=\begin{bmatrix} 1 & 1 & \cdots & 1\\ a_0 & a_1 & \cdots & a_n\\ \vdots & \vdots & & \vdots\\ a_0^n & a_1^n & \cdots & a_n^n \end{bmatrix} \begin{bmatrix} f_0\\ f_1\\ \vdots\\ f_n \end{bmatrix}\\ =\sum_{i=0}^nf_i\begin{bmatrix} 1\\ a_i\\ \vdots\\ a_i^n \end{bmatrix}=\begin{bmatrix} b'_0\\ b'_1\\ \vdots\\ b'_n \end{bmatrix}=B' \end{split} \]

令形式幂级数\(B'(x)=\sum_{i=0}^nb'_ix^i\)，则

\[\begin{split} B'(x)=\sum_{i=0}^n\frac{f_i}{1-a_ix}\mod x^{n+1}\\ =\frac{\sum_{i=0}^nf_i\prod_{j\neq i}(1-a_jx)}{\prod_{i=0}^n(1-a_ix)}\mod x^{n+1} \end{split} \]

\(B'(x)\)可通过分治FFT求出。令\(P_{i,i}(x)=f_i, Q_{i,i}(x)=1-a_ix\)，则

\[\begin{split} Q_{l,r}(x)=Q_{l,m}(x)\times Q_{m+1,r}(x)\\ P_{l,r}(x)=P_{l,m}(x)\times Q_{m+1,r}(x)+Q_{l,m}(x)\times P_{m+1,r}(x) \end{split} \]

最后\(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\)，求

\[\prod_{1\leq i,j\leq n}(a_i+a_j)\mod 998244353 \]

\(n\leq 10^5\)。

题解

放眼全局，形如\((a_i+a_i)\)的项只出现一次，而\((a_i+a_j)\)（\(i\neq j\)）出现了两次。因为原式形式不统一，不便于进行推导，所以对原式进行拆分

\[\prod_{1\leq i,j\leq n}(a_i+a_j)=\left(\prod_{i=1}^n2a_i\right)\left(\prod_{i=1}^n\prod_{j=i+1}^n(a_i+a_j)\right)^2 \]

设\(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\)

\[F(l,r)=F(l,m)\times F(m+1, r)\times\prod_{i=l}^m\prod_{j=m+1}^r(a_i+a_j) \]

考虑每个\(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\)。这或许可以衍生出一些乱搞做法……