多项式相关算法集成

发表于 | 分类于 | 0 comments
阅读次数:   |   字数统计: 3k(字)   |   阅读时长: 16(分)

原理

快速数论变换(FNT)

原理同FFT一样,用与弧度具有同样性质的一组值替代弧度进行计算,具体性质如下:

FFT 中能在 O(nlog⁡n) 时间内变换用到了单位根 的什么性质:

  1. $\omega_n^n = 1$

  2. , , ⋯, 是互不相同的，这样带入计算出来的点值才可以用来还原出系数

  3. , ，这使得在按照指数奇偶分类之后能够把带入的值也减半使得问题规模能够减半

        这点保证了能够使用相同的方法进行逆变换得到系数表示

参考:

  1. 从多项式乘法到快速傅里叶变换

  2. 快速数论变换 (NTT)

  3. FNT用到的各种素数

多项式逆元

参考:

  1. 多项式求逆元

多项式除法及求模

参考:

  1. 多项式除法及求模

多项式多点求值与拉格朗日插值

参考:

  1. 多项式的多点求值与快速插值
  2. 多项式多点求值和插值
  3. 维基百科:拉格朗日插值法

实现

展开代码 
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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
/*
 * 多项式相关算法模板集合
 */
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
// 2281701377是一个原根为3的指数,平方刚好不会爆long long
// 另外一个常用的998244353 = 119ll * (1 << 23) + 1
// 1004535809=479⋅221+1 加起来刚好不会爆 int 也不错
const ll mod_v = 17ll * (1 << 27) + 1;  
const int MAXL = 18, N = 1 << MAXL, M = 1 << MAXL;  // MAXL最大取mov_v的2次幂
struct polynomial {
​    ll coef[N]; // size设为多项式系数个数的两倍
​    ll a[N], b[N], d[N], r[N];
​    ll xcoor[N], ycoor[N], v[N], mcoef[N];  // size设为点个数的两倍
​    vector<vector<ll> > poly_divisor;
​     
​    static ll mod_pow(ll x,ll n,ll m) {
​        ll ans = 1;
​        while(n > 0){
​    	    if(n & 1)
​    	        ans = ans*x%m;
​    	    x = x*x%m;
​    	    n >>= 1;
​        }
​        return ans;
​    }
​     
​    struct FastNumberTheoreticTransform  {
​        ll omega[N], omegaInverse[N];
​        int range;
​    
​        // 初始化频率
​        void init (const int& n) {
​    	    range = n;
​        	ll base = mod_pow(3, (mod_v - 1) / n, mod_v);
​        	ll inv_base = mod_pow(base, mod_v - 2, mod_v);
​        	omega[0] = omegaInverse[0] = 1;
​        	for (int i = 1; i < n; ++i) {
​        	    omega[i] = omega[i - 1] * base % mod_v;
​        	    omegaInverse[i] = omegaInverse[i - 1] * inv_base % mod_v;
​        	}
​        }
​        
​        // Cooley-Tukey算法:O(n*logn)
​        void transform (ll *a, const ll *omega, const int &n) {
​        	for (int i = 0, j = 0; i < n; ++i) {
​        		    if (i > j) std::swap (a[i], a[j]);
​            		for(int l = n >> 1; ( j ^= l ) < l; l >>= 1);
​            }
​        
​        	for (int l = 2; l <= n; l <<= 1) {
​        	    int m = l / 2;
​        	    for (ll *p = a; p != a + n; p += l) {
​        	        for (int i = 0; i < m; ++i) {
​        	            ll t = omega[range / l * i] * p[m + i] % mod_v;
​        	            p[m + i] = (p[i] - t + mod_v) % mod_v;
​        	            p[i] = (p[i] + t) % mod_v;
​        	        }
​        	    }
​        	}
​        }
​        

        // 时域转频域
        void dft (ll *a, const int& n) {
    	    transform(a, omega, n);
        }
    
        // 频域转时域
        void idft (ll *a, const int& n) {
    	    transform(a, omegaInverse, n);
    	    for (int i = 0; i < n; ++i) a[i] = a[i] * mod_pow(n, mod_v - 2, mod_v) % mod_v;
        }
    } fnt ;
    
    // 与模比较转换值
    ll mod_trans(ll v) {
        return abs(v) <= mod_v / 2 ? v : (v < 0 ? v + mod_v : v - mod_v);
    }
    
