转置原理

转置原理

初等矩阵

众所周知，矩阵的初等变换有以下三种：

  1. 交换第 ii 行（列）和第 jj 行（列）。
  2. 将第 ii 行（列）乘上 kk 。
  3. 将第 ii 行（列）的 kk 倍对应加到第 jj 行（列）上。

同时，我们还知道单位矩阵 II ：

[10⋯001⋯0⋮⋮⋮⋮00⋯1]\begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}

对于任意矩阵 AA ，有 A×I=I×A=AA\times I=I\times A=A 。

若将初等变换作用在 II 上，就可以得到初等矩阵 EE ：

  • 对于初等行变换，则 E×AE \times A 为对 AA 进行初等行变换。
  • 对于初等列变换，则 A×EA \times E 为对 AA 进行初等列变换。

算法的转置

如果算法可以当作方阵 AA ,

如果算法的输入可以当作向量 vv ，输出为向量 AvAv ，

则称这个算法为线性算法。

而对于算法矩阵的转置矩阵 ATA^T ，

如果算法的输入可以当作向量 vv ，输出就为向量 ATvA^Tv 。

转置算法有一些重要的性质：

引理 1 （转置的性质）：

(∏i=1kAi)T=∏i=k1AiT(\prod_{i=1}^k A_i)^T=\prod_{i=k}^1 A_i^T

引理 2 （初等矩阵转置的性质）：

对于一个初等矩阵 EE ：

如果 EE 是 II 进行前两种初等变换，则 ET=EE^T=E ；

如果 EE 是 II 进行第三种初等变换，若第三种初等变换执行 Ii=Ii+kIjI_i=I_i+kI_j ，则有 ETE_T 为 II 执行 Ij=Ij+kIiI_j=I_j+kI_i 。

线性算法判定

以下为随便写的。

将输入，输出，算法流程中所使用的辅助变量拼成一个向量作为输入。

为保证行文统一，以下均称向量中的量为变量，矩阵中的量为常量。

则线性算法流程中应只包含如下三种操作：

  1. 交换两个变量。
  2. 将某个变量 xx 变为 mxmx ， mm 为常量。
  3. 对于变量 x,yx,y 和常量 mm ，将 xx 变为 x+myx+my 。

用这个判定，我们可以比较方便的做一些事。

线性算法与转置算法

对于一个线性算法 AA ，由判定可以保证其可分解为 A=∏i=1kEkA=\displaystyle\prod_{i=1}^k E_k ， EkE_k 均为初等矩阵。

然后其转置算法 ATA^T ，有 AT=∏i=k1EkTA_T=\displaystyle\prod_{i=k}^1 E_k^T 。

由转置算法性质，我们得到了转置算法的算法流程：

  1. 将线性算法的操作顺序倒序。
  2. 将其中形如 x=x+myx=x+my 的操作变为 y=y+mxy=y+mx 。

这是后面一切的基础。

多项式乘法（卷积）的转置

问题是对于长度 n+1n+1 的向量 aa 和长度 m+1m+1 的向量 bb ，求长度 n+m+1n+m+1 的向量 cc 满足：

ci=∑j=0iajbi−jc_i=\sum_{j=0}^i a_jb_{i-j}

可以很快的写出算法结果：

1
2
3
4
5
for(int i=0;i<n+m+1;i++){
    for(int j=min(0,i-m);j<=max(i,n);j++){
        c[i]=c[i]+a[j]*b[i-j];
    }
}

如果将向量 bb 视作常量，那么可以轻松写出转置算法：

1
2
3
4
5
for(int i=n+m;i>=0;i--){
    for(int j=max(i,n);j>=min(0,i-m);j--){
        a[j]=a[j]+c[i]*b[i-j];
    }
}

即： ai=∑kckbk−ia_i=\displaystyle\sum_k c_kb_{k-i}

这是一个减法卷积的形式。具体来说，将 bb 反转后正常卷积 cc 和 bb 之后，设结果为 dd ，则有 ai=di+m−1a_i=d_{i+m-1} 。

