Polynomial Division
コードについての説明
体 $K$ 上の多項式 $f, g$ ($\mathrm{deg}(f) = n$, $\mathrm{deg}(g) = m$, ($g$ の最高次の係数) $\neq 0$) を考える.
このとき $f = q \cdot g + r$ を満たす多項式 $q, r$ で $\mathrm{deg}(r) < m$ となるものを求めるアルゴリズムである.
体上の多項式環はユークリッド整域であることからそのような $q, r$ はかならず存在する.
可換環上の多項式についても同様の division を考えることができるがその場合は $g$ がモニック多項式である場合に well-defined となる.
以下の実装は $K = \mathbb{Z} / p \mathbb{Z}$ ($p$ は素数) の場合の実装で特に素数 $p$ について
$p = k * 2^l + 1 (k > 0, n + 1 \le 2^l)$ を満たすこと(数論変換に乗ること)を仮定している.
こちらの参考資料が分かりやすく, それを基に実装を行った.
以下では $n < m$ の場合は明らかなので $n \ge m$ の場合について示す.
$f(x) = q(x) \cdot g(x) + r(x)$ を満たす $q(x), r(x)$ を求めたいのだが, 同値変形すると
$x^n f(1/x) = \left( x^{n-m} q(1/x) \right) \cdot \left(x^m g(1 / x) \right) + x^{n-m+1} \cdot \left( x^{m-1} r(1/x) \right)$ となる.
今多項式 $h(x)$ について $h'_k(x)$ を $h'_k(x) = x^k h(1/x) (k \ge \mathrm{deg}(h))$ と定義すると, 上式は
$f'_n(x) = q'_{n-m}(x) \cdot g'_m(x) + x^{n-m+1} \cdot r'_{m-1}(x)$ と書き直せる. さらに,
$\Rightarrow f'_n(x) \equiv q'_{n-m}(x) \cdot g'_m(x)\ \mathrm{mod}\ x^{n-m+1}$
$\Rightarrow q'_{n-m}(x) \equiv f'_{n}(x) \cdot g'_m(x)^{-1}\ \mathrm{mod}\ x^{n-m+1}$ と変形することができる.
ここで $g'_m(x)^{-1}$ は $K[x] / (x^{n-m+1})$ 上の $g'_m(x)$ の乗法逆元を意味する.
つまり $ g'_m (x) \cdot g'_m(x)^{-1} \equiv 1\ \mathrm{mod}\ x^{n-m+1}$ を満たすような多項式とする.
これはこちらのアルゴリズムを用いることで $\O (n \log n)$ 時間で求めることができる.
最後に $\mathrm{deg} (q'_{n-m}) \le n - m$ であることからゆえに
$q'_{n-m}(x) = ( \left( f'_n(x) \cdot h_l(x) \right) \ \mathrm{mod}\ x^{n-m+1})$ となり, $q$ を求めることができる.
加えて $r$ も $r = f - q \cdot g$ より計算することができる.
時間計算量だが 逆元を求める部分にかかる計算量が $\O( n \log n)$ で, その他の部分の計算も $\O( n \log n)$ で行うことができるので全体で $\O( n \log n)$ となる.
(注) 便宜上長さ $n - m + 1$ の配列を解となる商の多項式として返しているが必ずしも $n - m$ 次多項式とは限らないことに注意する(最高次の係数が $\mathbb{Z} / p \mathbb{Z}$ 上で $0$ となる可能性があるので).
剰余の多項式の方は $0$ でない限り最高次係数が非ゼロとなるようになっている.
(関数)
polynomial_division$(a, b)$: それぞれ $a, b$ を係数に持つ多項式 $f, g$ について $f = q \cdot g + r$ を満たす多項式 $q, r$ で $\mathrm{deg}(r) < m$
を満たす多項式 $q, r$ の係数を ($q$ の係数, $r$ の係数) のように配列のペアとして返す.
時間計算量: $\O (n \log n)$
コード
#define MOD 998244353
#define root 3
unsigned int add(const unsigned int x, const unsigned int y)
{
return (x + y < MOD) ? x + y : x + y - MOD;
}
unsigned int sub(const unsigned int x, const unsigned int y)
{
return (x >= y) ? (x - y) : (MOD - y + x);
}
unsigned int mul(const unsigned int x, const unsigned int y)
{
return (unsigned long long)x * y % MOD;
}
unsigned int mod_pow(unsigned int x, unsigned int n)
{
unsigned int res = 1;
while(n > 0){
if(n & 1){ res = mul(res, x); }
x = mul(x, x);
n >>= 1;
}
return res;
}
unsigned int inverse(const unsigned int x)
{
return mod_pow(x, MOD - 2);
}
void ntt(vector<int>& a, const bool rev = false)
{
unsigned int i, j, k, l, p, q, r, s;
const unsigned int size = a.size();
if(size == 1) return;
vector<int> b(size);
r = rev ? (MOD - 1 - (MOD - 1) / size) : (MOD - 1) / size;
s = mod_pow(root, r);
vector<unsigned int> kp(size / 2 + 1, 1);
for(i = 0; i < size / 2; ++i) kp[i + 1] = mul(kp[i], s);
for(i = 1, l = size / 2; i < size; i <<= 1, l >>= 1){
for(j = 0, r = 0; j < l; ++j, r += i){
for(k = 0, s = kp[i * j]; k < i; ++k){
p = a[k + r], q = a[k + r + size / 2];
b[k + 2 * r] = add(p, q);
b[k + 2 * r + i] = mul(sub(p, q), s);
}
}
swap(a, b);
}
if(rev){
s = inverse(size);
for(i = 0; i < size; i++){ a[i] = mul(a[i], s); }
}
}
vector<int> convolute(const vector<int>& a, const vector<int>& b, int asize, int bsize, int _size)
{
const int size = asize + bsize - 1;
int t = 1;
while (t < size){ t <<= 1; }
vector<int> A(t, 0), B(t, 0);
for(int i = 0; i < asize; i++){ A[i] = (a[i] < MOD) ? a[i] : (a[i] % MOD); }
for(int i = 0; i < bsize; i++){ B[i] = (b[i] < MOD) ? b[i] : (b[i] % MOD); }
ntt(A), ntt(B);
for(int i = 0; i < t; i++){ A[i] = mul(A[i], B[i]); }
ntt(A, true);
A.resize(_size);
return A;
}
vector<int> polynomial_inverse(const vector<int>& a, int r)
{
vector<int> h = {(int)inverse(a[0])};
int t = 1;
for(int i = 0; t < r; ++i){
t <<= 1;
vector<int> res = convolute(a, convolute(h, h, t / 2, t / 2, t), min((int)a.size(), t), t, t);
for(int j = 0; j < t; ++j){
res[j] = sub(0, res[j]);
if(j < t / 2) res[j] = add(res[j], mul(2, h[j]));
}
swap(h, res);
}
h.resize(r);
return h;
}
pair<vector<int>, vector<int> > polynomial_division(const vector<int>& a, const vector<int>& b)
{
const int n = a.size() - 1, m = b.size() - 1;
assert(b[m] != 0);
if(n < m) return {vector<int>(n - m + 1, 0), a};
vector<int> reva(n + 1), revb(m + 1);
for(int i = 0; i <= n; ++i) reva[n - i] = a[i];
for(int i = 0; i <= m; ++i) revb[m - i] = b[i];
vector<int> inv_revb = polynomial_inverse(revb, n - m + 1);
vector<int> res = convolute(reva, inv_revb, n - m + 1, n - m + 1, n - m + 1);
vector<int> q(n - m + 1), r;
for(int i = 0; i < n - m + 1; ++i) q[i] = res[n - m - i];
vector<int> qb = convolute(q, b, n - m + 1, m + 1, n + 1);
bool first = false;
for(int i = n; i >= 0; --i){
const int val = sub(a[i], qb[i]);
if(!first && val > 0){
first = true, r.resize(i + 1);
}
if(first) r[i] = val;
}
return {q, r};
}