本文简单介绍一下基本多项式理论。
持续更新中。
因为是简单介绍，只说明做法，不给予严格证明，以及前人是如何想出的这个方法。
多项式乘法
多项式乘法与和卷积
\(A=\sum\limits_{i=0}^{n}a_ix^i,B=\sum\limits_{i=0}^{n}b_ix^i\)。
显然有 \(A\times
B=\sum\limits_{i=0}^n\sum\limits_{j=0}^na_ib_jx^{i+j}=\sum\limits_{d=0}^n\sum\limits_{i+j=d}^na_ib_jx^d\)。
\(C=A\times
B,C(d)=\sum\limits_{i+j=d}^na_ib_jx^d\)。
这玩意就是和卷积，可见和卷积与多项式乘法有着千丝万缕的关系。
快速傅里叶变换
众所周知，一个 \(n-1\)
次多项式（函数）可以用 \(n\)
个不同的点唯一确定出来。
而如果考虑两个多项式 \(f,g\)，当\(x=x_0\) 时，\(f\times g\ (x_0)=f(x_0)\times g(x_0)\)
。
这个很好理解。
所以点值表达乘积就是直接将相同横坐标的 \(y\) 乘。
那么如何将多项式 \(f\)
转成点值。
代入复数，用 \(n\) 次单位根做 \(x\)。
即 \(x^n=1\) 的 \(n\) 个复数解，记录以 \(x\) 轴正方向，逆时针第一个复根为 \(w\)，所有单位根就是 \(w^0,w^1,\cdots,w^{n-1}\)。
目标：求 \(\forall \
i,F(k)=\sum\limits_{j=0}^{n-1}a_j(w^i)^j\)。
我们将次方奇数偶数分开考虑 \(F(k)=\sum\limits_{j=0}^{n}a_j(w^k)^j=\sum\limits_{j=0}^{n/2-1}a_{2j}(w^k)^{2j}+w^k\sum\limits_{j=0}^{n/2-1}a_{2j+1}(w^k)^{2j}\)。
\(F(k)=\sum\limits_{j=0}^{n}a_j(w^k)^j=\sum\limits_{j=0}^{n/2-1}a_{2j}(w^{2k})^{j}+w^k\sum\limits_{j=0}^{n/2-1}a_{2j+1}(w^{2k})^{j}\)。
由于 \((w^k)^2=(w^{k+n/2})^2\)
（\(w^n=1\))。
所以 \(F(k)\) 与 \(F(k+n/2)\) 只差一点点，就是 \(w^k\) 项的正负号 （\(w^k=-w^{k+n/2}\)）。
原本你要求 \(F(w^0),F(w^1),\cdots
F(w^{n-1})\),现在只要求 \(F(w^0),F(w^1),\cdots F(w^{n/2-1})\)。
范围直接缩小一倍。那么对于接下来的 \(F_0(w^{2k}),F_1(w^{2k})\)（系数是奇数(1)还是偶数(0)）也可以递归下求。
递归效率缓慢，能够优化。
发现本质上每次将 \(k\),缩小一半的操作其实是想要最后把 \(w^0\) 算出，然后计算出别的。
可以如图按照一定顺序把 \(a_i\times
(w^0)^i\) 按顺序排开，然后合并。
大概。。。
这个顺序通过观察 发现是二进制反过来
，所以能快速求出。
代码解读一下。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| const pie=acos(-1);
inline void fft(complex *a)
{
for(int i=0;i<lim;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int mid=1;mid<lim;mid=mid<<1)
{
complex wn(cos(p/mid),sin(p/mid));
for(int r=mid<<1,j=0;j<lim;j+=r)
{
complex w(1.0,0.0);
for(int k=0;k<mid;k++,w=w*wn)
{
complex x=a[k+j],y=w*a[k+j+mid];
a[k+j]=x+y;
a[k+j+mid]=x-y;
}
}
}
}
|
一般人不用点值表达，所以必须要有逆傅里叶变换。
不难发现 由于范德蒙德矩阵在 \(A_{i,j}=w^{ij}\) 有着巨NB的性质，\(A^{-1}_{i,j}=w^{-ij}\)。 （前面那个 \(-1\) 的意思是逆矩阵）
所以只要把上文代码中的 wn(cos(p/mid),sin(p/mid)) 改成
wn(cos(-p/mid),sin(-p/mid)) 即可。
证明是个构造证明，我也不会从暴力手算 \(A_{i,j}\) 逆矩阵出发证明，而只会验证 \(A\times B=\epsilon\),\(B\) 即为你构造的那个逆矩阵。
快速数论变换
把单位根换成原根有相同性质。
再此贴出现在的Poly模板；
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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
| #include<iostream>
#include<cmath>
#include<vector>
#include<algorithm>
#include<assert.h>
using namespace std;
#define int long long
const int mod=998244353;
inline int qpow(int a,int b){
int ret=1;
while(b){
if(b&1) ret=ret*a%mod;
a=a*a%mod;b>>=1;
}return ret;
}
inline void up(int &x){
if(x>=mod) x-=mod;
}
int *rev[23],*wn[23];
inline void init(){
for(int mid=1,lg=0;mid<(1<<20);mid<<=1,lg++){
wn[lg]=new int[mid+1];
wn[lg][0]=1;
int w=qpow(3,(mod-1)/(mid<<1));
for(int k=1;k<mid;k++){
wn[lg][k]=wn[lg][k-1]*w%mod;
}
}
}
struct poly{
#define rep(i,l,r) for(int i=l;i<=r;i++)
vector<int> a;
inline poly(){}
inline poly(vector<int> v){a=v;}
inline poly(int n){a.resize(n+1);}
inline int size() const{
return a.size()-1;
}
inline void resize(int lim){
a.resize(lim+1);
}
inline void ntt(int lim){
resize(lim-1);
int len=__lg(lim);
if(rev[len]==nullptr){
rev[len]=new int [lim+1];
rev[len][0]=0;
for(int i=1;i<lim;i++) {
rev[len][i]=(rev[len][i>>1]>>1)|((i&1)?(lim>>1):0);
}
}
for(int i=0;i<lim;i++) if(rev[len][i]<i) swap(a[rev[len][i]],a[i]);
for(int mid=1,lg=0;mid<lim;mid<<=1,lg++){
for(int j=0;j<lim;j+=(mid<<1)){
for(int k=0;k<mid;k++){
int x=a[j+k],y=a[j+k+mid]*wn[lg][k]%mod;
a[j+k]=(x+y);if(a[j+k]>=mod) a[j+k]-=mod;
a[j+k+mid]=(x-y);if(a[j+k+mid]<0) a[j+k+mid]+=mod;
}
}
}
}
inline void intt(int lim){
resize(lim-1);
reverse(a.begin()+1,a.end());
ntt(lim);
int inv=qpow(lim,mod-2);
for(int i=0;i<lim;i++) a[i]=a[i]*inv%mod;
}
inline int& operator [](int x){
return a[x];
}
inline poly operator +(const poly&b){
int sza=size(),szb=b.size();
int lim=max(sza,szb);
poly c(lim);
rep(i,0,lim) {
int x=0,y=0;
if(i<=sza) x=a[i];
if(i<=szb) y=b.a[i];
c[i]=(x+y),up(c[i]);
}
return c;
}
inline poly operator -(const poly&b){
int sza=size(),szb=b.size();
int lim=max(sza,szb);
poly c(lim);
rep(i,0,lim) {
int x=0,y=0;
if(i<=sza) x=a[i];
if(i<=szb) y=b.a[i];
c[i]=(x-y+mod),up(c[i]);
}
return c;
}
inline poly operator *(const poly&b){
poly c=b,d=(*this);
int len=c.size()+(d.size());
int lim=1;
while(lim<=len) lim<<=1;
d.ntt(lim),c.ntt(lim);
rep(i,0,lim-1) c[i]=c[i]*d[i]%mod;
c.intt(lim);
c.resize(len);
return c;
}
inline poly operator *(int b){
poly c=(*this);
for(auto &x:c.a) x=(x*b)%mod;
return c;
}
inline poly inv(int lim){
if(lim==1){
poly x(0);x[0]=qpow(a[0],mod-2);
return x;
}
poly c=(*this),ans;
c.resize(lim-1);
int l=(lim+1)/2;
poly inv=c.inv(l);
ans=(inv*2)-(inv*inv)*c;
ans.resize(lim-1);
return ans;
}
inline poly Sqrt(int lim){
if(lim==1){
poly ret(0);ret[0]=1;
return ret;
}
poly c=(*this);c.resize(lim-1);
int l=(lim+1)/2;
poly F0=c.Sqrt(l);
poly INV=F0.inv(lim-l+1);
auto ret=F0-(F0*F0-c)*INV*((mod+1)>>1);
ret.resize(lim-1);
return ret;
}
};
int32_t main(){
ios::sync_with_stdio(false);
init();
}
|
对于EGF的阐释性的理解。
注意我们说的OGF,EGF都是对于一个数列说的。
\(f\rightarrow F\)。
\(f\) 的EGF定义为 \(\sum \frac{f_i}{i!}\cdot x^i\)。
EGF 的好的性质主要体现在多项式每个系数除了对应的 \(i!\)。
而CMD_BLK 说的二项卷积是对于 \(f,g\)
两个数列的卷积。
对于他们生成的函数/多项式，卷积依旧是经典的加法卷积。
如果考虑原本数列在卷积后变成什么了。
这么理解： \([x^i]1*[x^j]1=\binom{i+j}{i}\),把1变成原先每一项的含义即可。
讲道理prufer序列真的不是特别直观。