伪代码：

1
2
3
4
reverse(b);
d=poly_multi(b,c);
for i in {0,n-1}:
    a[i]=d[i+m-1];

转置原理

对于一个线性算法 AA ，如果对于输入 vv ，可以在 O(T(n))O(T(n)) 内计算 AvAv ，则必然可以在 O(T(n))O(T(n)) 内计算 ATvA^Tv 。

反之，如果对于输入 vv ，可以在 O(T(n))O(T(n)) 内计算 ATvA^Tv ，则必然可以在 O(T(n))O(T(n)) 内计算 AvAv 。

这就是转置原理，又称特勒根原理。

多项式多点求值

给定多项式 F(x)=∑i=0n−1fixiF(x)=\displaystyle\sum_{i=0}^{n-1} f_ix^i ，和 nn 个数 a0,a1,⋯ ,an−1a_0,a_1,\cdots,a_{n-1} ，求 F(a1),F(a2),F(a3),⋯ ,F(an−1)F(a_1),F(a_2),F(a_3),\cdots,F(a_{n-1}) 。

其实容易发现求的东西是一个点值范德蒙德矩阵 AA 乘系数向量 ff 。

[1a0⋯a0n−11a1⋯a1n−1⋮⋮⋮⋮1an−1⋯an−1n−1]×[f0f1⋮fn−1]\begin{bmatrix} 1 & a_0 & \cdots & a_0^{n-1} \\ 1 & a_1 & \cdots & a_1^{n-1} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & a_{n-1} & \cdots & a_{n-1}^{n-1} \end{bmatrix} \times \begin{bmatrix} f_0\\ f_1\\ \vdots\\ f_{n-1} \end{bmatrix}

将矩阵转置，令结果的第 kk 项为 yky_k ，容易求出 yy 的生成函数： ∑i=0n−1fi1−aix\displaystyle\sum_{i=0}^{n-1}\frac{f_i}{1-a_ix}

这个东西可以使用分治 FFT 求解。

具体来说，维护分子 FF 和分母 GG ，

则有：

F=FlGr+FrGlG=GlGr\begin{aligned} F& = F_lG_r+F_rG_l\\ G &= G_lG_r \end{aligned}

然后我们将这个的线性算法形式表示出来（这里的其实是多项式多点求值的转置问题的线性算法，当然由于转置两次等于没转置，所以之后我们表示的这个算法的转置算法就是原问题）：

  • 初始化 GG （由于 aa 是常量，所以这个不受到转置影响）。
  • 初始化 Fl,lF_{l,l} 。
  • 自下而上：
    • 区间中点设为 mm 。
    • Fl,m=Fl,m×Gm+1,rF_{l,m}=F_{l,m}\times G_{m+1,r}
    • Fm+1,r=Fm+1,r×Gl,mF_{m+1,r}=F{m+1,r}\times G_{l,m}
    • Fl,r=Fl,r+Fl,mF_{l,r}=F_{l,r}+F{l,m}
    • Fl,r=Fl,r+Fm+1,rF_{l,r}=F_{l,r}+F_{m+1,r}
  • F0,n−1=F0,n−1×G0,n−1−1F_{0,n-1}=F_{0,n-1} \times G_{0,n-1}^{-1}

转置这个算法：

  • 初始化 GG 。
  • 初始化 F0,n−1=f×TG0,n−1−1F_{0,n-1}=f\times^T G_{0,n-1}^{-1} ，其中 ×T\times^T 为多项式乘法的转置。
  • 自上而下：
    • 区间中点设为 mm 。
    • Fl,m=Fl,m+Fl,rF_{l,m}=F_{l,m}+F_{l,r}
    • Fm+1,r=Fm+1,r+Fl,rF_{m+1,r}=F_{m+1,r}+F_{l,r}
    • Fl,m=Fl,m×TGm+1,rF_{l,m}=F_{l,m}\times^T G_{m+1,r}
    • Fm+1,r=Fl,m×TGl,mF_{m+1,r}=F_{l,m}\times^T G_{l,m}
  • 赋值答案 ansl=Fl,lans_l=F_{l,l} 。