    // 二分求\prod{x-xi}的所有二分子多项式的系数
    void binary_subpoly(int l, int r, int idx) {
        if (l == r - 1) {
    	    poly_divisor[idx].push_back(-xcoor[l]);
    	    poly_divisor[idx].push_back(1);        
        } else {
        	int lidx = (idx << 1) + 1, ridx = lidx + 1;
        	binary_subpoly(l, (l + r) / 2, lidx);
        	binary_subpoly((l + r) / 2, r, ridx);
        	int t = poly_divisor[lidx].size() + poly_divisor[ridx].size() - 1;
        	int p = 1;
        	while(p < t) p <<= 1;
        	copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a);
        	fill(a + poly_divisor[lidx].size(), a + p, 0);
        	copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b);
        	fill(b + poly_divisor[ridx].size(), b + p, 0);
        	fnt.dft(a, p);
        	fnt.dft(b, p);
        	for (int i = 0; i < p; i++)  a[i] *= b[i]; 
        	fnt.idft(a, p);
        	for (int i = 0; i < t; i++) poly_divisor[idx].push_back(mod_trans(a[i]));
        }
    }
    
    // 模x^deg,a为要求逆元的多项式系数,结果存放在b[0~deg]中
    // T(deg) = T(deg/2) + deg*log(deg)，复杂度O(deg*log(deg))
    void polynomial_inverse(int deg, ll* a, ll* b, ll* tmp) {
    	if(deg == 1) {
    		b[0] = mod_pow(a[0], mod_v - 2, mod_v);
    	} else {
    		polynomial_inverse((deg + 1) >> 1, a, b, tmp);
    		int p = 1;
    		while(p < (deg << 1) - 1) p <<= 1;
    		copy(a, a + deg, tmp);
    		fill(tmp + deg, tmp + p, 0);
    	    fill(b + ((deg + 1) >> 1), b + p, 0);
    	    //fnt.init(p);
    	    fnt.dft(tmp, p);
    	    fnt.dft(b, p);
    		for(int i = 0; i != p; ++i) {
    			b[i] = (2 - tmp[i] * b[i] % mod_v) * b[i] % mod_v;
    			if(b[i] < 0) b[i] += mod_v;
    		}
    	    fnt.idft(b, p);
    		fill(b + deg, b + p, 0);
    	}
    }
    
    // A = D*B + R,A为n项n-1次幂,B为m项m-1次幂,D为n-m+1项n-m次幂,R为m-1项m-2次幂
    // 要求a,b中系数以低次到多次顺序排列
    // n >= m，复杂度O(n*logn)；n < m，复杂度O(n)
    int polynomial_division(int n, int m, ll *A, ll *B, ll *D, ll *R) {
        if (n < m) {
    	    copy(A, A + n, R);
    	    return n;
        } else {
    	    static ll A0[N], B0[N], tmp[N]; //数组太大会爆栈,添加到全局区
    
    	    int p = 1, t = n - m + 1;
    	    while(p < (t << 1) - 1) p <<= 1;
    
    	    fill(A0, A0 + p, 0);
    	    reverse_copy(B, B + m, A0);
    	    polynomial_inverse(t, A0, B0, tmp);
    	    fill(B0 + t, B0 + p, 0);
    	    fnt.dft(B0, p);
    
    	    reverse_copy(A, A + n, A0);
    	    fill(A0 + t, A0 + p, 0);
    	    fnt.dft(A0, p);
    
    	    for(int i = 0; i != p; ++i)
    	    	A0[i] = A0[i] * B0[i] % mod_v;
    	    fnt.idft(A0, p);
    	    reverse(A0, A0 + t);
    	    copy(A0, A0 + t, D);
    
    	    for(p = 1; p < n; p <<= 1);
    	    fill(A0 + t, A0 + p, 0);
    	    fnt.dft(A0, p);
    	    copy(B, B + m, B0);
    	    fill(B0 + m, B0 + p, 0);
    	    fnt.dft(B0, p);
    	    for(int i = 0; i != p; ++i)
    	    	A0[i] = A0[i] * B0[i] % mod_v;
    	    fnt.idft(A0, p);
    	    for(int i = 0; i != m - 1; ++i)
    	    	R[i] = ((A[i] - A0[i]) % mod_v + mod_v) % mod_v;
    	    //fill(R + m - 1, R + p, 0);
    	    return m - 1;
        }
    }
    
    // 多项式的点值计算
    // l和r为存储要求的点数组的左右边界(左闭右开), idx为除数多项式索引(初始０)
    // polycoef为用于计算点的多项式系数，num为其系数个数
    // 设多项式项数x=num，点数y=r-l，n=max(x,y)，复杂度O(n(logn)^2)
    void polynomial_calculator(int l, int r, int idx, int num, ll *polycoef) {
        int mid = (l + r) / 2;
        int lidx = (idx << 1) + 1, ridx = lidx + 1;
        int lsize = poly_divisor[lidx].size(), rsize = poly_divisor[ridx].size();
        ll *lmod_poly = new ll[lsize - 1], *rmod_poly = new ll[rsize - 1];
        copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a);
        copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b);
        int lplen = polynomial_division(num, lsize, polycoef, a, d, lmod_poly);
        int rplen = polynomial_division(num, rsize, polycoef, b, d, rmod_poly);
        if (l == mid - 1) {
    	    v[l] = mod_trans(lmod_poly[0]);
        } else {
    	    polynomial_calculator(l, (l + r) / 2, lidx, lplen, lmod_poly);
        }
        if (r == mid + 1) {
    	    v[(l + r) / 2] = mod_trans(rmod_poly[0]);
        } else {
    	    polynomial_calculator((l + r) / 2, r, ridx, rplen, rmod_poly);
        }
        delete []lmod_poly;
        delete []rmod_poly;
    }
    
    // 拉格朗日插值：二分＋快速数论变换
    // l和r为存储要求的点数组的左右边界(左闭右开), idx为由点二分构造出多项式的索引(初始０)
    // polycoef为插值得到的多项式结果
    // 设点个数为n，复杂度O(n*(logn)^2)，结果polycoef为n-1次多项式
    void polynomial_interpolate(int l, int r, int idx, ll *polycoef) {
        if (l == r - 1) {
    	    polycoef[0] = ycoor[l] * mod_pow(v[l], mod_v - 2, mod_v) % mod_v;
        } else {
        	int mid = (l + r) >> 1;
        	int lidx = (idx << 1) + 1, ridx = lidx + 1;
        	int sz = poly_divisor[idx].size() - 1;
        	int lsize = poly_divisor[lidx].size(), rsize = poly_divisor[ridx].size();
        	int p = 1;
        	while (p < sz) p <<= 1;
        	ll *leftpoly = new ll[p], *rightpoly = new ll[p];
        	polynomial_interpolate(l, mid, lidx, leftpoly);
        	polynomial_interpolate(mid, r, ridx, rightpoly);
        	copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a);
        	copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b);
        	fill(leftpoly + lsize - 1, leftpoly + p, 0);
        	fill(rightpoly + rsize - 1, rightpoly + p, 0);
        	fill(a + lsize, a + p, 0);
        	fill(b + rsize, b + p, 0);
        	fnt.dft(leftpoly, p);
        	fnt.dft(b, p);
        	for (int i = 0; i < p; i++) leftpoly[i] = leftpoly[i] * b[i] % mod_v;
        	fnt.idft(leftpoly, p);
        	fnt.dft(rightpoly, p);
        	fnt.dft(a, p);
        	for (int i = 0; i < p; i++) rightpoly[i] = rightpoly[i] * a[i] % mod_v;
        	fnt.idft(rightpoly, p);
        	for (int i = 0; i < sz; i++) polycoef[i] = mod_trans((leftpoly[i] + rightpoly[i]) % mod_v); 
        	delete []leftpoly;
        	delete []rightpoly; 
        }
    }
    
    // 初始化点个数到二分子多项式个数
    // 调用polynomial_calculator和polynomial_interpolate前调用
    void init(int vnum) {
        int vnum2 = 1;
        while (vnum2 < vnum) vnum2 <<= 1;
        for (int i = 0; i < 2 * vnum2 - 1; i++)  poly_divisor.push_back(vector<ll>());
    }
} 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
int main() {
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);
    int n, m;
    cin >> n >> m;
    for (int i = 0; i < n; i++) cin >> poly.a[i];
    for (int i = 0; i < m; i++) cin >> poly.b[i];
    cout << "a*b之后的真实系数:" << endl;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
            poly.r[i + j] += poly.a[i] * poly.b[j];
    for (int i = 0; i < n + m - 1; i++)
        cout << poly.r[i] << " ";
    cout << endl;
    cout << "利用fft计算得出的系数:" << endl;
    int p = 1;
    while (p < n + m - 1) p <<= 1;  // 只要p大于多项式结果中的最大次幂即可
    poly.fnt.dft(poly.a, p);
    poly.fnt.dft(poly.b, p);
    for (int i = 0; i < p; i++) {
        poly.d[i] = poly.a[i] * poly.b[i] % mod_v;
    }
    poly.fnt.idft(poly.d, p);
    for (int i = 0; i < n + m - 1; i++) cout << poly.d[i] << " ";
	return 0;
}

