多项式一环魔法咏唱
FFT
经典板子，采用三步并两步优化。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long double ld;
typedef pair<ll, ll> PII;
typedef pair<ll, PII> PPI;
const int N = 1.5e6 + 10;
const double Pi = acos(-1);
int n, m, rev[N << 1];
struct Complex {
double x, y;
Complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
Complex operator+ (Complex const &C) const {
return Complex(x + C.x, y + C.y);
}
Complex operator- (Complex const &C) const {
return Complex(x - C.x, y - C.y);
}
Complex operator* (Complex const &C) const {
return Complex(x * C.x - y * C.y, x * C.y + y * C.x);
}
}f[N << 1];
void fft(Complex *f, int op) {
for (int i = 0; i < n; i ++ ) {
if (i < rev[i]) {
swap(f[i], f[rev[i]]);
}
}
for (int k = 1; k < n; k <<= 1) {
Complex w1 = Complex(cos(Pi / k), op * sin(Pi / k));
for (int i = 0; i < n; i += k << 1) {
Complex wk = Complex(1, 0);
for (int j = 0; j < k; j ++ , wk = wk * w1) {
Complex x = f[i + j], y = wk * f[i + j + k];
f[i + j] = x + y, f[i + j + k] = x - y;
}
}
}
if (op == -1) for (int i = 0; i < n; i ++ ) f[i].x /= n, f[i].y /= n;
}
void solve() {
cin >> n >> m;
for (int i = 0; i <= n; i ++ ) cin >> f[i].x;
for (int i = 0; i <= m; i ++ ) cin >> f[i].y;
for (m += n, n = 1; n <= m; n <<= 1);
for (int i = 0; i < n; i ++ ) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
fft(f, 1);
for (int i = 0; i < n; i ++ ) f[i] = f[i] * f[i];
fft(f, -1);
for (int i = 0; i <= m; i ++ ) cout << (int)(f[i].y / 2 + 0.5) << " ";
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
int T = 1;
// cin >> T;
while (T -- ) solve();
return 0;
}
NTT
封装一下我们的 NTT。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define clr(f, n) memset(f, 0, sizeof(int) * (n))
#define cpy(f, g, n) memcpy(f, g, sizeof(int) * (n))
const int N = 1050000, _G = 3, mod = 998244353;
ll qpow(ll a, ll k = mod - 2) {
ll res = 1;
while (k) {
if (k & 1) res = res * a % mod;
k >>= 1;
a = a * a % mod;
}
return res;
}
const int invG = qpow(_G);
int rev[N << 1], rev_len;
void rev_init(int n) {
if (n == rev_len) return;
for (int i = 0; i < n; i ++ ) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
rev_len = n;
}
void NTT(int *g, int op, int n) {
rev_init(n);
static ull f[N << 1], Gk[N << 1] = {1};
for (int i = 0; i < n; i ++ ) f[i] = g[rev[i]];
for (int k = 1; k < n; k <<= 1) {
int G1 = qpow(~op ? _G : invG, (mod - 1) / (k << 1));
for (int i = 1; i < k; i ++ ) Gk[i] = Gk[i - 1] * G1 % mod;
for (int i = 0; i < n; i += k << 1) {
for (int j = 0; j < k; j ++ ) {
int tmp = Gk[j] * f[i | j | k] % mod;
f[i | j | k] = f[i | j] + mod - tmp;
f[i | j] += tmp;
}
}
if (k == (1 << 10)) for (int i = 0; i < n; i ++ ) f[i] %= mod;
}
if (~op) for (int i = 0; i < n; i ++ ) g[i] = f[i] % mod;
else {
int invn = qpow(n);
for (int i = 0; i < n; i ++ ) g[i] = f[i] % mod * invn % mod;
}
}
void px(int *f, int *g, int n) {
for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i] * g[i] % mod;
}
void covolution(int *f, int *g, int len, int lim) {
static int sav[N << 1];
int n; for (n = 1; n < len << 1; n <<= 1);
clr(sav, n); cpy(sav, g, n);
NTT(sav, 1, n); NTT(f, 1, n);
px(f, sav, n); NTT(f, -1, n);
clr(f + lim, n - lim), clr(sav, n);
}
int F[N << 1], G[N << 1];
void solve() {
int n, m;
cin >> n >> m;
for (int i = 0; i <= n; i ++ ) cin >> F[i];
for (int i = 0; i <= m; i ++ ) cin >> G[i];
covolution(F, G, n + m + 1 >> 1, n + m + 1);
for (int i = 0; i <= n + m; i ++ ) cout << F[i] << " ";
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
int T = 1;
// cin >> T;
while (T -- ) solve();
return 0;
}
多项式求逆
考虑求 \(F(x)G(x) \equiv 1 \pmod{x^n}\) 中的 \(G(x)\)。
若考虑分治，设我们已经知道了 \(G(x)\) 的低 \(x^{\frac{n}{2}}\) 项系数多项式 \(G_*(x)\)，看看能不能求出 \(G(x)\)
\[G(x) = F^{-1}(x) \pmod{x^n}
\]
\[G_*(x) = F^{-1}(x) \pmod{x^{\frac{n}{2}}}
\]
由于 \(G(x) = G_*(x) \pmod{x^{\frac{n}{2}}}\)，考虑作差有
\[G(x) - G_*(x) = 0 \pmod{x^{\frac{n}{2}}}
\]
为了实现扩域，两边同时平方有
\[\left(G(x) - G_*(x)\right)^2 = 0 \pmod{x^n}
\]
也就是说
\[G(x)^2 - 2G(x)G_*(x) + G_*(x)^2 = 0 \pmod{x^n}
\]
两边乘 \(F(x)\) 有
\[G(x) = 2G_*(x) - G_*(x)^2F(x) \pmod{x^n}
\]
由于 \(G_1(x) = F_1^{-1}(x)\)，那么之后的项可以通过递推得到。
我们发现 \(G_*(x)^2F(x)\) 是不超过 \(2n\) 次的多项式，但是我们计算一定会用到 \(2n\) 长度的卷积，那有没有什么办法卡掉这个常数呢？
当然是有的，考虑 \(G(x)\) 与 \(F(x)\) 乘积的意义，\(G_*(x)F_n(x) = 1 + Q(x) \pmod{x^n}\)，其中 \(Q(x)\) 为 \(x^{\frac{n}{2}} \sim x^{n - 1}\) 的项的组合，那么 \(G_*(x)^2F(x) = G_*(x) + Q(x)G_*(x) \pmod{x^n}\)，其中 \(G_*(x)\) 对 \(\frac{n}{2}\) 以后的无贡献，也就是说我们实际上只需要 \(Q(x)\) 与 \(G_*(x)\) 卷积后的后 \(\frac{n}{2}\) 项即可，前面的项对我们是无意义的，所以我们利用这个性质可以通过下面的操作来实现卡常。
由于 NTT 是定义在 \(n = 2^k\) 序列上的循环卷积，我们用 \(F(x)\) 与 \(G_*(x)\) 线性卷积后，得到的会是 \(0 \sim x^{\frac{3n}{2} - 1}\) 的卷积序列，由于我们进行 \(n\) 长度的循环卷积，那么属于 \(x^{n} \sim x^{\frac{3n}{2} - 1}\) 部分的贡献会计算到 \(x^{0} \sim x^{\frac{n}{2} - 1}\) 上，由于 \(1\) 和卷积溢出部分的贡献我们都不需要，所以可以直接把前置位清空进行下一次卷积，这样卷积出来的结果取后 \(\frac{n}{2}\) 就是正确的。
void Poly_inv(int *f, int m) {
int n; for (n = 1; n < m; n <<= 1);
static int g1[N << 1], g2[N << 1], sav[N << 1];
g1[0] = qpow(f[0]);
for (int len = 2; len <= n; len <<= 1) {
cpy(g2, g1, len >> 1), cpy(sav, f, len);
NTT(g2, 1, len), NTT(sav, 1, len);
px(g2, sav, len); NTT(g2, -1, len);
clr(g2, len >> 1); cpy(sav, g1, len);
NTT(sav, 1, len); NTT(g2, 1, len);
px(g2, sav, len); NTT(g2, -1, len);
for (int i = len >> 1; i < len; i ++ ) g1[i] = (2ll * g1[i] - g2[i] + mod) % mod;
}
cpy(f, g1, m), clr(g1, n), clr(g2, n), clr(sav, n);
}
多项式牛顿迭代
已知 \(G(F(x)) = 0\)，\(G(x)\) 的系数已知，求 \(F(x)\)。
考虑倍增，若我们已知 \(F_*(x)\) 使得 \(G(F_*(x)) = 0 \pmod{x^{\frac{n}{2}}}\)，想要求得 \(G(F(x)) \equiv 0 \pmod{x^n}\) 的下一组解，我们将 \(G(F(x))\) 在 \(F_*(x)\) 处泰勒展开：
\[G(F(x)) = G(F_*(x)) + \frac{G'(F_*(x))}{1!}\left(F(x) - F_*(x)\right) + \frac{G''(F_*(x))}{2!}\left(F(x) - F_*(x)\right)^2 + \cdots
\]
由于 \(F_n(x)\) 与 \(F_{\frac{n}{2}}(x)\) 前 \(\frac{n}{2}\) 项相同，因此 \(\left(F_n(x) - F_{\frac{n}{2}}(x)\right)^2\) 至少是 \(n\) 次项的，在模 \(x^n\) 意义下等于 \(0\)，于是
\[G(F(x)) = G(F_*(x)) + \frac{G'\left(F_*(x)\right)}{1!}\left(F(x) - F*(x)\right) \pmod{x^n}
\]
那么 \(F(x) = F_*(x) - \frac{G\left(F_*(x)\right)}{G'\left(F_*(x)\right)} \pmod{x^n}\)
其中 \(G'\left(F_*(x)\right)\) 的精度达到 \(\frac{n}{2}\) 就可以了，因为 \(G\left(F_*(x)\right)\) 的前 \(\frac{n}{2}\) 都为 \(0\)。
这里的 \(G(x)\) 通常都是好算的，并且我们通常把 \(G\) 当做常量。
-
多项式求逆
考虑 \(A(x)B(x) \equiv 1 \pmod{x^n}\)，\(A(x)\) 已知，考虑转换其形式为 \(A(x)B(x) - 1 \equiv 0 \pmod{x^n}\)。设多项式 \(G(B(x)) = A(x)B(x) - 1\) 在 \(B(x)\) 处的取值为 \(0\)，那么 \(G'(B(x)) = A(x)\)，这里我们把 \(A(x)\) 看作常量，\(B(x)\) 看作主元，利用牛顿迭代公式，我们可以得到
\[B(x) \equiv B_*(x) - \frac{G\left(B_*(x)\right)}{G'\left(B_*(x)\right)} \equiv B_*(x) - \frac{A(x)B_*(x) - 1}{A(x)} \pmod{x^n} \]实际上就是 \(B(x) \equiv B_*(x) - B(x)\left(A(x)B_*(x) - 1\right) \pmod{x^n}\)，需要注意到 \(A(x)B_*(x)\) 前 \(\frac{n}{2}\) 项为 \(1\)，这决定了 \(\left(A(x)B_*(x) - 1\right)\) 最低项是 \(x^{\frac{n}{2}}\) 所以 \(B(x)\) 的精度为 \(x^{\frac{n}{2}}\)，因此
\[B(x) \equiv 2B_*(x) - B_*(x)^2A(x) \pmod{x^n} \]这和我们直接推导的公式是一致的。
-
多项式开根
考虑 \(A(x) \equiv B(x)^2\)，考虑转换其形式为 \(B(x)^2 - A(x) = 0\)。设 \(G(B(x)) = B(x)^2 - A(x)\)，那么 \(G'(B(x)) = 2B(x)\)，利用公式可以得到
\[B(x) = B_*(x) - \frac{G(B_*(x))}{G'(B_*(x))} = B_*(x) - \frac{B_*(x)^2 - A(x)}{2B_*(x)} = \frac{A(x) + B_*^2(x)}{2B_*(x)} \]所有精度都需为 \(x^n\)。
void Poly_sqrt(int *f, int m) {
int n; for (n = 1; n < m; n <<= 1);
static int b1[N << 1], b2[N << 1];
b1[0] = 1;
for (int len = 2; len <= n; len <<= 1) {
for (int i = 0; i < len >> 1; i ++ ) b2[i] = (b1[i] << 1) % mod;
Poly_inv(b2, len);
NTT(b1, 1, len); px(b1, b1, len); NTT(b1, -1, len);
for (int i = 0; i < len; i ++ ) b1[i] = (f[i] + b1[i]) % mod;
covolution(b1, b2, len, len);
}
cpy(f, b1, m); clr(b1, n << 1); clr(b2, n << 1);
}
多项式带余除法
给定 \(F(x), \, G(x)\)，其中 \(F(x)\) 为 \(n\) 次多项式，\(G(x)\) 为 \(m\) 次多项式，请求出满足 \(F(x) = G(x)Q(x) + R(x)\) 的 \(Q(x)\) 与 \(R(x)\)。
考虑到 \(Q(x), \, R(x)\) 的最高次较小，我们可以选择翻转系数的方式来消去该部分。
因为
\[F\left(\frac{1}{x}\right) = G\left(\frac{1}{x}\right)Q\left(\frac{1}{x}\right) + Q\left(\frac{1}{x}\right)
\]
即
\[x^nF\left(\frac{1}{x}\right) = x^{m}G\left(\frac{1}{x}\right)x^{n - m}Q\left(\frac{1}{x}\right) + x^nQ\left(\frac{1}{x}\right)
\]
记 \(F^T(x)\) 表示 \(F(x)\) 将高次项与低次项系数交换的翻转，也就是
\[F^T(x) = G^T(x)Q^T(x) + x^{n - m + 1}R^T(x)
\]
在模 \(x^{n - m + 1}\) 意义下，可以消去余数，即
\[F^T(x) \equiv G^T(x)Q^T(x) \pmod{x^{n - m + 1}}
\]
那么考虑求多项式逆元就可以求出
\[Q^T(x) \equiv \frac{F^T(x)}{G^T(x)} \pmod{x^{n - m + 1}}
\]
求出 \(Q(x)\) 后，通过简单的移项操作就可以求出
\[R(x) = F(x) - G(x)Q(x)
\]
void Poly_div(int *f, int *g, int n, int m) {
static int Q[N << 1], R[N << 1];
int len = n - m + 1;
rev(f, n); cpy(R, f, len); rev(f, n);
rev(g, m); cpy(Q, g, len); rev(g, m);
Poly_inv(Q, len); covolution(Q, R, len, len); rev(Q, len);
covolution(g, Q, n, n);
for (int i = 0; i < m - 1; i ++ ) g[i] = (f[i] - g[i] + mod) % mod;
clr(g + m - 1, len); cpy(f, Q, len); clr(f + len, n - len);
}
多项式积分/求导
这个很容易，设 \(F[i]\) 表示 \(i\) 次项系数，也就是
\[f'(x) = \sum_{i = 0} iF[i]x^{i - 1}
\]
\[\int f(x)dx = \sum_{i = 0} \frac{F[i]}{i + 1}x^{i + 1} + C
\]
int inv[N << 1], inv_len;
void inv_init(int n) {
if (n <= inv_len) return;
if (!inv_len) inv[0] = inv[1] = 1, inv_len = 1;
for (int i = inv_len + 1; i <= n; i ++ ) inv[i] = 1ll * inv[mod % i] * (mod - mod / i) % mod;
inv_len = n;
}
void Poly_d(int *f, int n) {
for (int i = 1; i < n; i ++ ) f[i - 1] = 1ll * f[i] * i % mod;
f[n - 1] = 0;
}
void Poly_int(int *f, int n) {
for (int i = n; i; i -- ) f[i] = 1ll * f[i - 1] * inv[i] % mod;
f[0] = 0;
}
多项式 \(\ln\)
考虑 \(\ln A(x) = B(x)\)，两边同时求导，根据链式法则可得
\[A'(x)\frac{1}{A(x)} = B'(x)
\]
两边积分可得
\[B(x) = \int \frac{A'(x)}{A(x)}dx
\]
代码比较好写
void Poly_ln(int *f, int n) {
static int sav[N << 1];
inv_init(n); cpy(sav, f, n);
Poly_d(sav, n); Poly_inv(f, n);
covolution(f, sav, n, n); Poly_int(f, n - 1);
clr(sav, n);
}
多项式 \(\exp\)
由于 \(e^{A(x)} = B(x) \Leftrightarrow A(x) = \ln B(x) \Leftrightarrow \ln B(x) - A(x) = 0\)
根据牛顿迭代法 \(F(B(x)) = \ln B(x) - A(x)\)，\(F'(B(x)) = \frac{1}{B(x)}\)，则有
\[B(x) = B_*(x) - \frac{F(B_*(x))}{F'(B_*(x))} = B_*(x)\left(1 - \ln B_*(x) + A(x)\right)
\]
必须保证 \(A[0] = 0\)，则 \(B[0] = 1\)，否则常数项为无穷级数。
注意，\(\ln B_*(x)\) 的精度需要扩展到 \(x^n\)。
void Poly_exp(int *f, int m) {
static int b1[N << 1], b2[N << 1];
int n; for (n = 1; n < m; n <<= 1);
b1[0] = 1;
for (int len = 2; len <= n; len <<= 1) {
cpy(b2, b1, len >> 1); Poly_ln(b2, len);
for (int i = 0; i < len; i ++ ) b2[i] = (f[i] - b2[i] + mod) % mod;
b2[0] = (b2[0] + 1) % mod;
covolution(b1, b2, len, len);
}
cpy(f, b1, m); clr(b1, n); clr(b2, n);
}
多项式快速幂
容易想到 \(A^k(x) = e^{k\ln A(x)}\)，如果保证了 \(A[0] = 1\)，那么可以快速写出代码。
void Poly_pow(int *f, int n, int k) {
Poly_ln(f, n);
for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i] * k % mod;
Poly_exp(f, n);
}
如果不保证 \(A[0] = 1\) 呢？我们也可以通过乘一个单项式 \(cx^p\)，使得
\[A^k(x) = \left(cx^pB(x)\right)^k = c^kx^{pk} \cdot B^k(x)
\]
其中 \(B[0] = 1\)，然后利用降幂定理快速计算即可。
void Poly_pow(int *f, int n, string k) {
int k1 = 0, k2 = 0, p = 0, c;
while (!f[p]) p ++ ;
for (int i = 0; k[i]; i ++ ) {
k1 = (10ll * k1 + k[i] - '0') % mod;
k2 = (10ll * k2 + k[i] - '0') % (mod - 1);
if (1ll * k1 * p >= n) return clr(f, n), void();
}
n -= p * k1; c = qpow(f[p]);
for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i + p] * c % mod;
clr(f + n, p * k1); Poly_ln(f, n);
for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i] * k1 % mod;
Poly_exp(f, n); c = qpow(c, mod - 1 - k2);
for (int i = n - 1; i >= 0; i -- ) f[p * k1 + i] = 1ll * f[i] * c % mod;
clr(f, p * k1);
}
然后我们进行一个超级拼装
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define clr(f, n) memset(f, 0, sizeof(int) * (n))
#define cpy(f, g, n) memcpy(f, g, sizeof(int) * (n))
#define rev(f, n) reverse(f, f + (n))
const int N = 150000, _G = 3, mod = 998244353;
ll qpow(ll a, ll k = mod - 2) {
ll res = 1;
while (k) {
if (k & 1) res = res * a % mod;
k >>= 1;
a = a * a % mod;
}
return res;
}
const int invG = qpow(_G);
int rev[N << 1], rev_len;
void rev_init(int n) {
if (rev_len == n) return;
for (int i = 0; i < n; i ++ ) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
rev_len = n;
}
void NTT(int *g, int op, int n) {
rev_init(n);
static ull f[N << 1], Gk[N << 1] = {1};
for (int i = 0; i < n; i ++ ) f[i] = g[rev[i]];
for (int k = 1; k < n; k <<= 1) {
int G1 = qpow(~op ? _G : invG, (mod - 1) / (k << 1));
for (int i = 1; i < k; i ++ ) Gk[i] = Gk[i - 1] * G1 % mod;
for (int i = 0; i < n; i += k << 1) {
for (int j = 0; j < k; j ++ ) {
int tmp = Gk[j] * f[i | j | k] % mod;
f[i | j | k] = f[i | j] + mod - tmp;
f[i | j] += tmp;
}
}
if (k == (1 << 10)) for (int i = 0; i < n; i ++ ) f[i] %= mod;
}
if (~op) for (int i = 0; i < n; i ++ ) g[i] = f[i] % mod;
else {
int invn = qpow(n);
for (int i = 0; i < n; i ++ ) g[i] = f[i] % mod * invn % mod;
}
}
void px(int *f, int *g, int n) {
for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i] * g[i] % mod;
}
int inv[N << 1], inv_len;
void inv_init(int n) {
if (n <= inv_len) return;
if (!inv_len) inv[0] = inv[1] = 1, inv_len = 1;
for (int i = inv_len + 1; i <= n; i ++ ) inv[i] = 1ll * inv[mod % i] * (mod - mod / i) % mod;
inv_len = n;
}
void Poly_d(int *f, int n) {
for (int i = 1; i < n; i ++ ) f[i - 1] = 1ll * f[i] * i % mod;
f[n - 1] = 0;
}
void Poly_int(int *f, int n) {
for (int i = n; i; i -- ) f[i] = 1ll * f[i - 1] * inv[i] % mod;
f[0] = 0;
}
void covolution(int *f, int *g, int len, int lim) {
static int sav[N << 1];
int n; for (n = 1; n < len << 1; n <<= 1);
clr(sav, n); cpy(sav, g, n);
NTT(sav, 1, n); NTT(f, 1, n);
px(f, sav, n); NTT(f, -1, n);
clr(f + lim, n - lim), clr(sav, n);
}
void Poly_inv(int *f, int m) {
static int g1[N << 1], g2[N << 1], sav[N << 1];
int n; for (n = 1; n < m; n <<= 1);
g1[0] = qpow(f[0]);
for (int len = 2; len <= n; len <<= 1) {
cpy(g2, g1, len >> 1), cpy(sav, f, len);
NTT(g2, 1, len), NTT(sav, 1, len);
px(g2, sav, len); NTT(g2, -1, len);
clr(g2, len >> 1); cpy(sav, g1, len);
NTT(sav, 1, len); NTT(g2, 1, len);
px(g2, sav, len); NTT(g2, -1, len);
for (int i = len >> 1; i < len; i ++ ) g1[i] = (2ll * g1[i] - g2[i] + mod) % mod;
}
cpy(f, g1, m), clr(g1, n), clr(g2, n), clr(sav, n);
}
void Poly_sqrt(int *f, int m) {
static int b1[N << 1], b2[N << 1];
int n; for (n = 1; n < m; n <<= 1);
b1[0] = 1;
for (int len = 2; len <= n; len <<= 1) {
for (int i = 0; i < len >> 1; i ++ ) b2[i] = (b1[i] << 1) % mod;
Poly_inv(b2, len);
NTT(b1, 1, len); px(b1, b1, len); NTT(b1, -1, len);
for (int i = 0; i < len; i ++ ) b1[i] = (f[i] + b1[i]) % mod;
covolution(b1, b2, len, len);
}
cpy(f, b1, m); clr(b1, n << 1); clr(b2, n << 1);
}
void Poly_div(int *f, int *g, int n, int m) {
static int Q[N << 1], R[N << 1];
int len = n - m + 1;
rev(f, n); cpy(R, f, len); rev(f, n);
rev(g, m); cpy(Q, g, len); rev(g, m);
Poly_inv(Q, len); covolution(Q, R, len, len); rev(Q, len);
covolution(g, Q, n, n);
for (int i = 0; i < m - 1; i ++ ) g[i] = (f[i] - g[i] + mod) % mod;
clr(g + m - 1, len); cpy(f, Q, len); clr(f + len, n - len);
}
void Poly_ln(int *f, int n) {
static int sav[N << 1];
inv_init(n); cpy(sav, f, n);
Poly_d(sav, n); Poly_inv(f, n);
covolution(f, sav, n, n); Poly_int(f, n - 1);
clr(sav, n);
}
void Poly_exp(int *f, int m) {
static int b1[N << 1], b2[N << 1];
int n; for (n = 1; n < m; n <<= 1);
b1[0] = 1;
for (int len = 2; len <= n; len <<= 1) {
cpy(b2, b1, len >> 1); Poly_ln(b2, len);
for (int i = 0; i < len; i ++ ) b2[i] = (f[i] - b2[i] + mod) % mod;
b2[0] = (b2[0] + 1) % mod;
covolution(b1, b2, len, len);
}
cpy(f, b1, m); clr(b1, n); clr(b2, n);
}
void Poly_pow(int *f, int n, string k) {
int k1 = 0, k2 = 0, p = 0, c;
while (!f[p]) p ++ ;
for (int i = 0; k[i]; i ++ ) {
k1 = (10ll * k1 + k[i] - '0') % mod;
k2 = (10ll * k2 + k[i] - '0') % (mod - 1);
if (1ll * k1 * p >= n) return clr(f, n), void();
}
n -= p * k1; c = qpow(f[p]);
for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i + p] * c % mod;
clr(f + n, p * k1); Poly_ln(f, n);
for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i] * k1 % mod;
Poly_exp(f, n); c = qpow(c, mod - 1 - k2);
for (int i = n - 1; i >= 0; i -- ) f[p * k1 + i] = 1ll * f[i] * c % mod;
clr(f, p * k1);
}
浙公网安备 33010602011771号