Polynomial

12 / 01 / 2020 | 最后修改于 12 / 01 / 2020

OI Math

多项式乘法 NTT

前置知识

原根

若存在 gg 且使得 gk≡1(modp)g^k \equiv 1 \pmod p 的最小正整数为 φ(p)\varphi(p)，则 gg 为 pp 的原根

以下讨论中选择 pp 为质数

nn 次单位根 wn=gp−1nw_n = g^{\frac {p-1} n}

一些性质

主过程

f=[A0,A1,A2,⋯ ,An−1]f(wni)=A0+A1wni+A2wni2+⋯+An−1wni(n−1)f1=[A0,A2,A4,⋯ ,An−2]f2=[A1,A3,A5,⋯ ,An−1]f(wni)=f1(wni2)+wnif2(wni2)f(wni)=f1(wn2i)+wnif2(wn2i)f(wni+n2)=f1(wn2i)−wnif2(wn2i)\begin{aligned} f &= [A_0, A_1, A_2, \cdots, A_{n-1}] \\ f(w_n^i) &= A_0 + A_1w_n^i + A_2w_n^{i2} + \cdots + A_{n-1}w_n^{i(n-1)} \\ f_1 &= [A_0, A_2, A_4, \cdots, A_{n-2}] \\ f_2 &= [A_1, A_3, A_5, \cdots, A_{n-1}] \\ f(w_n^i) &= f_1(w_n^{i2}) + w_n^if_2(w_n^{i2}) \\ f(w_n^i) &= f_1(w_{\frac{n}{2}}^{i}) + w_n^if_2(w_{\frac{n}{2}}^{i}) \\ f(w_n^{i+\frac{n}{2}}) &= f_1(w_{\frac{n}{2}}^{i}) - w_n^if_2(w_{\frac{n}{2}}^{i}) \end{aligned}

多项式乘法逆 Inv

倍增处理，假设已经求出 (modxn2)\pmod {x^\frac{n}{2}} 意义下的逆元记为 inv0inv_0

f∗inv0≡1(modxn2)f∗inv≡1(modxn)f∗(inv−inv0)≡0(modxn2)inv−inv0≡0(modxn2)\begin{aligned} f * inv_0 &\equiv 1 \pmod {x^\frac{n}{2}} \\ f * inv &\equiv 1 \pmod {x^n} \\ f * (inv - inv_0) &\equiv 0 \pmod {x^\frac{n}{2}} \\ inv - inv_0 &\equiv 0 \pmod {x^\frac{n}{2}} \end{aligned}

因为 inv−inv0inv - inv_0 最低系数不为零的项次数为 xn2x^\frac{n}{2} 平方后 最低系数不为零的项次数为 xnx^n 所以

(inv−inv0)2≡0(modxn)inv2+inv02−2inv∗inv0≡0(modxn)\begin{aligned} (inv - inv_0)^2 &\equiv 0 \pmod {x^n} \\ inv^2 + inv_0^2 - 2inv * inv_0 &\equiv 0 \pmod {x^n} \end{aligned}

两边同时乘 ff

inv+inv02∗f−2inv0≡0(modxn)inv≡2inv0−inv02∗f(modxn)\begin{aligned} inv + inv_0^2*f - 2inv_0 &\equiv 0 \pmod {x^n} \\ inv &\equiv 2inv_0 - inv_0^2*f \pmod {x^n} \end{aligned}

多项式对数函数 Ln

先微分再积分