多项式逆元

展开代码 
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
int main() {
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);
    int n;
    cin >> n;
    for(int i = 0; i != n; ++i)
        cin >> poly.a[i];
    int m;
    cin >> m;  // 输入模x^m
    poly.polynomial_inverse(m, poly.a, poly.b, poly.d);
    cout << "inverse: " << endl;
    for(int i = 0; i != m; ++i)
        cout << (poly.b[i] + mod_v) % mod_v << " ";
    cout << endl;
    cout << "a*b相乘取模验证取模后的前m项系数:" << endl;
    memset(poly.d, 0, sizeof(poly.d));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            poly.d[i + j] = (poly.d[i + j] + poly.a[i] * poly.b[j] % mod_v)% mod_v;
        }   
    }   
    cout << m << endl;
    for (int i = 0; i < m; i++) {
        cout << poly.d[i] << " ";
    }   
    return 0;
}

输入数据得到输出后的结果如下图

多项式除法以及取模

展开代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
int main() {
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i != n; ++i)
    cin >> poly.a[i]; // 0次幂系数开始输入,缺失的幂系数输入0
    for (int i = 0; i < m; i++)
    cin >> poly.b[i]; // 0次幂系数开始输入
    int rlen = poly.polynomial_division(n, m, poly.a, poly.b, poly.d, poly.r);
    cout << "输出商多项式:" << endl;
    for (int i = 0; i < n - m + 1; i++) 
    cout << poly.d[i] << " ";
    cout << endl;
    cout << "输出余数多项式:" << endl;
    for (int i = 0; i < rlen; i++)
    cout << poly.r[i] << " ";
	return 0;
}

