标签:fft swap ios oid stat 要求 除法 algorithm cst
http://blog.miskcoo.com/2015/05/polynomial-division
给出一个 \(n\) 次多项式 \(A(x)\) 和一个 \(m(m \le n)\) 次多项式 \(B(x)\) ,要求求出两个多项式 \(D(x), R(x)\) , 满足
其中 \(degD \le degA - degB = n - m, degR < m\)
可以在 \(O(n \log n)\) 的时间求解
首先,我们先想办法消除 \(R(x)\) 的影响,我们定义
实际上就是将 \(A(x)\) 的系数翻转, 例如
接下来,我们将 \((1)\) 中的 \(x\) 用 \(\dfrac 1x\) 替换,并在两边同乘 \(x^n\)
观察发现此时 \(x^{n - m + 1}R^R(x)\) 的非零项都在 \(n - m + 1\) 上,而 \(D^R(x)\) 的最高次项为 \(n - m\) , 那么我们有
那么我们就可以用一次求出逆元的复杂度求出 \(D(x)\) ,再代回原式就可以得到 \(R(x)\)
不打换行真的会死.......
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
inline char nc() {
static char buf[100000], *l = buf, *r = buf;
return l==r&&(r=(l=buf)+fread(buf,1,100000,stdin),l==r)?EOF:*l++;
}
template<class T> void read(T &x) {
x = 0; int f = 1, ch = nc();
while(ch<‘0‘||ch>‘9‘){if(ch==‘-‘)f=-1;ch=nc();}
while(ch>=‘0‘&&ch<=‘9‘){x=x*10-‘0‘+ch;ch=nc();}
x *= f;
}
#define inv(a) quick_power(a, mod - 2)
typedef long long ll;
const int g = 3;
const int mod = 998244353;
const int phi = mod - 1;
const int maxn = 400000 + 5;
const int maxlog = 20;
int n, m; ll F[maxn], G[maxn], D[maxn], R[maxn];
inline ll sum(ll x) {
return x >= mod ? x - mod : x;
}
inline ll dec(ll x) {
return x < 0 ? x + mod : x;
}
ll quick_power(ll x, ll y) {
ll re = 1;
while(y) {
if(y & 1) re = re * x % mod;
x = x * x % mod;
y >>= 1;
}
return re;
}
namespace polynomial {
int rev[maxn];
ll w[maxlog][maxn][2];
void init() {
for(int i = 1, s = 0; i < maxn; i <<= 1, ++s) {
ll wn0 = quick_power(g, phi / (i * 2));
ll wn1 = quick_power(g, -phi / (i * 2) + phi);
w[s][0][0] = w[s][0][1] = 1;
for(int k = 1; k < i; ++k) {
w[s][k][0] = w[s][k - 1][0] * wn0 % mod;
w[s][k][1] = w[s][k - 1][1] * wn1 % mod;
}
}
}
void init_rev(int n, int L) {
for(int i = 1; i < n; ++i) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (L - 1));
}
}
void FFT(ll *A, int n, int f) {
int d = f == -1;
for(int i = 0; i < n; ++i) if(i > rev[i]) {
swap(A[i], A[rev[i]]);
}
for(int i = 1, s = 0; i < n; i <<= 1, ++s) {
for(int j = 0, p = i << 1; j < n; j += p) {
ll *u = A + j, *v = A + j + i;
for(int k = 0; k < i; ++k, ++u, ++v) {
ll x = *u, y = *v * w[s][k][d] % mod;
*u = sum(x + y);
*v = dec(x - y);
}
}
}
if(f == -1) {
ll r = inv(n);
for(int i = 0; i < n; ++i) {
A[i] = A[i] * r % mod;
}
}
}
void inverse(int step, ll *A, ll *B) {
static ll T[maxn];
if(step == 1) {
B[0] = inv(A[0]);
return;
}
inverse((step + 1) >> 1, A, B);
int n = 1, L = 0, m = step << 1;
for(n = 1; n <= m; n <<= 1) ++L;
init_rev(n, L);
copy(A, A + step, T), fill(T + step, T + n, 0);
FFT(T, n, 1);
FFT(B, n, 1);
for(int i = 0; i < n; ++i) {
B[i] = dec(2 - T[i] * B[i] % mod) * B[i] % mod;
}
FFT(B, n, -1);
fill(B + step, B + n, 0);
}
void division(ll *A, ll *B, int n, int m, ll *D, ll *R) {
static ll A0[maxn], B0[maxn];
memset(A0, 0, sizeof(A0));
memset(B0, 0, sizeof(B0));
int len = n - m + 1;
reverse_copy(B, B + m + 1, A0);
inverse(len, A0, B0);
int p, L = 0, tmp = len << 1;
for(p = 1; p <= tmp; p <<= 1) ++L;
init_rev(p, L);
reverse_copy(A, A + n + 1, A0), fill(A0 + len, A0 + p, 0);
FFT(A0, p, 1);
FFT(B0, p, 1);
for(int i = 0; i < p; ++i) {
A0[i] = A0[i] * B0[i] % mod;
}
FFT(A0, p, -1);
reverse(A0, A0 + len);
copy(A0, A0 + len, D);
L = 0, tmp = n;
for(p = 1; p <= tmp; p <<= 1) ++L;
init_rev(p, L);
copy(B, B + m + 1, B0), fill(B0 + m + 1, B0 + p, 0);
fill(A0 + len, A0 + p, 0);
FFT(A0, p, 1);
FFT(B0, p, 1);
for(int i = 0; i < p; ++i) {
A0[i] = A0[i] * B0[i] % mod;
}
FFT(A0, p, -1);
for(int i = 0; i < m; ++i) {
R[i] = dec(A[i] - A0[i]);
}
}
}
int main() {
// freopen("testdata.in", "r", stdin);
read(n), read(m);
for(int i = 0; i <= n; ++i) read(F[i]);
for(int i = 0; i <= m; ++i) read(G[i]);
polynomial :: init();
polynomial :: division(F, G, n, m, D, R);
for(int i = 0; i <= n - m; ++i) {
if(i) printf(" ");
printf("%lld", D[i]);
}
printf("\n");
for(int i = 0; i < m; ++i) {
if(i) printf(" ");
printf("%lld", R[i]);
}
printf("\n");
return 0;
}
标签:fft swap ios oid stat 要求 除法 algorithm cst
原文地址:https://www.cnblogs.com/ljzalc1022/p/12909522.html