emm……NTT左转去看各位神犇的博客吧QwQ
这里只贴代码及部分操作的推导过程。
首先是喜闻乐见的 NTT 多项式乘法板子（这个就不解释了）
namespace NTT{
ll lim, len;
inline ll qpow(ll a, ll b){
ll res = 1;
while(b){
if(b & 1) res = res * a % mod;
a = a * a % mod, b >>= 1;
}
return res;
}
inline void get_rev(cl n){
lim = 1, len = 0;
while(lim < n) lim <<= 1, ++len;
for(int i = 0; i < lim; ++i) p[i] = (p[i >> 1] >> 1) | ((i & 1) << (len - 1));
}
inline void ntt(ll A[], cl lim, cl type){
for(int i = 0; i < lim; ++i)
if(i < p[i]) swap(A[i], A[p[i]]);
for(int mid = 1; mid < lim; mid <<= 1){
ll Wn = qpow(type == 1 ? G : Gi, (mod - 1) / (mid << 1));
for(int i = 0; i < lim; i += (mid << 1)){
ll w = 1;
for(int j = 0; j < mid; ++j, w = w * Wn % mod){
ll x = A[i + j], y = w * A[i + j + mid] % mod;
A[i + j] = (x + y) % mod;
A[i + j + mid] = (x - y + mod) % mod;
}
}
}
if(type == 1) return;
ll inv = qpow(lim, mod - 2);
for(int i = 0; i < lim; ++i) A[i] = A[i] * inv % mod;
}
inline void Mul(cl n, cl m, ll a[], ll b[]){
get_rev(n + m);
ntt(a, lim, 1), ntt(b, lim, 1);
for(int i = 0; i < lim; ++i) a[i] = a[i] * b[i] % mod;
ntt(a, lim, -1);
}
}
using namespace NTT;
紧跟着的是求逆
已知多项式 \(F(x)\)，要求 \(G(x)F(x) \equiv 1\ \ (mod\ \ x ^ n)\)
设 \(H(x)F(x) \equiv 1 \ \ (mod\ x^{⌈ \frac{n}{2} ⌉})\)
那么 $$(G(x) - H(x)) \equiv 0\ \ (mod\ \ x ^ {⌈ \frac{n}{2} ⌉})$$
两边同时平方：
\[(G(x) - H(x))^2 \equiv 0\ \ (mod\ \ x ^ {⌈ \frac{n}{2} ⌉})
\]
\[G(x)^2 - 2G(x)H(x) + H(x)^2 \equiv 0\ \ (mod\ \ x ^ n)
\]
两边再同时乘上 \(F(x)\)，由于 \(G(x)F(x) \equiv 1\ \ (mod\ \ x ^ n)\)：
\[G(x) - 2H(x) + H(x)^2F(x) \equiv 0\ \ (mod\ \ x ^ n)
\]
移一下项并提取公因多项式（：
\[G(x) = H(x)(2 - H(x)F(x))
\]
inline void Inv(ll n, ll a[], ll b[]){
if(n == 1) return b[0] = qpow(a[0], mod - 2), void();
Inv((n + 1) >> 1, a, b);
get_rev(n << 1);
for(int i = 0; i < n; ++i) c[i] = a[i];
for(int i = n; i < lim; ++i) c[i] = 0;
ntt(c, lim, 1), ntt(b, lim, 1);
for(int i = 0; i < lim; ++i)
b[i] = (2ll - c[i] * b[i] % mod + mod) * b[i] % mod;
ntt(b, lim, -1);
for(int i = n; i < lim; ++i) b[i] = 0;
}
然后是多项式开根
我们已知 \(F(x)\)，要求 \(G(x)\) 使得 \(G(x)^2 \equiv F(x)\ \ (mod\ \ x^n)\)
设 \(H(x)^2 \equiv F(x)\ \ (mod\ x^{⌈ \frac{n}{2} ⌉})\)
那么
\[G(x) \equiv H(x)\ \ (mod\ x^{⌈ \frac{n}{2} ⌉})
\]
\[(G(x) - H(x)) \equiv 0\ \ (mod\ x^{⌈ \frac{n}{2} ⌉})
\]
还是两边同时平方：
\[(G(x) - H(x))^2 \equiv 0\ \ (mod\ x^n)
\]
\[G(x)^2 - 2G(x)H(x) + H(x)^2 \equiv 0\ \ (mod\ \ x ^ n)
\]
由于 \(G(x)^2 \equiv F(x)\ \ (mod\ \ x^n)\) 得：
\[F(x) - 2G(x)H(x) + H(x)^2 \equiv 0 \ \ (mod\ \ x ^ n)
\]
移项：
\[G(x) = \frac{F(x) + H(x)^2}{2H(x)}
\]
那么事实上我的代码里的式子是这样的：
\[G(x) = \frac{F(x)H(x)^{-1} + H(x)}{2}
\]
\(b(x) \rightarrow H(x)\)
\(d(x) \rightarrow F(x)\)
\(e(x) \rightarrow H(x)^{-1}\)
inline void Sqrt(ll n, ll a[], ll b[]){
if(n == 1) return b[0] = 1, void();
Sqrt((n + 1) >> 1, a, b);
get_rev(n << 1);
memset(e, 0, sizeof(e));
Inv(n, b, e);
for(int i = 0; i < n; ++i) d[i] = a[i];
for(int i = n; i < lim; ++i) d[i] = 0;
ntt(d, lim, 1), ntt(e, lim, 1), ntt(b, lim, 1);
for(int i = 0; i < lim; ++i) b[i] = (b[i] + d[i] * e[i] % mod) * inv2 % mod;
ntt(b, lim, -1);
for(int i = n; i < lim; ++i) b[i] = 0;
}
接着自然是是求 ln
这个就需要通道一点微积分的知识了，虽然我也不太会QwQ，但是我会贺代码(｡･ω･｡)
还是先来推推式子吧。
已知 \(F(x)\)，求 \(G(x) = lnF(x)\)
众所周知，\(g(x) = ln\,f(x)\) 的导数为 \(g'(x) = \frac{f'(x)}{f(x)}\)，所以：
\[G'(x) = \frac{F'(x)}{F(x)}
\]
也就是说，\(G'(x)\) 我们可以直接求出来了，下面给出不定积分的一个运算公式：
\[\int x^a\,dx = \frac{1}{a + 1} x ^ {a + 1}
\]
由于我们的 \(G(x) = a_0 + a_1x + a_2x^2 + ···+ a_nx^n\)
所以对每一项求个积分就完了，具体看代码吧。
namespace Ln{
inline void Diff(cl n, ll a[], ll b[]){//微分求导
for(int i = 1; i < n; ++i) b[i - 1] = i * a[i];
b[n - 1] = 0;
}
inline void Integral(cl n, ll a[], ll b[]){//积分
for(int i = 1; i < n; ++i) b[i] = a[i - 1] * qpow(i, mod - 2) % mod;
b[0] = 0;
}
inline void ln(cl n, ll a[], ll b[]){
memset(f1, 0, sizeof(f1));
memset(f2, 0, sizeof(f2));
Diff(n, a, f1), Inv(n, a, f2);//f1(x) = a'(x)，f2(x) = a(x)^(-1)
Mul(n, n, f1, f2);//f1(x) = f1(x) * f2(x)
Integral(n, f1, b);//g(x) = f(x) 的积分
}
}
using namespace Ln;
exp 当然也不能少
已知 \(F(x)\)，求 \(G(x) = e^{F(x)}\)
emm……要用到泰勒展开，牛顿迭代什么的，然鹅我太蒻了，还不会QwQ。
所以这里就只有结论了，像我这种蒟蒻还是全文背诵吧。
\[F(x) = F_0(x) - \frac{G(F_0(x))}{G'(F_0(x)}
\]
再经过一番清新简单的推导过程之后……
\[G(x) \equiv G_0(x)(1 - lnG_0(x) + F(x))\ \ (mod\ \ x ^ n)
\]
所以就可以算了……
inline void Exp(cl n, ll a[], ll b[]){
if(n == 1) return b[0] = 1, void();
Exp((n + 1) >> 1, a, b);
get_rev(n << 1);
ln(n, b, d);// d(x) = ln(b(x))
for(int i = 0; i < n; ++i) e[i] = a[i];
for(int i = n; i < lim; ++i) e[i] = 0;
ntt(d, lim, 1), ntt(e, lim, 1), ntt(b, lim, 1);
for(int i = 0; i < lim; ++i) b[i] = (1ll - d[i] + e[i] + mod) * b[i] % mod;
ntt(b, lim, -1);
for(int i = n; i < lim; ++i) b[i] = 0;
}
最后是多项式除法
除法由于其优秀的边界故不和上面放到一起，并且这里给出完整代码。
以及我这里是计算出 \(Q(x)\) 和 \(R(x)\) 后一块输出，由于计算 \(R(x)\) 时需要多项式卷积一下，因此对 \(Q(x)\) 做了一遍 \(NTT\)，所以卷积完之后一定记得要 \(NTT\) 回来啊。
来看看推导过程。
已知 \(F(x)\)，\(G(x)\)，求 \(F(x) = Q(x)G(x) + R(x)\)
首先就要用到一个神奇的操作，设 \(F_R(x)\) 是 \(F(x)\) 把系数反过来之后的多项式，即
\(F_R(x) = a_n + a_{n - 1}x + a_{n - 2}x^2 +···+a_0x^n\)
容易发现，\(F_R(x) = x^nF(\frac{1}{x})\)
然后就可以愉快的推式子啦 :）
\[F(x) = Q(x)G(x) + R(x)
\\
F(\frac{1}{x}) = Q(\frac{1}{x})G(\frac{1}{x}) + R(\frac{1}{x})
\]
\(F(x)\) 是 \(n\) 次的，\(G(x)\) 是 \(m\) 次的，\(Q(x)\) 是 \(n - m\) 次的，\(R(x)\) 是 \(n - m - 1\) 次的，所以两边同乘 \(x^n\)：
\[x^nF(x) = x^{n - m}Q(x) * x^{m}G(x) + x^{n - m - 1}R(x) * x^{m + 1}
\\
x^nF(x) \equiv x^{n - m}Q(x) * x^{m}G(x)\ \ (mod\ \ x^{n - m - 1})
\\
F_R(x) \equiv Q_R(x) * G_R(x)\ \ (mod\ \ x^{n - m - 1})
\\
Q_R(x) \equiv F_R(x) * G_R(x)^{-1}\ \ (mod\ \ x^{n - m - 1})
\]
\(F_R(x)\) 和 \(G_R(x)\) 都是已知的（预处理一下即可），\(Q(x)\) 就是 \(Q_R(x)\) 的系数倒过来。
计算出 \(Q_(x)\) 之后 \(R(x)\) 也就简单了：
\[R(x) = F(x) - Q(x)G(x)
\]
一定要注意边界啊啊啊啊！
#include <bits/stdc++.h>
#define ll long long
using namespace std;
namespace IO{
inline ll read(){
ll x = 0;
char ch = getchar();
while(!isdigit(ch)) ch = getchar();
while(isdigit(ch)) x = (x << 3) + (x << 1) + ch - '0', ch = getchar();
return x;
}
template <typename T> inline void write(T x){
if(x > 9) write(x / 10);
putchar(x % 10 + '0');
}
inline void print(ll a[], ll n){
for(int i = 0; i <= n; ++i) write(a[i]), putchar(' ');
puts("");
}
}
using namespace IO;
const ll N = 3e5 + 10;
const ll mod = 998244353;
const ll G = 3, Gi = 332748118;
ll n, m;
ll c[N];
ll f[N], g[N], fr[N], gr[N], ig[N];
ll q[N], r[N], p[N];
namespace NTT{
ll lim, len;
inline ll qpow(ll a, ll b){
ll res = 1;
while(b){
if(b & 1) res = res * a % mod;
a = a * a % mod, b >>= 1;
}
return res;
}
inline void get_rev(int n){
lim = 1, len = 0;
while(lim <= n) lim <<= 1, ++len;
for(int i = 0; i <= lim; ++i) p[i] = (p[i >> 1] >> 1) | ((i & 1) << (len - 1));
}
inline void ntt(ll A[], ll lim, int type){
for(int i = 0; i <= lim; ++i)
if(i < p[i]) swap(A[i], A[p[i]]);
for(int mid = 1; mid < lim; mid <<= 1){
ll Wn = qpow(type == 1 ? G : Gi, (mod - 1) / (mid << 1));
for(int i = 0; i < lim; i += (mid << 1)){
ll w = 1;
for(int j = 0; j < mid; ++j, w = w * Wn % mod){
ll x = A[i + j], y = w * A[i + j + mid] % mod;
A[i + j] = (x + y) % mod;
A[i + j + mid] = (x - y + mod) % mod;
}
}
}
if(type == 1) return;
ll inv = qpow(lim, mod - 2);
for(int i = 0; i <= lim; ++i) A[i] = A[i] * inv % mod;
}
inline void Mul(ll n, ll m, ll a[], ll b[]){
get_rev(n + m);
ntt(a, lim, 1), ntt(b, lim, 1);
for(int i = 0; i <= lim; ++i) a[i] = a[i] * b[i] % mod;
ntt(a, lim, -1), ntt(b, lim, -1);//--------------------------NTT回来，NTT回来，NTT回来！！！
}
inline void Inv(ll n, ll a[], ll b[]){
if(n == 0) return b[0] = qpow(a[0], mod - 2), void();
Inv(n >> 1, a, b);
get_rev(n << 1);
for(int i = 0; i <= n; ++i) c[i] = a[i];
for(int i = n + 1; i <= lim; ++i) c[i] = 0;
ntt(c, lim, 1), ntt(b, lim, 1);
for(int i = 0; i <= lim; ++i) b[i] = (2ll - c[i] * b[i] % mod + mod) * b[i] % mod;
ntt(b, lim, -1);
for(int i = n + 1; i <= lim; ++i) b[i] = 0;
}
}
using namespace NTT;
inline void Division(){
for(int i = n - m + 1; i <= m; ++i) gr[i] = 0;
Inv(n - m, gr, ig);
Mul(n, n - m, fr, ig);
for(int i = 0; i <= n - m; ++i) q[i] = fr[n - m - i];
Mul(n - m, m, g, q);
for(int i = 0; i < m; ++i) r[i] = (f[i] - g[i] + mod) % mod;
print(q, n - m);
print(r, m - 1);
}
signed main(){
n = read(), m = read();
for(int i = 0; i <= n; ++i) f[i] = read(), fr[n - i] = f[i];
for(int i = 0; i <= m; ++i) g[i] = read(), gr[m - i] = g[i];
Division();
return 0;
}
最后再来一发 P5273 【模板】多项式幂函数 (加强版) 的代码，怒调 3h（
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int mod = 998244353;
const int G = 3, Gi = 332748118;
const int N = 3e5 + 10;
ll n, k, flag, phik;
namespace IO{
inline ll read(){
ll x = 0;
char ch = getchar();
while(!isdigit(ch)) ch = getchar();
while(isdigit(ch)) x = (x << 3) + (x << 1) + ch - '0', ch = getchar();
return x;
}
inline ll readk(){
ll k = 0;
char ch = getchar();
while(!isdigit(ch)) ch = getchar();
while(isdigit(ch)){
if((k << 3) + (k << 1) + ch - '0' > n) flag = 1;
k = ((k << 3) + (k << 1) + ch - '0') % mod;
phik = ((phik << 3) + (phik << 1) + ch - '0') % (mod - 1);
ch = getchar();
}
return k;
}
template <typename T> inline void write(T x){
if(x > 9) write(x / 10);
putchar(x % 10 + '0');
}
}
using namespace IO;
ll A[N], B[N], c[N], lnb[N];
ll lowa[N], lna[N];
int p[N];
ll f1[N], f2[N], f[N], g[N];
inline ll qpow(ll a, int b){
ll res = 1;
while(b){
if(b & 1) res = res * a % mod;
a = a * a % mod, b >>= 1;
}
return res;
}
inline void clear(ll a[], int l, int r) {for(int i = l; i < r; ++i) a[i] = 0;}
inline void clone(ll a[], ll b[], int n) {for(int i = 0; i < n; ++i) a[i] = b[i];}
inline void print(ll a[], int n) {for(int i = 0; i < n; ++i) write(a[i]), putchar(' '); puts("");}
namespace NTT{
int lim, len;
inline void get_rev(int n){
lim = 1, len = 0;
while(lim < n) lim <<= 1, ++len;
for(int i = 0; i < lim; ++i) p[i] = (p[i >> 1] >> 1) | ((i & 1) << (len - 1));
}
inline void ntt(ll A[], int lim, int type){
for(int i = 0; i < lim; ++i)
if(i < p[i]) swap(A[i], A[p[i]]);
for(int mid = 1; mid < lim; mid <<= 1){
ll Wn = qpow(type == 1 ? G : Gi, (mod - 1) / (mid << 1));
for(int i = 0; i < lim; i += (mid << 1)){
ll w = 1;
for(int j = 0; j < mid; ++j, w = w * Wn % mod){
ll x = A[i + j], y = w * A[i + j + mid] % mod;
A[i + j] = (x + y) % mod;
A[i + j + mid] = (x - y + mod) % mod;
}
}
}
if(type == 1) return;
ll inv = qpow(lim, mod - 2);
for(int i = 0; i < lim; ++i) A[i] = A[i] * inv % mod;
}
inline void Mul(int n, int m, ll a[], ll b[]){
get_rev(n + m);
clear(A, 0, lim), clear(B, 0, lim);
clone(A, a, n), clone(B, b, m);
ntt(A, lim, 1), ntt(B, lim, 1);
for(int i = 0; i < lim; ++i) A[i] = A[i] * B[i] % mod;
ntt(A, lim, -1);
clone(a, A, lim);
}
inline void Inv(int n, ll a[], ll b[]){
if(n == 1) return b[0] = qpow(a[0], mod - 2), void();
Inv((n + 1) >> 1, a, b);
get_rev(n << 1);
clone(c, a, n), clear(c, n, lim);
ntt(b, lim, 1), ntt(c, lim, 1);
for(int i = 0; i < lim; ++i) b[i] = (2ll - c[i] * b[i] % mod + mod) * b[i] % mod;
ntt(b, lim, -1);
clear(b, n, lim);
}
}
using namespace NTT;
namespace Poly{
inline void Diff(int n, ll a[], ll b[]){
for(int i = 1; i < n; ++i) b[i - 1] = i * a[i] % mod;
b[n - 1] = 0;
}
inline void Integral(int n, ll a[], ll b[]){
for(int i = 1; i < n; ++i) b[i] = a[i - 1] * qpow(i, mod - 2) % mod;
b[0] = 0;
}
inline void Ln(int n, ll a[], ll b[]){
get_rev(n << 1);
clear(f1, 0, lim), clear(f2, 0, lim);
Diff(n, a, f1), Inv(n, a, f2);
Mul(n, n, f1, f2);
Integral(n, f1, b);
}
inline void Exp(int n, ll a[], ll b[]){
if(n == 1) return b[0] = 1, void();
Exp((n + 1) >> 1, a, b);
get_rev(n << 1);
Ln(n, b, lnb);
for(int i = 0; i < n; ++i) lnb[i] = (a[i] - lnb[i] + mod) % mod;
lnb[0]++;
Mul(n, n, b, lnb), clear(b, n, lim);
}
inline void Qpow(int n, int k, ll a[], ll b[], ll phik){
clear(b, 0, n);
int shift = 0;
while(!a[shift]) shift++;
if((ll)shift * k >= n) return;
n -= shift;
for(int i = 0; i < n; ++i) lowa[i] = a[i + shift];
int low0 = lowa[0], inv0 = qpow(low0, mod - 2);
for(int i = 0; i < n; ++i) lowa[i] = lowa[i] * inv0 % mod;
Ln(n, lowa, lna);
for(int i = 0; i < n; ++i) lna[i] = lna[i] * k % mod;
Exp(n, lna, b);
shift *= k;
for(int i = n + shift - 1; i >= shift; --i) b[i] = b[i - shift] * qpow(low0, phik) % mod;
clear(b, 0, shift);
}
}
using namespace Poly;
int main(){
n = read(), k = readk();
for(int i = 0; i < n; ++i) f[i] = read();
if(!flag || f[0]) Qpow(n, k, f, g, phik);
print(g, n);
return 0;
}
\[\_EOF\_
\]