多项式多点求值

展开代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
int main() {
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);
    int n, vnum; 
    // 输入要计算的点
    cin >> n;
    for (int i = 0; i < n; i++) {
        cin >> poly.coef[i];
    }
    cin >> vnum;
    for (int i = 0; i < vnum; i++) {
        cin >> poly.xcoor[i];        
    }
    // 二分求子多项式                                                                               
    poly.init(vnum);
    poly.binary_subpoly(0, vnum / 2, 1);
    poly.binary_subpoly(vnum / 2, vnum, 2); 
    // 将输入点代入插值得到的多项式中进行验证，输出计算得到的结果
    cout << "计算结果:" << endl;
    poly.polynomial_calculator(0, vnum, 0, n, poly.coef);
    for (int i = 0; i < vnum; i++) cout << poly.v[i] << " ";
    cout << endl;
}

输入数据后得到结果如下

拉格朗日快速插值计算

展开代码 
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
int main() {
    // 多项式系数默认低次到高次排列
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);   
    // 多项式插值
    int vnum; 
    // 输入要计算的点
    cin >> vnum;
    for (int i = 0; i < vnum; i++) {
        cin >> poly.xcoor[i] >> poly.ycoor[i];        
    }
    // 二分求子多项式
    poly.init(vnum);
    poly.binary_subpoly(0, vnum, 0); 
    for (unsigned int i = 1; i < poly.poly_divisor[0].size(); i++) {
        poly.mcoef[i - 1] = poly.poly_divisor[0][i] * i % mod_v;
    }
    // 遍历i计算所有\sum_{j!=i}{xi-xj} 
    poly.polynomial_calculator(0, vnum, 0, poly.poly_divisor[0].size() - 1, poly.mcoef);
    // 拉格朗日插值计算多项式
    poly.polynomial_interpolate(0, vnum, 0, poly.coef);
    // 输出插值得到的多项式系数 
    for (int i = 0; i < vnum; i++) {
        cout << poly.coef[i] << " ";
    }
    cout << endl;
    // 将输入点代入插值得到的多项式中进行验证，输出计算得到的结果
    poly.polynomial_calculator(0, vnum, 0, vnum, poly.coef);
    for (int i = 0; i < vnum; i++) cout << poly.v[i] << " ";
    cout << endl;
    return 0;
}

输入数据得到的结果如下

---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址: https://brooksj.com/2019/07/24/多项式相关算法集成/
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢！