多项式运算
现在大概有：
-
多项式乘法逆
-
多项式开根
-
多项式带余除法
-
多项式牛顿迭代
-
多项式 \(\exp\)
-
多项式 \(\ln\)
-
多项式快速幂
多项式复合函数一类的高端内容后面再补。
很多地方的严谨证明因为水平不足略过，等之后数学水平提高再回来修正。
一些有用的资料：
NTT与多项式全家桶 - command_block 的博客
前言
- 运算定义
这里的多项式运算都可以用加法和乘法（加法卷积）定义。
- 界
通常对于无限项的多项式，只考虑多项式的前 \(n\) 项，也就是在 \(\pmod {x^n}\) 意义下运算。
相应地有运算性质：
-
\(F(x) \bmod x^n + G(x) \bmod x^n \equiv (F(x) + G(x)) \bmod x^n \pmod {x^n}\)
-
\(F(x) \bmod x^n \cdot G(x) \bmod x^n \equiv (F(x) \cdot G(x)) \bmod x^n \pmod {x^n}\)
所以可以在运算的过程中限定多项式的界，从而避免高次项。
当然这个也可以用来卡常。
- 循环卷积
普通的卷积是线性卷积，另外还有周期卷积和循环卷积。
循环卷积的定义是：
\((F * G)(x) = \sum\limits_{k = 0}^m x^k \sum\limits_{i + j \bmod m = k} F_i G_j\)
更直观的形式是 \((F * G)(x) = \sum\limits_{i = 0}^{m} \sum\limits_{j = 0}^m F_i G_j x^{(i + j) \bmod m}\)
封装好的 NTT 模板：
void calc_rev(int k) { for (int i = 1; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }
ll qpow(ll base, ll power)
{
ll res = 1;
while (power)
{
if (power & 1) res = res * base % mod;
base = base * base % mod;
power >>= 1;
}
return res;
}
void NTT(ll *A, int n)
{
calc_rev(n);
for (int i = 1; i < n; i++)
if (rev[i] > i) swap(A[i], A[rev[i]]);
for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
{
ll wn = qpow(g, (mod - 1) / len);
wp[0] = 1;
for (int i = 1; i <= len; i++) wp[i] = wp[i - 1] * wn % mod;
for (int l = 0, r = len - 1; r <= n; l += len, r += len)
{
int w = 0;
for (int p = l; p < l + m; p++, w++)
{
ll x = A[p], y = wp[w] * A[p + m] % mod;
A[p] = (x + y) % mod, A[p + m] = (x - y + mod) % mod;
}
}
}
}
void INTT(ll *A, int n)
{
NTT(A, n);
reverse(A + 1, A + n);
int inv = qpow(n, mod - 2);
for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * inv % mod;
}
多项式乘法逆
求多项式在模意义下的逆元。
朴素做法
这个虽然有一种复杂度为 \(O(n^2)\) 的递推做法，但是有更优的 \(O(n \log n)\) NTT 倍增做法。
当界为 \(\pmod {x^1}\) （只有常数项）时直接求自然数的逆元就行。
反之，假设已经有 \(R_* (x) \equiv F(x)^{-1} \pmod {x^\frac{n}{2}}\)，考虑倍增求出 \(R(x) \equiv F(x)^{-1} \pmod {x^n}\).
因为 \(R(x) \equiv R_*(x) \pmod {x^{\frac{n}{2}}}\)，所以 \(R(x) - R_*(x) \equiv 0 \pmod {x^{\frac{n}{2}}}\).
考虑对两边平方得到 \((R(x) - R_*(x))^2 \equiv 0 \pmod {x^n}\)
展开得到 \(R(x)^2 - 2 R(x) R_*(x) + R_*(x)^2 \equiv 0 \pmod {x^n}\)
两边同乘 \(F(x)\) 得到 \(R(x) - 2 R_*(x) + R_*(x)^2 F(x) \equiv 0 \pmod {x^n}\)
整理即 \(R(x) \equiv 2 R_*(x) - R_*(x)^2 F(x) \pmod {x^n}\)
时间复杂度为 \(T(n) = T(\frac{n}{2}) + O(n \log n) = O(n \log n)\)，下面的倍增法同理。
优化
NTT 代入的是原根，它具有单位根的性质，所以实际上 NTT 求的是循环卷积。
观察倍增的式子：\(R(x) \equiv 2 R_*(x) - R_*(x)^2 F(x) \pmod {x^n}\).
需要用到多项式乘法的项只有 \(R_*(x)^2 F(x)\)，其中 \(R_*(x)\) 的次数是 \(\frac{n}{2}\)，\(F(x)\) 的次数是 \(n\)，并且我们只关心 \(R(x)\) 的后 \(\frac{n}{2}\) 项。
可以考虑先算 \(R_*(x) F(x)\)，此时由 \(R_*(x)\) 的定义得其前 \(\frac{n}{2}\) 项中只有第一项为 \(1\)，其余都为 \(0\). 于是循环卷积只会影响前 \(\frac{n}{2}\) 项。
然后再考虑乘上 \(R_*(x)\)，正好还是影响前 \(\frac{n}{2}\) 项。
void invp(ll *f, ll *r, int n)
{
int k = 1;
while (k < n) k <<= 1;
r[0] = qpow(f[0], mod - 2);
for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
{
for (int i = 0; i < len; i++) Rt[i] = r[i], Ft[i] = f[i];
NTT(Ft, len), NTT(Rt, len);
for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
INTT(Rt, len);
for (int i = 0; i < m; i++) Rt[i] = 0; Rt[0] = 1;
for (int i = 0; i < len; i++) Ft[i] = r[i];
NTT(Ft, len), NTT(Rt, len);
for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
INTT(Rt, len);
for (int i = m; i < len; i++) r[i] = (r[i] * 2ll - Rt[i] + mod) % mod;
}
memset(Ft, 0, k * sizeof(ll));
memset(Rt, 0, k * sizeof(ll));
for (int i = n; i < k; i++) r[i] = 0;
}
多项式开根
还是考虑倍增。
令 \(B(x)^2 \equiv A(x) \pmod {x^n}\)，假设现在已知 \(B_*(x) \equiv A(x) \pmod {x^{\frac{n}{2}}}\)，考虑倍增求出 \(B(x)\).
同求逆有 \(B(x) - B_*(x) \equiv 0 \pmod {x^{\frac{n}{2}}}\).
平方再展开得 \(B(x)^2 - 2 B(x) B_*(x) + B_*(x)^2 \equiv 0 \pmod {x^n}\).
又因为 \(B(x)^2 \equiv A(x) \pmod {x^n}\)，代入得：
\(A(x) - 2 B(x) B_*(x) + B_*(x)^2 \equiv 0 \pmod {x^n}\).
整理得 \(B(x) \equiv \frac{\frac{A(x)}{B_*(x)} + B_*(x)}{2}\)
考虑上全家桶就行。
注意常数项是二次剩余，当 \(A(x)\) 常数项为 \(1\) 时显然 \(B(x)\) 常数项也为 \(1\).
void sqrtp(ll *f, ll *s, int n)
{
s[0] = mod_sqrt(f[0]);
if (s[0] == -1) return;
int k = 1;
while (k < (n << 1)) k <<= 1;
ll inv2 = qpow(2, mod - 2, mod);
for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
{
invp(s, rt, m);
for (int i = 0; i < m; i++) ft[i] = f[i];
NTT(ft, len), NTT(rt, len);
for (int i = 0; i < len; i++) rt[i] = rt[i] * ft[i] % mod;
INTT(rt, len);
for (int i = 0; i < m; i++) rt[i] = (rt[i] + s[i]) % mod;
for (int i = 0; i < len; i++) s[i] = rt[i] * inv2 % mod, rt[i] = ft[i] = 0;
}
for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
for (int i = n; i < k; i++) s[i] = 0;
}
多项式带余除法
整除的情况比较好做，考虑把余数去掉。
舍弃高次项可以直接用界做，所以考虑把整个多项式反转（对称的项互换系数）。
定义多项式 \(F(x)\) 的反转 \(F_R(x) = x^n F(\frac{1}{x})\).
考虑题目所求：\(F(x) = Q(x) G(x) + R(x)\)，其中 \(F\) 是 \(n\) 次已知多项式，\(G\) 是 \(m\) 次已知多项式，\(Q(x)\) 是 \(n - m\) 次未知多项式，\(Q\) 是 \(m - 1\) 次未知多项式（不够用 \(0\) 补足）。
考虑换元，根据上式有 \(F(\frac{1}{x}) = Q(\frac{1}{x}) G(\frac{1}{x}) + R(\frac{1}{x})\).
两边同乘 \(x^n\) 得 \(x^n F(\frac{1}{x}) = x^n Q(\frac{1}{x}) G(\frac{1}{x}) + x^n R(\frac{1}{x})\)
又因为 \(x^n F(\frac{1}{x}) = F_R(x), x^n Q(\frac{1}{x}) G(\frac{1}{x}) = Q_R(x) G_R(x), x^n R(\frac{1}{x}) = x^{n - m + 1} R_R(x)\)，代入原式得：
\(F_R(x) = Q_R(x) G_R(x) + x^{n - m + 1} R_R(x)\)
注意到 \(R(x)\) 是 \(n - m\) 次多项式，不妨给上式加上界：
\(F_R(x) = Q_R(x) G_R(x) + x^{n - m + 1} R_R(x) \pmod {x^{n - m + 1}}\)
于是可以求出 \(Q_R(x)\)，从而求出 \(Q(x)\) 和 \(R(x)\).
void divp(ll *f, ll *g, ll *quo, ll *rem, int n, int m)
{
int len = n - m + 1;
for (int i = 0; i < min(len, m); i++) rt[i] = g[m - i - 1];
invp(rt, ft, len);
for (int i = 0; i < len; i++) rt[i] = f[n - i - 1];
int k = 1;
while (k < (len << 1)) k <<= 1;
NTT(rt, k), NTT(ft, k);
for (int i = 0; i < k; i++) rt[i] = rt[i] * ft[i] % mod;
INTT(rt, k);
for (int i = 0; i < len; i++) quo[i] = rt[len - i - 1];
for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
for (int i = 0; i < m; i++) rt[i] = g[i];
for (int i = 0; i < len; i++) ft[i] = quo[i];
k = 1;
while (k < (m + len)) k <<= 1;
NTT(rt, k), NTT(ft, k);
for (int i = 0; i < k; i++) rem[i] = rt[i] * ft[i] % mod;
INTT(rem, k);
for (int i = 0; i < m - 1; i++) rem[i] = (f[i] - rem[i] + mod) % mod;
for (int i = m - 1; i < k; i++) rem[i] = 0;
for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
}
多项式牛顿迭代
需要简单的微积分知识，只需要知道求导和积分的算法即可。
首先定义多项式的求导：\(F^{\prime}(x) = \sum\limits_{i = 1}^n a_{i - 1} i x^i\)
void diffp(ll *f, ll *der, int n)
{
for (int i = 1; i < n; i++) der[i - 1] = f[i] * i % mod;
der[n - 1] = 0;
}
以及多项式的积分：\(\int F(x) dx = \sum\limits_{i = 0}^{n - 1} \frac{a_i x^{i + 1}}{i}\)
void intep(ll *f, ll *inte, int n)
{
for (int i = 1; i < n; i++)
if (!inv[i]) inv[i] = (i == 1 ? 1 : inv[mod % i] * (mod - mod / i) % mod);
for (int i = 1; i < n; i++) inte[i] = f[i - 1] * inv[i] % mod;
inte[0] = 0;
}
也就是对每一项分别求导再求和，符合求导和积分的经典形式。
同时多项式的求导和积分也是互为逆运算。
牛顿迭代：已知函数 \(G(F(x)) = 0\)，求 \(F(x) \bmod x^n\)，这里 \(G(x)\) 通常为根据需要构造的函数。
结论：令 \(F_*(x) \equiv F(x) \pmod {x^{\frac{n}{2}}}\)，\(F(x) \equiv F_*(x) - \frac{G(F_*(x))}{G^{\prime}(F_*(x))} \pmod {x^n}\).
具体证明见 多项式计数杂谈 - command_block 的博客
多项式对数函数
求：\(\ln(A(x)) \equiv B(x) \pmod {x^n}\).
两边求导得 \(ln^{\prime}(A(x)) A^{\prime}(x) \equiv B^{\prime}(x) \pmod {x^n}\)
又因为 \(ln^{\prime}(A(x)) = \frac{1}{A(x)}\)，所以得：
\(\frac{A^{\prime}(x)}{A(x)} \equiv B^{\prime}(x) \pmod {x^n}\)
重新积分得 \(B(x) \equiv \frac{A^{\prime}(x)}{A(x)} dx\)
void lnp(ll *f, ll *ln, int n)
{
diffp(f, ft, n), invp(f, rt, n);
int k = 1;
while (k < (n << 1)) k <<= 1;
NTT(ft, k), NTT(rt, k);
for (int i = 0; i < k; i++) ft[i] = ft[i] * rt[i] % mod;
INTT(ft, k);
intep(ft, ln, n);
for (int i = 0; i < k; i++) ft[i] = rt[i] = 0;
}
多项式指数函数
求 \(B(x) \equiv \exp(A(x)) \pmod {x^n}\).
考虑牛顿迭代。设 \(G(F(x)) = \ln(F(x)) - A(x)\)，那么 \(G(B(x)) \equiv 0 \pmod {x^n}\).
根据牛顿迭代结论有 \(B(x) \equiv B_*(x) - \frac{G(B_*(x))}{\frac{1}{B_*(x)}} = B_*(x) - G(B_*(x)) B_*(x) \pmod {x^n}\).
整理得 \(B(x) \equiv (1 - \ln(B_*(x) + A(x)) B_*(x) \pmod {x^n}\).
继续考虑倍增求。
边界为 \(B_0 = 1\)（\(\exp(0) = 1\)）
void expp(ll *f, ll *exp, int n)
{
int k = 1;
while (k < n) k <<= 1;
exp[0] = 1;
for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
{
for (int i = 0; i < len; i++) Fn[i] = exp[i];
lnp(Fn, Rn, len);
for (int i = 0; i < len; i++) Rn[i] = (f[i] - Rn[i] + mod) % mod;
Rn[0] = (Rn[0] + 1) % mod;
NTT(Fn, len << 1), NTT(Rn, len << 1);
for (int i = 0; i < (len << 1); i++) Fn[i] = Fn[i] * Rn[i] % mod;
INTT(Fn, len << 1);
for (int i = 0; i < len; i++) exp[i] = Fn[i];
}
for (int i = 0; i < (k << 1); i++) Fn[i] = Rn[i] = 0;
for (int i = n; i < k; i++) exp[i] = 0;
}
多项式快速幂
当 \(A_0 = 1\) 时，根据 \((A(x))^k = \exp(k \ln(A(x)))\) 直接做。
void powp(ll *f, ll *pw, int n, int k)
{
lnp(f, fn, n);
for (int i = 0; i < n; i++) fn[i] = fn[i] * k % mod;
expp(fn, pw, n);
for (int i = 0; i < n; i++) fn[i] = 0;
}
\(A_0 \neq 1\) 时求不到 \(\exp\)，考虑把 \(A_0\) 化成 \(1\).
直接给整个多项式乘上 \(\frac{1}{A_0}\)，最后再乘回 \(A_0^k\).
如果遇到系数为 \(0\) 的项也没关系，直接把整个多项式移一下 /xk
不过因为是多项式乘法，最后 \(0\) 的个数会变成原来的 \(k\) 倍，可以直接特判。
void powp2(ll *f, ll *pw, int n, int k, int k2)
{
int pos = 0;
while ((pos < n) && (f[pos] == 0)) pos++;
if (1ll * pos * k >= n) return;
int len0 = pos * k, len = n - len0;
ll inv = qpow(f[pos], mod - 2, mod), val = qpow(f[pos], k2, mod);
for (int i = 0; i < len; i++) fn[i] = f[pos + i] * inv % mod;
lnp(fn, rn, len);
for (int i = 0; i < len; i++) rn[i] = rn[i] * k % mod, fn[i] = 0;
expp(rn, fn, len);
for (int i = 0; i < len0; i++) pw[i] = 0;
for (int i = len0; i < n; i++) pw[i] = fn[i - len0] * val % mod;
for (int i = 0; i < len; i++) fn[i] = rn[i] = 0;
}
全家桶模板：
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <map>
#include <algorithm>
using namespace std;
typedef long long ll;
const int sz = 3e6 + 5;
const int mod = 998244353;
const int g = 3;
int n;
int rev[sz];
ll F[sz], G[sz], inv[sz], wp[sz];
ll Ft[sz], Rt[sz], ft[sz], rt[sz], Fn[sz], Rn[sz], fn[sz], rn[sz];
ll quo[sz], rem[sz];
map<ll, ll> mp;
void readk(int &k, int &k2)
{
k = k2 = 0;
char ch = getchar();
while ((ch < '0') || (ch > '9')) ch = getchar();
while ((ch >= '0') && (ch <= '9'))
{
ll tmp = 10ll * k + ch - '0', tmp2 = 10ll * k2 + ch - '0';
k = tmp % mod, k2 = tmp2 % (mod - 1);
if (tmp >= mod) k += mod;
if (tmp2 >= mod - 1) k2 += mod - 1;
ch = getchar();
}
}
void calc_rev(int k) { for (int i = 1; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }
ll qpow(ll base, ll power, ll mod)
{
ll res = 1;
while (power)
{
if (power & 1) res = res * base % mod;
base = base * base % mod;
power >>= 1;
}
return res;
}
ll gcd(ll a, ll b) { return (!b ? a : gcd(b, a % b)); }
ll bsgs(ll a, ll n, ll p, ll ad)
{
ll m = ceil(sqrt(p));
for (ll i = 0, t = n; i <= m; i++, t = t * a % p) mp[t] = i;
for (ll i = 0, pw = qpow(a, m, p), t = ad; i <= m; i++, t = t * pw % p)
if (mp.count(t) && (i * m - mp[t] >= 0)) return i * m - mp[t];
return -1;
}
ll exbsgs(ll a, ll n, ll p)
{
a %= p, n %= p;
if ((n == 1) || (p == 1)) return 0;
if (a == 0) return (n == 0 ? 1 : -1);
ll cnt = 0, ad = 1, d;
while ((d = gcd(a, p)) != 1)
{
if (n % d != 0) return -1;
cnt++, n /= d, p /= d, ad = ad * (a / d) % p;
if (ad == n) return cnt;
}
ll ans = bsgs(a, n, p, ad);
if (ans == -1) return -1;
return ans + cnt;
}
ll mod_sqrt(ll x)
{
ll base = exbsgs(g, x, mod);
if ((base == -1) || (base & 1)) return -1;
return qpow(g, base / 2, mod);
}
void NTT(ll *A, int n)
{
calc_rev(n);
for (int i = 1; i < n; i++)
if (rev[i] > i) swap(A[i], A[rev[i]]);
for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
{
ll wn = qpow(g, (mod - 1) / len, mod);
wp[0] = 1;
for (int i = 1; i <= len; i++) wp[i] = wp[i - 1] * wn % mod;
for (int l = 0, r = len - 1; r <= n; l += len, r += len)
{
int w = 0;
for (int p = l; p < l + m; p++, w++)
{
ll x = A[p], y = wp[w] * A[p + m] % mod;
A[p] = (x + y) % mod, A[p + m] = (x - y + mod) % mod;
}
}
}
}
void INTT(ll *A, int n)
{
NTT(A, n);
reverse(A + 1, A + n);
int inv = qpow(n, mod - 2, mod);
for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * inv % mod;
}
void invp(ll *f, ll *r, int n)
{
int k = 1;
while (k < n) k <<= 1;
r[0] = qpow(f[0], mod - 2, mod);
for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
{
for (int i = 0; i < len; i++) Rt[i] = r[i], Ft[i] = f[i];
NTT(Ft, len), NTT(Rt, len);
for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
INTT(Rt, len);
for (int i = 0; i < m; i++) Rt[i] = 0; Rt[0] = 1;
for (int i = 0; i < len; i++) Ft[i] = r[i];
NTT(Ft, len), NTT(Rt, len);
for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
INTT(Rt, len);
for (int i = m; i < len; i++) r[i] = (r[i] * 2ll - Rt[i] + mod) % mod;
}
memset(Ft, 0, k * sizeof(ll));
memset(Rt, 0, k * sizeof(ll));
for (int i = n; i < k; i++) r[i] = 0;
}
void divp(ll *f, ll *g, ll *quo, ll *rem, int n, int m)
{
int len = n - m + 1;
for (int i = 0; i < min(len, m); i++) rt[i] = g[m - i - 1];
invp(rt, ft, len);
for (int i = 0; i < len; i++) rt[i] = f[n - i - 1];
int k = 1;
while (k < (len << 1)) k <<= 1;
NTT(rt, k), NTT(ft, k);
for (int i = 0; i < k; i++) rt[i] = rt[i] * ft[i] % mod;
INTT(rt, k);
for (int i = 0; i < len; i++) quo[i] = rt[len - i - 1];
for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
for (int i = 0; i < m; i++) rt[i] = g[i];
for (int i = 0; i < len; i++) ft[i] = quo[i];
k = 1;
while (k < (m + len)) k <<= 1;
NTT(rt, k), NTT(ft, k);
for (int i = 0; i < k; i++) rem[i] = rt[i] * ft[i] % mod;
INTT(rem, k);
for (int i = 0; i < m - 1; i++) rem[i] = (f[i] - rem[i] + mod) % mod;
for (int i = m - 1; i < k; i++) rem[i] = 0;
for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
}
void sqrtp(ll *f, ll *s, int n)
{
s[0] = mod_sqrt(f[0]);
if (s[0] == -1) return;
int k = 1;
while (k < (n << 1)) k <<= 1;
ll inv2 = qpow(2, mod - 2, mod);
for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
{
invp(s, rt, m);
for (int i = 0; i < m; i++) ft[i] = f[i];
NTT(ft, len), NTT(rt, len);
for (int i = 0; i < len; i++) rt[i] = rt[i] * ft[i] % mod;
INTT(rt, len);
for (int i = 0; i < m; i++) rt[i] = (rt[i] + s[i]) % mod;
for (int i = 0; i < len; i++) s[i] = rt[i] * inv2 % mod, rt[i] = ft[i] = 0;
}
for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
for (int i = n; i < k; i++) s[i] = 0;
}
void diffp(ll *f, ll *der, int n)
{
for (int i = 1; i < n; i++) der[i - 1] = f[i] * i % mod;
der[n - 1] = 0;
}
void intep(ll *f, ll *inte, int n)
{
for (int i = 1; i < n; i++)
if (!inv[i]) inv[i] = (i == 1 ? 1 : inv[mod % i] * (mod - mod / i) % mod);
for (int i = 1; i < n; i++) inte[i] = f[i - 1] * inv[i] % mod;
inte[0] = 0;
}
void lnp(ll *f, ll *ln, int n)
{
diffp(f, ft, n), invp(f, rt, n);
int k = 1;
while (k < (n << 1)) k <<= 1;
NTT(ft, k), NTT(rt, k);
for (int i = 0; i < k; i++) ft[i] = ft[i] * rt[i] % mod;
INTT(ft, k);
intep(ft, ln, n);
for (int i = 0; i < k; i++) ft[i] = rt[i] = 0;
}
void expp(ll *f, ll *exp, int n)
{
int k = 1;
while (k < n) k <<= 1;
exp[0] = 1;
for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
{
for (int i = 0; i < len; i++) Fn[i] = exp[i];
lnp(Fn, Rn, len);
for (int i = 0; i < len; i++) Rn[i] = (f[i] - Rn[i] + mod) % mod;
Rn[0] = (Rn[0] + 1) % mod;
NTT(Fn, len << 1), NTT(Rn, len << 1);
for (int i = 0; i < (len << 1); i++) Fn[i] = Fn[i] * Rn[i] % mod;
INTT(Fn, len << 1);
for (int i = 0; i < len; i++) exp[i] = Fn[i];
}
for (int i = 0; i < (k << 1); i++) Fn[i] = Rn[i] = 0;
for (int i = n; i < k; i++) exp[i] = 0;
}
void powp(ll *f, ll *pw, int n, int k)
{
lnp(f, fn, n);
for (int i = 0; i < n; i++) fn[i] = fn[i] * k % mod;
expp(fn, pw, n);
for (int i = 0; i < n; i++) fn[i] = 0;
}
void powp2(ll *f, ll *pw, int n, int k, int k2)
{
int pos = 0;
while ((pos < n) && (f[pos] == 0)) pos++;
if (1ll * pos * k >= n) return;
int len0 = pos * k, len = n - len0;
ll inv = qpow(f[pos], mod - 2, mod), val = qpow(f[pos], k2, mod);
for (int i = 0; i < len; i++) fn[i] = f[pos + i] * inv % mod;
lnp(fn, rn, len);
for (int i = 0; i < len; i++) rn[i] = rn[i] * k % mod, fn[i] = 0;
expp(rn, fn, len);
for (int i = 0; i < len0; i++) pw[i] = 0;
for (int i = len0; i < n; i++) pw[i] = fn[i - len0] * val % mod;
for (int i = 0; i < len; i++) fn[i] = rn[i] = 0;
}
int main()
{
return 0;
}
