多项式 - 除法与取模

Posted 2021-02-08 nanjoqin

一类问题：给定一个 (n) 次多项式 (F(x)) 和一个 (m) 次多项式 (G(x))，请求出多项式 (Q(x))，(R(x))，满足以下条件：

• (Q(x)) 次数为 (n?m)，(R(x)) 次数小于 (m)
• (F(x)=Q(x)?G(x)+R(x))

所有的运算在模 (998244353) 意义下进行。
　

多项式求逆

放一篇博客。直接做多项式的除法是无从下手的，或者说没有很优秀的方法可以优化它的复杂度。

然而逆元是个很方便的东西，乘法和除法是一对互逆运算，于是想到可以找到多项式 (G(x)) 的逆元 (G^{-1}(x)) 而转变成 (NTT) 解决的多项式乘问题。

我们设 (B(x)) 是 (A(x)) 在模 (x^n) 意义下的逆元，有 (A(x)B(x) equiv 1 (mod x^n))，下面列出几条规律。

当 (n=1) 时，(A(x)=c)，(c) 是一个常数，则 (B(x)=c^{-1})。

当 (n>1) 时，假设在 ((mod x^{lceil frac{n}{2} ceil)}) 意义下 (A(x)) 的逆元是 (B^′(x))，并且我们已经求出，那么：(A(x)B^′(x) equiv 1(mod x^{lceil frac{n}{2} ceil}))

而且在 (mod x^{lceil frac{n}{2} ceil}) 下也有：(A(x)B(x) equiv 1(mod x^{lceil frac{n}{2} ceil}))

然后就可以得到：(B(x)?B^′(x) equiv 0(mod x^{lceil frac{n}{2} ceil}))

两边平方 (B^2(x)?2B^′(x)B(x)+B^{′2}(x) equiv 0(mod x^n))

然后同时乘上 (A(x))，移项可以得到 (B(x) equiv 2B^′(x)?A(x)B^{′2}(x)(mod x^n))，由此式即可得到逆元。
　

(Code：)

#include <cmath>
#include <queue>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long u64;

const int mod = 998244353;
const int maxn = 600000 + 10;
int n, r[maxn], bit[maxn]; u64 f[maxn], inv[maxn], tmp[maxn];

inline int read() {
  register char ch = 0; register int w = 0, x = 0;
  while( !isdigit(ch) ) w |= (ch == ‘-‘), ch = getchar();
  while( isdigit(ch) ) x = (x * 10) + (ch ^ 48), ch = getchar();
  return w ? -x : x;
}

inline u64 Fast_pow(u64 a, int p) {
  u64 x = a, ans = 1ll;
  for( ; p; x = x * x % mod, p = p >> 1) if( p & 1 ) ans = x * ans % mod;
  return ans;
}