g′=(ln⁡f)′=f′fg' = (\ln f)' = \frac{f'}{f}

多项式指数函数 Exp

g=efg = e^f \\

设 h(g)=ln⁡g−fh(g) = \ln g - f 其中 gg 为多项式

求 h(g)≡0(modxn)h(g) \equiv 0 \pmod {x^n} 的根，考虑牛顿迭代法

假设已经求出 h(g0)≡0(modxn2)h(g_0) \equiv 0 \pmod {x^\frac{n}{2}}

h(g)=h(g0)+h′(g0)1!(g−g0)+h′′(g0)2!(g−g0)2+⋯h(g) = h(g_0) + \frac{h'(g_0)}{1!}(g-g_0) + \frac{h''(g_0)}{2!}(g-g_0)^2 + \cdots

因为 g−g0g - g_0 最低系数不为零的项次数为 xn2x^\frac{n}{2} 平方后 最低系数不为零的项次数为 xnx^n 所以

h(g)≡h(g0)+h′(g0)∗(g−g0)(modxn)0≡h(g0)+h′(g0)∗(g−g0)(modxn)g≡g0−h(g0)h′(g0)(modxn)g≡g0−g0∗(ln⁡g0−f)(modxn)g≡g0∗(f−ln⁡g0+1)(modxn)\begin{aligned} h(g) &\equiv h(g_0) + h'(g_0) * (g-g_0) \pmod {x^n} \\ 0 &\equiv h(g_0) + h'(g_0) * (g-g_0) \pmod {x^n} \\ g &\equiv g_0 - \frac{h(g_0)}{h'(g_0)} \pmod {x^n} \\ g &\equiv g_0 - g_0 * (\ln g_0 - f) \pmod {x^n} \\ g &\equiv g_0 * (f - \ln g_0 + 1) \pmod {x^n} \end{aligned}

多项式快速幂 Pow

ans=fkln⁡ans=ln⁡fkln⁡ans=kln⁡fans=ekln⁡f (ignore the constant)\begin{aligned} ans &= f^k \\ \ln ans &= \ln f^k \\ \ln ans &= k\ln f \\ ans &= e^{k\ln f} \ \text{(ignore the constant)} \end{aligned} 

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
/// @tags: NTT Polynomial
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>

typedef long long ll;

int const N = 1 << 18, P = 998244353, g = 3, ig = 332748118, inv2 = 499122177;

template <typename T>
inline T &read(T &x) {
  x = 0;
  bool F = false;
  short ch = getchar();
  while (!isdigit(ch)) {
    if (ch == '-') F = true;
    ch = getchar();
  }
  while (isdigit(ch)) x = x * 10ll % P + (ch ^ '0'), ch = getchar();
  if (F) x = -x;
  return x;
}

template <typename T>
inline T &readmod(T &x) {
  x = 0;
  bool F = false;
  short ch = getchar();
  while (!isdigit(ch)) {
    if (ch == '-') F = true;
    ch = getchar();
  }
  while (isdigit(ch)) x = 1ll * x * 10 % P + (ch ^ '0'), ch = getchar();
  if (F) x = -x;
  return x;
}

class Polynomial {
 private:
  int F[N];
  static int Cvt[N << 2];

 public:
  void NTT(bool const typ, int const n);
  void inv(Polynomial &Res, int const n) const;
  void ln(Polynomial &Res, int const n) const;
  void exp(Polynomial &Res, int const n) const;
  void pow(Polynomial &Res, int const n, int const k) const;
  void sqrt(Polynomial &res, int const n) const;
  void div(Polynomial &Q, Polynomial &R, Polynomial const &D, int const n,
           int const m) const;
  int *operator&() { return F; }
  int &operator[](int index) { return F[index]; };
  int const &operator[](int index) const { return F[index]; };
  inline void pre(int n) {
    for (int i = 1, len = 2; len <= n; ++i, len <<= 1)
      for (int j = 1, *const cvt = Cvt + len - 1; j < len; ++j)
        cvt[j] = cvt[j >> 1] >> 1 | ((j & 1) << (i - 1));
  }
  inline void clear(int n) {
    int maxl = 1;
    while (maxl < n) maxl <<= 1;
    memset(F, 0, sizeof(int) * maxl);
  }
} F, G;

int n, m;
int Polynomial::Cvt[N << 2];

inline ll qpow(ll base, int exp) {
  ll Res = 1;
  while (exp) {
    if (exp & 1) Res = Res * base % P;
    base = base * base % P;
    exp >>= 1;
  }
  return Res;
}

/// @param typ 正/逆向 @param n 项数 必须是2的整数次幂
inline void Polynomial::NTT(bool const typ, int const n) {
  for (int i = 1, *const cvt = Cvt + n - 1; i < n; ++i)
    if (i < cvt[i]) std::swap(F[i], F[cvt[i]]);
  for (int i = 2; i <= n; i <<= 1) {
    int mid = i >> 1, wn = qpow(typ ? g : ig, (P - 1) / i);
    for (int j = 0; j < n; j += i) {
      ll wk = 1;
      for (int k = 0; k < mid; ++k, (wk *= wn) %= P) {
        ll t = wk * F[j + k + mid] % P;
        if ((F[j + k + mid] = F[j + k] - t) < 0) F[j + k + mid] += P;
        if ((F[j + k] += t) >= P) F[j + k] -= P;
      }
    }
  }
  if (!typ) {
    ll inv = qpow(n, P - 2);
    for (int i = 0; i < n; ++i) F[i] = inv * F[i] % P;
  }
}

/// @param Res 应清空 @param n 模的次数（项数）
inline void Polynomial::inv(Polynomial &Res, int const n) const {
  static Polynomial tmp;
  if (n == 1) return Res[0] = qpow(F[0], P - 2), void();
  inv(Res, (n + 1) >> 1);
  int maxl = 1;
  while (maxl < n << 1) maxl <<= 1;  // n - 1 次多项式卷 ⌈n / 2⌉ - 1 次多项式
  tmp.clear(n << 1);
  memcpy(&tmp, F, sizeof(int) * n);
  tmp.NTT(true, maxl), Res.NTT(true, maxl);
  for (int i = 0; i < maxl; ++i)
    Res[i] = static_cast<ll>(Res[i]) *
             ((2ll - static_cast<ll>(Res[i]) * tmp[i] % P + P) % P) % P;
  Res.NTT(false, maxl);
  for (int i = n; i < maxl; ++i) Res[i] = 0;  // mod x^n
}

inline void Polynomial::ln(Polynomial &Res, int const n) const {
  static Polynomial df, invf;
  df.clear(n << 1), invf.clear(n << 1);
  int maxl = 1;
  while (maxl < n << 1) maxl <<= 1;
  for (int i = 0; i + 1 < n; ++i) df[i] = static_cast<ll>(F[i + 1]) * (i + 1) % P;
  inv(invf, n);
  invf.NTT(true, maxl), df.NTT(true, maxl);
  for (int i = 0; i < maxl; ++i) Res[i] = static_cast<ll>(df[i]) * invf[i] % P;
  Res.NTT(false, maxl);
  for (int i = n - 1; i; --i) Res[i] = static_cast<ll>(Res[i - 1]) * qpow(i, P - 2) % P;
  Res[0] = 0;
}

inline void Polynomial::exp(Polynomial &Res, int const n) const {
  static Polynomial lnres;
  if (n == 1) { return Res[0] = 1, void(); }
  exp(Res, (n + 1) >> 1);
  Res.ln(lnres, n);
  int maxl = 1;
  while (maxl < n << 1) maxl <<= 1;
  for (int i = 0; i < n; ++i)
    if ((lnres[i] = F[i] - lnres[i]) < 0) lnres[i] += P;
  ++lnres[0];
  if (lnres[0] >= P) lnres[0] -= P;
  lnres.NTT(true, maxl), Res.NTT(true, maxl);
  for (int i = 0; i < maxl; ++i) Res[i] = static_cast<ll>(Res[i]) * lnres[i] % P;
  Res.NTT(false, maxl);
  for (int i = n; i < maxl; ++i) Res[i] = 0;
}

inline void Polynomial::pow(Polynomial &Res, int const n, int const k) const {
  static Polynomial tmp;
  ln(tmp, n);
  for (int i = 0; i < n; ++i) tmp[i] = 1ll * tmp[i] * k % P;
  tmp.exp(Res, n);
}

/// @param Q Quotient @param R Remainder @param D divider
inline void Polynomial::div(Polynomial &Q, Polynomial &R, Polynomial const &D,
                            int const n, int const m) const {
  static Polynomial Dr, iDr, Fr;
  for (int i = 0; i <= m; ++i) Dr[i] = D[m - i];
  for (int i = 0; i <= n; ++i) Fr[i] = F[n - i];
  Dr.inv(iDr, n - m + 1);
  int maxl = 1;
  while (maxl <= n + m) maxl <<= 1;
  iDr.NTT(true, maxl << 1), Fr.NTT(true, maxl << 1);
  for (int i = 0; i < maxl << 1; ++i) Q[i] = 1ll * iDr[i] * Fr[i] % P;
  Q.NTT(false, maxl << 1);
  for (int i = n - m + 1; i < maxl << 1; ++i) Q[i] = 0;
  std::reverse(&Q, &Q + n - m + 1);
  std::reverse(&Dr, &Dr + m + 1);
  Q.NTT(true, maxl << 1), Dr.NTT(true, maxl << 1);
  for (int i = 0; i < maxl << 1; ++i) R[i] = 1ll * Q[i] * Dr[i] % P;
  R.NTT(false, maxl << 1), Q.NTT(false, maxl << 1);
  for (int i = 0; i < m; ++i)
    if ((R[i] = F[i] - R[i]) < 0) R[i] += P;
}

/// @param Res 应清空 @param n 模的次数（项数）
inline void Polynomial::sqrt(Polynomial &Res, int const n) const {
  static Polynomial Resinv, copyF;
  if (n == 1) { return Res[0] = 1, void(); }
  sqrt(Res, (n + 1) >> 1);
  for (int i = 0; i <= n << 1; ++i) Resinv[i] = 0;
  Res.inv(Resinv, n);
  int maxl = 1;
  while (maxl < n << 1) maxl <<= 1;
  memcpy(&copyF, F, sizeof(int) * n);
  for (int i = n; i < maxl; ++i) copyF[i] = 0;
  copyF.NTT(true, maxl), Resinv.NTT(true, maxl), Res.NTT(true, maxl);
  for (int i = 0; i < maxl; ++i)
    /// @warning        P 过大时可能会炸 ↓
    if ((Res[i] = 1ll * inv2 * (Res[i] + 1ll * copyF[i] * Resinv[i] % P) % P) >= P)
      Res[i] -= P;
  Res.NTT(false, maxl);
  for (int i = n; i < maxl; ++i) Res[i] = 0;
}

int main() {
#ifndef ONLINE_JUDGE
#ifdef LOCAL
  freopen("/tmp/CodeTmp/testdata.in", "r", stdin);
  freopen("/tmp/CodeTmp/testdata.out", "w", stdout);
#else
  freopen("Polynomial.in", "r", stdin);
  freopen("Polynomial.out", "w", stdout);
#endif
#endif

  read(n);
  for (int i = 0; i < n; ++i) read(F[i]);
  F.pre(n << 2);
  F.sqrt(G, n);
  for (int i = 0; i < n; ++i) printf("%d ", G[i]);
  return 0;
}