这样就可在 O(nlog⁡2n)O(n\log^2n) 内计算多项式多点求值了。

题目： 多项式多点求值 。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include <bits/stdc++.h>
#define poly vector<int>
using namespace std;
const int N=65010<<4,P=998244353;
int n,m,r[N];
poly a,b,ans;
int kpow(int a,int b){
    int t=1;
    while(b){
        if(b&1)t=1ll*t*a%P;
        a=1ll*a*a%P;
        b>>=1;
    }
    return t;
}
void NTT(poly &f){
    int n=f.size();
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)*n>>1);
    for(int i=0;i<n;i++)if(i<r[i])swap(f[i],f[r[i]]);
    for(int i=2;i<=n;i<<=1){
        int G=kpow(3,(P-1)/i);
        for(int j=0;j<n;j+=i){
            int g=1;
            for(int k=j;k<j+(i>>1);k++){
                int x=f[k],y=1ll*f[k+(i>>1)]*g%P;
                f[k]=(x+y)%P;
                f[k+(i>>1)]=(x-y+P)%P;
                g=1ll*g*G%P;
            }
        }
    }
}
void INTT(poly &f){
    NTT(f),reverse(f.begin()+1,f.end());
    int in=kpow(f.size(),P-2);
    for(int i=0;i<f.size();i++)f[i]=1ll*f[i]*in%P;
}
poly MTT(poly a,poly b){
    int n=a.size()+b.size()-1,m=1;
    while(m<n)m<<=1;
    a.resize(m),b.resize(m);NTT(a);NTT(b);
    for(int i=0;i<a.size();i++)a[i]=1ll*a[i]*b[i]%P;
    INTT(a);
    a.resize(n);
    return a;
}
poly IMTT(poly a,poly b){
    int n=a.size(),m=b.size();
    reverse(b.begin(),b.end());b=MTT(a,b);
    for(int i=0;i<n;i++)a[i]=b[i+m-1];
    return a;
}
poly ta,tb;
poly inv(poly &a){
    int n=a.size();
    poly b(1);b[0]=kpow(a[0],P-2);
    for(int tt=2;(tt>>1)<n;tt<<=1){
        ta.resize(tt);b.resize(tt);
        for(int i=0;i<tt;i++)ta[i]=a[i];
        tb=MTT(b,b);ta=MTT(ta,tb);
        for(int i=0;i<tt;i++)b[i]=(2ll*b[i]%P-ta[i]+P)%P;
    }
    b.resize(n);
    return b;
}
poly F[N<<2];
void build(int p,int l,int r){
    if(l==r){
        F[p].resize(2);
        F[p][0]=1,F[p][1]=P-b[l];
        return ;
    }
    int mid=(l+r)>>1;
    build(p<<1,l,mid);
    build(p<<1|1,mid+1,r);
    F[p]=MTT(F[p<<1],F[p<<1|1]);
}
void PME(int p,int l,int r,poly f,poly &g){
    f.resize(r-l+1);
    if(l==r)return (void)(g[l]=f[0]);
    int mid=(l+r)>>1;
    PME(p<<1,l,mid,IMTT(f,F[p<<1|1]),g);
    PME(p<<1|1,mid+1,r,IMTT(f,F[p<<1]),g);
}
int main(){
    cin>>n>>m;
    n++;
    int M=max(n,m);
    a.resize(M+1),b.resize(M),ans.resize(M);
    for(int i=0;i<n;i++)cin>>a[i];
    for(int i=0;i<m;i++)cin>>b[i];
    build(1,0,M-1);
    PME(1,0,M-1,IMTT(a,inv(F[1])),ans);
    for(int i=0;i<m;i++)cout<<ans[i]<<"\n";
    return 0;
}