inline void Number_theory_transform(u64 *a, int limit, int type) {
  int rado = bit[limit];
  for(int i = 0; i < limit; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (rado - 1));
  for(int i = 0; i < limit; ++i) if( i < r[i] ) swap(a[i], a[r[i]]);
  for(int mid = 1; mid < limit; mid = mid << 1) {
    u64 Base_p = Fast_pow(3ll, (mod - 1) / (mid << 1));
    if( type == -1 ) Base_p = Fast_pow(Base_p, mod - 2);
    for(int l = 0, length = mid << 1; l < limit; l = l + length) {
      for(int k = 0, p = 1; k < mid; ++k, p = p * Base_p % mod) {
        u64 x = a[l + k], y = p * a[l + mid + k] % mod;
        a[l + k] = (x + y) % mod, a[l + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if( type == -1 ) for(int i = 0; i < limit; ++i) a[i] = a[i] * Fast_pow(limit, mod - 2) % mod;
}

inline void Inverse_transform(u64 *a, u64 *b, int limit) {
  if( limit == 1 ) { b[0] = Fast_pow(a[0], mod - 2); return ; }
  Inverse_transform(a, b, (limit + 1) >> 1);
  int len = limit << 1;
  for(int i = 0; i < limit; ++i) tmp[i] = a[i];
  for(int i = limit; i < len; ++i) tmp[i] = 0;
  Number_theory_transform(tmp, len, 1);
  Number_theory_transform(b, len, 1);
  for(int i = 0; i < len; ++i) b[i] = (2 * b[i] % mod - (b[i] * b[i] % mod) * tmp[i] % mod + mod) % mod;
  Number_theory_transform(b, len, -1);
  for(int i = limit; i < len; ++i) b[i] = 0;
}

int main(int argc, char const *argv[])
{
  freopen("..\nanjolno.in", "r", stdin);
  freopen("..\nanjolno.out", "w", stdout);

  scanf("%d", &n);
  for(int i = 0; i < n; ++i) f[i] = read();
  int limit = 1, rado = 0;
  while( limit <= n ) limit = limit << 1, ++rado;
  for(int i = 0; i < 20; ++i) bit[1 << i] = i;
  Inverse_transform(f, inv, limit);
  for(int i = 0; i < n; ++i) printf("%lld ", inv[i]);

  fclose(stdin), fclose(stdout);
  return 0;
}

　

多项式除法与取模

有了逆元，事情就简单了很多，但是依然只是个开始。

考虑最初的式子 (F(x)=Q(x)?G(x)+R(x))，要求得 (Q(x)) 最大的障碍就是 (R(x))，于是就想把 (R(x)) 通过一定的方法试试能不能消掉。

把上式的 (x) 全部用 (frac{1}{x}) 来替代，然后等式两边同时乘上 (x^n)，这意味着对于一个多项式其系数会全部反转，等式变成了：

(F^R(x)=Q^R(x)G^R(x)+x^{n?m+1}R^R(x))

这样，(Q(x)) 是要求的元素之一，反转后次数仍然不高于 (n?m)，而 (x^{n?m+1}R(x)) 这个部分的最低次项高于 (n?m)，放到 (mod x^{n?m+1}) 意义下，(R(x)) 的影响被消除了。

于是 (F^R(x)=Q^R(x)G^R(x)(mod x^{n?m+1}))。

这样就只需要求一个 (G^R(x)) 的逆元，就可以利用多项式常用的倍增方法得到 (Q^R(x))，那么显然，(R(x)=F(x)-G(x)Q(x))。
　

(Code：)

inline int Limit(int len) { for(int i = 1; ; i = i << 1) if( i > len ) return i; }

inline void Multiply_transform(u64 *a, u64 *b, u64 *c, int len_1, int len_2) {
  int len = Limit(len_1 + len_2);
  for(int i = 0; i <= len_1; ++i) tmp_1[i] = a[i];
  for(int i = 0; i <= len_2; ++i) tmp_2[i] = b[i];
  Number_theory_transform(tmp_1, len, 1);
  Number_theory_transform(tmp_2, len, 1);
  for(int i = 0; i < len; ++i) c[i] = tmp_1[i] * tmp_2[i] % mod, tmp_1[i] = tmp_2[i] = 0;
  Number_theory_transform(c, len, -1);
}

int main(int argc, char const *argv[])
{
  for(int i = 0; i < 20; ++i) bit[1 << i] = i;
  scanf("%d%d", &n, &m);
  for(int i = 0; i <= n; ++i) f_1[i] = read();
  for(int i = 0; i <= m; ++i) f_2[i] = read();
  reverse(f_1, f_1 + n + 1), reverse(f_2, f_2 + m + 1);
  Inverse_transform(f_2, inv, Limit(n - m));
  Multiply_transform(f_1, inv, ans, n - m, n - m);
  reverse(ans, ans + n - m + 1);
  for(int i = 0; i <= n - m; ++i) printf("%lld ", ans[i]); printf("
");
  reverse(f_1, f_1 + n + 1), reverse(f_2, f_2 + m + 1);
  Multiply_transform(ans, f_2, mul, n - m, m);
  for(int i = 0; i < m; ++i) printf("%lld ", (f_1[i] - mul[i] + mod) % mod); printf("
");

  return 0;
}

　
　　　　　　　　　　　　　　　　　　　　　与君相遇，乃思长生。