「UOJ500」任意基DFT - mem 的小站

「UOJ500」任意基DFT

「UOJ500」任意基DFT

2020 年 3 月 27 日

给定 $n$ 次多项式
$f(x) = \sum_{i=0}^n a_i x^i$
$Q$ 次询问，第 $i$ 次询问 $f(q_i)$ 对 $998244353$ 取模的值。

其中 $q_i$ 是一个一阶线性递推，给定 $q_0, x, y$ ，满足
$q_n = x q_{n-1} + y$
$1 \leq n \leq 2.5 \times 10^5, \ 1 \leq Q \leq 10^6, \ 2 \leq x < 998244353, \ 0 \leq q_0, y < 998244353$ 。

Solution Part1

首先这玩意儿肯定没法两个 $\log$ 多点求值，除非你是钱哥哥。

Solution Part2

虽然这东西好像小学生都会，但我自己推的时候怎么莫名搞到特征多项式那个方向去了，越学越傻逼了。。。

\left\{ \begin{aligned} g_{n+1} = x g_n + y \\ g_n = x g_{n-1} + y \\ \end{aligned} \right.

两式相减得到

(g_{n+1}-g_n) = x (g_n - g_{n-1})

不妨设

h_n = g_{n+1} - g_n

可以得到 $h$ 的通项公式

h_n = x^n h_0 = x^n ((x-1) g_0 + y)

从而得到 $g$ 的通项公式

g_n = g_0 + \sum_{i=0}^{n-1} h_i = g_0 + \frac {x^n-1} {x-1} ((x-1)g_0 + y)

不妨表示成 $ax^n + b$ 的形式，其中

\begin{cases} a = \dfrac {(x-1) g_0 + y} {x-1} \\ b = \dfrac {y} {x-1} \\ \end{cases}

Solution Part3

也就是说我们要对于所有 $i \in [1, Q]$ 求出 $f(a x^i + b)$ 的值，不妨设多项式 $g(x) = f(ax + b)$ ，也就是一个多项式平移。

\begin{aligned} g(x) &= f(ax + b) \\ &= \sum_{i=0}^n a_i (ax + b)^i \\ &= \sum_{i=0}^n a_i \sum_{j=0}^i \frac {i!} {j! (i-j)!} (ax)^j b^{i-j} \\ &= \sum_{i=0}^n a_i i! \sum_{j=0}^i \frac {(ax)^j} {j!} \cdot \frac {b^{i-j}} {(i-j)!} \\ \end{aligned}

容易发现可以卷积维护，需要反转一个多项式。

Solution Part4

回到之前的式子，现在我们只需要对于所以 $i \in [1, Q]$ 求出 $g(x^i)$ 的值，相当于求

\begin{aligned} G(z) &= \sum_{i=0}^Q z^i \sum_{j=0}^n a_j x^{ij} \\ &= \sum_{i=0}^Q z^i \sum_{j=0}^n a_j x^{\binom {i+j} 2 - \binom i 2 - \binom j 2} \\ &= \sum_{i=0}^Q z^i x^{-\binom i2} \sum_{j=0}^n a_j x^{-\binom j2} \cdot x^{\binom {i+j} 2} \\ \end{aligned}

容易发现可以卷积维护，需要反转一个多项式。

题外话：

一开始想过把 $ij$ 拆成 $\binom i 2 + \binom j 2 - \binom {i-j} 2$ ，结果 $i-j$ 的值域有正有负，反而难处理，尽量还是避免这种情况。

Code

#include <bits/stdc++.h>
using namespace std;
const int N = 2.1e6 + 10, mod = 998244353;
int n, q, rev[N];
struct z {
  int x;
  z(int x = 0) : x(x) {}
  friend inline z operator*(z a, z b) { return (long long)a.x * b.x % mod; }
  friend inline z operator-(z a, z b) { return (a.x -= b.x) < 0 ? a.x + mod : a.x; }
  friend inline z operator+(z a, z b) { return (a.x += b.x) >= mod ? a.x - mod : a.x; }
} q0, qx, qy, qc, w[N], f[N], g[N], t1[N], t2[N], t3[N], t4[N], fac[N], ifac[N];
inline z fpow(z a, int b) {
  z s = 1;
  for (; b; b >>= 1, a = a * a)
    if (b & 1) s = s * a;
  return s;
}
void initFac(int n) {
  fac[0] = ifac[0] = ifac[1] = 1;
  for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i;
  for (int i = 2; i <= n; i++) ifac[i] = (mod - mod / i) * ifac[mod % i];
  for (int i = 2; i <= n; i++) ifac[i] = ifac[i - 1] * ifac[i];
}
void initPow(z c, int n, z *a, z *b) {
  a[0] = b[0] = 1;
  for (int i = 1; i < n; i++) a[i] = a[i - 1] * c;
  for (int i = 1; i < n; i++) b[i] = b[i - 1] * a[i - 1];
}
int init(int n) {
  static int wk = 1;
  int k = -1, l = 1;
  while (l < n) l <<= 1, ++k;
  for (int i = 0; i < l; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << k);
  for (int &k = wk; k < l; k <<= 1) {
    z c = fpow(3, (mod - 1) / k / 2);
    w[k] = 1;
    for (int i = 1; i < k; i++) w[i + k] = w[i + k - 1] * c;
  }
  return l;
}
void ntt(z *a, z *b, int l) {
  for (int i = 0; i < l; i++) b[rev[i]] = a[i];
  for (int k = 1; k < l; k <<= 1)
    for (int i = 0; i < l; i += (k << 1))
      for (int j = 0; j < k; j++) {
        z x = b[i + j], y = b[i + j + k] * w[j + k];
        b[i + j] = x + y, b[i + j + k] = x - y;
      }
}
void mul(z *a, z *b, z *c, int n, int m) {
  static z t[N];
  int l = init(n + m - 1);
  z v = fpow(l, mod - 2);
  ntt(a, c, l), ntt(b, t, l);
  for (int i = 0; i < l; i++) t[i] = t[i] * c[i];
  reverse(t + 1, t + l), ntt(t, c, l);
  for (int i = 0; i < l; i++) c[i] = c[i] * v;
}
int main() {
#ifdef memset0
  freopen("1.in", "r", stdin);
#endif
  cin.tie(0)->sync_with_stdio(0);
  cin >> n >> q;
  initFac(max(n, q));
  for (int i = 0; i <= n; i++) cin >> f[i].x;
  cin >> q0.x >> qx.x >> qy.x;
  qc = (q0 * (qx - 1) + qy) * fpow(qx - 1, mod - 2);
  const auto move = [&](z *a, z c, z k, int n) {
    int i;
    z t;
    for (i = 0, t = 1; i < n; i++) t1[i] = a[i] * fac[i], t2[n - i - 1] = ifac[i] * t, t = t * c;
    mul(t1, t2, t3, n, n);
    for (i = 0, t = 1; i < n; i++) a[i] = t3[n - 1 + i] * ifac[i] * t, t = t * k;
  };
  memset(t1, 0, sizeof(t1)), memset(t2, 0, sizeof(t2));
  move(f, q0 - qc, qc, n + 1);
  memset(t1, 0, sizeof(t1)), memset(t2, 0, sizeof(t2));
  const auto czt = [&](z *a, z *b, z c, int n, int q) {
    initPow(c, n + q, t3, t1);
    reverse(t1, t1 + n + q);
    initPow(fpow(c, mod - 2), max(n, q), t3, t4);
    for (int i = 0; i < n; i++) t2[i] = a[i] * t4[i];
    mul(t1, t2, t3, n + q, n);
    for (int i = 0; i < q; i++) b[i] = t3[n + q - 1 - i] * t4[i];
  };
  czt(f, g, qx, n + 1, q + 1);
  cout << accumulate(g + 1, g + q + 1, 0, [](int a, z b) { return a ^ b.x; }) << endl;
}

评论

TABLE OF CONTENTS

Solution Part1

Solution Part2

Solution Part3

Solution Part4

Code