闲话 22.8.19
闲话
终于到Aster了啊！
算起来我把Aster()当做枚举前向星出边的宏定义已经半年了
Aster
(Aster, Aster, My Lover, and Death, "PINK KISS" in my heart)
「…309へ集合…」
「...在309集合...」
ドクターが撃たれたんだ
医生被枪击干掉了
アンチトロイの契約なんて遂行されちゃいなかった
那签过的防木马病毒条约并没有起到任何作用
100年前のシティ?ポップがネオン街のテーマソング
100年前流行的CityPop可是这条霓虹街的主题曲
すっかりNPCと私たちだけになった
然而这里只剩下NPC和我们了
ネオン管の夢色
霓虹管的梦色
パパとママに気付いてほしくて
多希望爸爸妈妈能察觉到
心配ないよアスター
没什么可担心的，Aster
私元気なセブンティーンさ
十七岁朝气蓬勃的我可是
この街に囚われたネオン
被这座城市囚禁住了的霓虹
つまりはファンタジー おとぎの話
总而言之一切只是幻想，是童话故事
ノンフィクションなのがたまにキズさ
就算不是虚构的偶尔也还是会有伤口
まぼろしのファンタジー 飽きたらおしまい
如幻影一般的幻想若厌倦了就结束了
夢の見過ぎに御注意を
小心不要太沉溺于梦境
(Aster, Aster, My Lover, and Death, "PINK KISS" in my heart)
作り物のワイヤーフレーム
人造品的线框图
特別なプレゼンテイション
特别的简述报告
永遠の甘い夢
永恒的美梦
人は何にも見えなくなった
人们变得什么都看不见了
仕事を持ってママと出会ってお家を買って心を買って
好好工作和妈妈相遇，买了房子买了心
真っ白な身体に植え付けられただけ
然后只需把一切植入到纯白的身体里
不安定な瞳は
不安分的瞳孔
パパとママの笑顔も忘れる
连爸爸妈妈的笑容都忘记了
幻想だったアスター
这幻想的都市Aster
全部どうやら紛い物だ
看来一切都只是伪造品
誰一人息を呑むネオン
没有人为之屏息的霓虹
ひとりのファンタジー 涙はいらない
只身一人的幻想不需要掉眼泪
センセーションだけが私なのだ
我是个只有精神感官的存在而已
夢色のファンタジー パパは気付けない
梦色的幻想爸爸并没有注意到
これで最後の夜になる
这将会成为最后的夜晚
壊れたファンタジー
崩坏了的幻想
壊れたファンタジー
崩坏了的幻想
さよならファンタジー 孤独な夢も
再见了幻想，孤独的梦也
ノンフィクションだから寂しくなる
因为不是虚构的才会感到如此寂寞吧
夢もケも同じようで偽物
梦境和现实一样都是假的
黄泉の国まで御注意を
到达黄泉之前请多注意
(Aster, Aster, My Lover, and Death, "PINK KISS" in my heart)
霓虹色的世界啊……
有一种淡淡的忧伤
Aster有紫菀的意思 而Aster的pv里主颜色就是霓虹紫
德莉莎也有紫菀这个名字的装甲嗯
那明天就放《银河录》好了！
春卷饭早期歌里我心目中的前三！
是宇宙的世界呢
在和jjdw聊今天该推什么番
jjdw：成神之日！
我：您是在什么精神状态下说出这句话的
jjdw：这番确实挺好的，展现了高超的麻将技巧
我：和棒球技巧
所以今天的推番不是《成神之日》
大家一致认为不要去看《成神之日》这番
又在和jjdw聊今天该推什么番
jjdw：珈百璃的堕落
我：您终于阳间了一会
所以今天的推番是《珈百璃的堕落》！
又是无厘头搞怪日常！
【模板】多项式
多项式是啥？
一个 \(n\) 阶多项式，是形如
\[F(x) = \sum_{i = 0}^n f_i x^i
\]
的式子，其中 \(f_i\) 是系数，\(x\) 是变量，\(n\) 也可以叫次数。
记 \(F(x)\) 的第 \(i\) 项系数为 \(f_i\)，其余类似记法同理。
我们以下式定义 \(F\) 与 \(G\) 的卷积 \(H\)：
\[h_i = \sum_{j = 0}^i f_j g_{i-j}
\]
查看定义可知这玩意就是 \(H(x) = F(x) \times G(x)\)。
多项式乘法
本处使用NTT，具体不讲。详见gtm的博客
typedef long long ll;
const int N = 3e6 + 10, mod = 998244353, g = 3, invg = 332748118;
const int siz_ll = sizeof(ll);
int btrs[N]; // butterfly_transform
inline int initrs(int k) { // 每次带入一个k作为目标多项式的长度
int limit = 1;
while (limit < k) limit <<= 1;
for (register int i = 0 ; i < limit; i ++)
btrs[i] = ((btrs[i >> 1] >> 1) | ((i & 1) ? limit >> 1 : 0)); // 线性递推 当然换写法也可以
return limit;
}
inline ll qp(ll a, ll b) {
ll ret = 1;
while (b) {
if (b & 1) ret = ret * a % mod;
a = a * a % mod;
b >>= 1;
} return ret;
}
int L, w[2][1<<19];
inline int __INITILIZE__UNIT__ROOT__() { // 预处理单位根，听说会快
L = 1<<19;
w[0][0] = w[1][0] = 1;
int wn = qp(g, (mod-1) / L);
for (register int i = 1; i < L; i++) w[0][L - i] = w[1][i] = 1ll * w[1][i - 1] * wn % mod;
return 1;
} int __INITIALIZER__ = __INITILIZE__UNIT__ROOT__();
inline void NTT (ll * A, const int limit, const int type) { // 显然的ntt，对长为limit的A系数进行处理，type=1为dft，type=0为idft
for (register int i = 0; i < limit; i++) if (i < btrs[i]) swap(A[i], A[btrs[i]]);
for (register int mid = 1; mid < limit; mid <<= 1) {
for (register int i = L / (mid<<1), j = 0; j < limit; j += (mid << 1)) {
for (register int k = 0; k < mid; k++) {
int x = A[j + k], y = w[type][i * k] * A[j + k + mid] % mod;
A[j + k] = (x + y) % mod;
A[j + k + mid] = (x - y + mod) % mod;
}
}
} if (type == 1) return;
ll inv = qp(limit, mod - 2);
for (register int i = 0; i < limit; i++) A[i] = A[i] * inv % mod;
}
inline void Mul(ll * F, ll * G, ll * Ret, int n, int m) { // Ret = F * G, len(F) = n, len(G) = m
static ll __F[N], __G[N];
int __n = m + n;
int limit = initrs(__n);
memcpy(__F, F, n * siz_ll);
memcpy(__G, G, m * siz_ll);
memset(__F + n, 0, siz_ll * (limit - n));
memset(__G + m, 0, siz_ll * (limit - m));
NTT(__F, limit, 1), NTT(__G, limit, 1);
for (register int i = 0; i < limit; i++) Ret[i] = __F[i] * __G[i] % mod;
NTT(Ret, limit, 0);
}
多项式乘法逆
给定 \(n-1\) 次多项式 \(F\) 的系数。输出多项式 \(G\) 的系数，满足 \(F(x) \times G(x) \equiv 1\ ( \text{mod} \ x^n)\)。
我们进行一波式子的推导。
假设我们已知 \(F(x) \times H(x) \equiv 1\ ( \text{mod} \ x^{\lceil \frac n2 \rceil})\)，又显然有 \(F(x) \times G(x) \equiv 1\ ( \text{mod} \ x^{\lceil \frac n2 \rceil})\)，则
\[\begin{aligned}
F(x) \times G(x) - F(x) \times H(x) &\equiv 0 & \pmod{x^{\lceil \frac n2 \rceil}} \\
F(x) \times (G(x) - H(x)) &\equiv 0 & \pmod{x^{\lceil \frac n2 \rceil}} \\
G(x) - H(x) &\equiv 0 & \pmod{x^{\lceil \frac n2 \rceil}} \\
\end{aligned}
\]
两边同时平方。由于 \(G(x) - H(x)\) 在 \(\text{mod} \ x^{\lceil \frac n2 \rceil}\) 意义下为 0，因此其第 0 项至第 \(\lceil \frac x2 \rceil - 1\) 项都为0。平方后显然有对于任意n以内项的系数都是由 \(\lceil \frac n2 \rceil - 1\) 以内项生成的，因此这些项均为0。我们有
\[\begin{aligned}
(G(x) - H(x))^2 & \equiv 0 &\pmod{ x^n}\\
G^2 (x) + H^2 (x)) - 2 G(x) H(x) &\equiv 0 &\pmod{ x^n}\\
F(x)G^2 (x) + F(x)H^2 (x)) - 2 F(x) G(x) H(x) &\equiv 0 & \pmod{ x^n}\\
G (x) + F(x)H^2 (x) - 2 H(x) &\equiv 0 &\pmod{x^n}\\
G (x) &\equiv 2 H(x) - F(x)H^2 (x) & \pmod{ x^n}
\end{aligned}
\]
因此我们可以通过递归的方式首先求得 \(\text{mod} \ x^{\lceil \frac n2 \rceil}\) 意义下的 \(G\)，随后带回，多项式乘法算出当前解。当 \(n = 1\) 时有边界为 \(g_0 = f_0 ^ {p-2}\)。
inline void Inv(ll * A, ll * Iv, int deg) {
if (deg == 1) { Iv[0] = qp(A[0], mod-2); return; }
static ll tmp[N];
Inv(A, Iv, (deg+1) >> 1);
int limit = initrs(deg<<1);
memset(tmp + deg, 0, siz_ll * (limit - deg));
memcpy(tmp, A, deg * siz_ll);
NTT(tmp, limit, 1), NTT(Iv, limit, 1);
for (register int i = 0; i < limit; i++) {
Iv[i] = (2 - tmp[i] * Iv[i] % mod + mod) * Iv[i] % mod;
} NTT(Iv, limit, 0);
for (register int i = deg; i < limit; i ++) Iv[i] = 0;
}
多项式除法
给定 \(n,m\) 阶多项式 \(F\)，\(G\) 的系数。输出多项式 \(Q\)，\(R\) 的系数，满足 \(F(x) = G(x) \times Q(x) + R(x)\)。
定义 \(F^{R}\) 为把 \(F\) 的系数翻转后得到的多项式。容易得到若 \(F\) 为 \(k\) 阶则 \(F^{R} = F(\frac 1x) x^k\)。
我们将 \(x = \frac 1x\) 带入原式，有
\[\begin{aligned}
F(\frac 1x) &= G(\frac 1x) \times Q(\frac 1x) + R(\frac 1x) &\\
x^n\times F(\frac 1x) &= x^n\times G(\frac 1x) \times Q(\frac 1x) + x^n\times R(\frac 1x) &\\
F(\frac 1x) x^n &= G(\frac 1x) x^m \times Q(\frac 1x) x^{n - m} + \times R(\frac 1x) x^{m-1} \times x^{n-m+1} &\\
F^R(x) &= G^R(x) \times Q^R(x)+ \times R^R(x) \times x^{n-m+1} &\\
F^R(x)&\equiv G^R(x) \times Q^R(x) &\pmod{x^{n-m+1}} \\
Q^R(x) &\equiv \frac {F^R(x)} {G^R(x)} &\pmod{x^{n-m+1}} \\
\end{aligned}
\]
由于 \(Q\) 是 \(n-m\) 次的，因此在\(\pmod {x^{n-m+1}}\) 意义下不会影响求解。
求解完后回带入原式求出 \(R(x) = F(x) - G(x) \times Q(x)\) 求解即可。
进行了部分工业化修饰，意义已填入注释。
inline void Div(ll * F, ll * G, ll * Quot, ll * Remd, int n, int m) { // F = G * Quot + Remd, len(F) = n, len(G) = m
static ll __F[N], __G[N], __Inv_G[N], __Q[N], __Q_G[N], __Res_G[N];
int __n = (n<<2);
memset(__F, 0, __n * siz_ll);
memset(__G, 0, __n * siz_ll);
memset(__Inv_G, 0, __n * siz_ll);
memset(__Q_G, 0, __n * siz_ll);
memset(__Q, 0, __n * siz_ll);
memset(__Res_G, 0, __n * siz_ll);
// for (register int i = 0; i < __n; i++) _F[i] = _G[i] = __Inv_G[i] = __Q_G[i] = __Q[i] = __Res_G[i] = 0;
memcpy(__F, F, n * siz_ll); // for (register int i = 0; i < n; i++) _F[i] = F[i];
memcpy(__G, G, m * siz_ll);
memcpy(__Res_G, G, m * siz_ll); // for (register int i = 0; i < m; i++) _G[i] = __Res_G[i] = G[i];
reverse(__F, __F + n); reverse(__G, __G + m);
memset(__G + n - m + 1, 0, siz_ll * (__n - n + m - 1)); // for (register int i = n - m + 1; i < __n; i++) _G[i] = 0;
Inv(__G, __Inv_G, n - m + 1);
Mul(__F, __Inv_G, __Q, n, n - m + 1);
memset(__Q + n - m + 1, 0, siz_ll * (__n - n + m - 1)); // for (register int i = n - m + 1; i < __n; i++) __Q[i] = 0;
reverse(__Q, __Q + n - m + 1);
memcpy(Quot, __Q, siz_ll * (n - m + 1)); // for (register int i = 0; i < n - m + 1; i++) Quot[i] = __Q[i];
Mul(__Q, __Res_G, __Q_G, n - m + 1, m);
for (register int i = 0; i < m-1; i++) Remd[i] = (F[i] - __Q_G[i] + mod) % mod;
}
多项式ln
给定 \(n-1\) 次多项式 \(F\) 的系数。输出多项式 \(G\) 的系数，满足 \(G(x) \equiv \ln(F(x))( \text{mod} \ x^n)\)。
两边求导，我们有
\[G'(x) \equiv \frac{F'(x)}{F(x)} \pmod {x^n}
\]
直接模拟多项式求导和多项式不定积分即可
这个小二基础吧
inline void Deri(ll * F, ll * Ret, int n) { // ret = F', len(F) = n
for (register int i = 1; i < n; i++) Ret[i-1] = 1ll * i * F[i] % mod;
Ret[n - 1] = 0;
}
inline void Intg(ll * F, ll * Ret, int n) { // ret' = F, len(F) = n
for (register int i = 1; i < n; i++) Ret[i] = 1ll * F[i-1] * qp(i, mod-2) % mod;
Ret[0] = 0;
}
inline void Ln(ll * F, ll * Ret, int n) { // ret = ln(F), len(F) = n
static ll __Dri[N], __Inv_F[N];
int limit = initrs(n<<1);
memset(__Dri, 0, siz_ll * limit);
memset(__Inv_F, 0, siz_ll * limit);
Deri(F, __Dri, n); Inv(F, __Inv_F, n);
memset(__Dri + n, 0, siz_ll * (limit - n));
memset(__Inv_F + n, 0, siz_ll * (limit - n));
Mul(__Dri, __Inv_F, __Dri, limit, limit);
Intg(__Dri, Ret, n);
memset(Ret + n, 0, siz_ll * (limit - n));
}
多项式exp
给定 \(n-1\) 次多项式 \(F\) 的系数。输出多项式 \(G\) 的系数，满足 \(G(x) \equiv e^{F(x)}( \text{mod} \ x^n)\)。
两边取自然对数，我们有 \(\ln(G(x)) - F(x) \equiv 0 \pmod{x ^n}\)
我们设 \(T(G(x)) = \ln G(x) - F(x)\)，于是我们就需要找 \(T\) 的零点。
牛爵爷已经给了我们一个优秀的方法求解如上问题。我们设目前解为 \(G_0(x)\)，则可以有如下迭代式：
\[G(x) \leftarrow G_0(x) - \frac{T(G_0(x))}{T'(G_0(x))}
\]
每一次迭代都可以使精度翻倍。这种方法被称为牛顿迭代法。
回到原题，我们视 \(F(x)\) 为常数，对 \(T\) 求导得 \(\frac{G'(x)} {G(x)}\)
带入牛迭式子可以得到
\[G(x)= G_0(x) - \frac{\ln(G_0(x)) - F(x) }{\frac{G_0'(x)} {G_0(x)}} = \frac{G_0(x) (1 - \ln (G_0(x)) + F(x)}{G_0'(x)}
\]
迭代求解即可。
inline void Exp(ll * F, ll * Ret, int n) { // F = ln(Ret), len(F) = n
if (n == 1) { Ret[0] = 1; return; }
Exp(F, Ret, (n + 1) >> 1);
static ll __Ln_R[N], __Res_F[N];
memset(__Ln_R, 0, siz_ll * ((n<<1) + 1));
memcpy(__Res_F, F, siz_ll * ((n<<1) + 1));
Ln(Ret, __Ln_R, n);
int limit = initrs(n<<1);
memset(__Res_F + n, 0, siz_ll * (limit - n));
for (register int i = 0; i < limit; i++)
__Res_F[i] = (__Res_F[i] - __Ln_R[i] + mod) % mod;
__Res_F[0]++;
Mul(Ret, __Res_F, Ret, limit, limit);
memset(Ret + n, 0, siz_ll * (limit - n));
}
先写到这里（
最后放个全的板子：
多项式1/4家桶
typedef long long ll;
const int N = 3e6 + 10, mod = 998244353, g = 3, invg = 332748118;
const int siz_ll = sizeof(ll);
int n, m;
ll a[N], b[N], s;
int btrs[N]; // butterfly_transform
inline int initrs(int k) {
int limit = 1;
while (limit < k) limit <<= 1;
for (register int i = 0 ; i < limit; i ++)
btrs[i] = ((btrs[i >> 1] >> 1) | ((i & 1) ? limit >> 1 : 0));
return limit;
}
inline ll qp(ll a, ll b) {
ll ret = 1;
while (b) {
if (b & 1) ret = ret * a % mod;
a = a * a % mod;
b >>= 1;
} return ret;
}
int L, w[2][1<<19];
inline int __INITILIZE__UNIT__ROOT__() {
L = 1<<19;
w[0][0] = w[1][0] = 1;
int wn = qp(g, (mod-1) / L);
for (register int i = 1; i < L; i++) w[0][L - i] = w[1][i] = 1ll * w[1][i - 1] * wn % mod;
return 1;
} int __INITIALIZER__ = __INITILIZE__UNIT__ROOT__();
inline void NTT (ll * A, const int limit, const int type) {
for (register int i = 0; i < limit; i++) if (i < btrs[i]) swap(A[i], A[btrs[i]]);
for (register int mid = 1; mid < limit; mid <<= 1) {
for (register int i = L / (mid<<1), j = 0; j < limit; j += (mid << 1)) {
for (register int k = 0; k < mid; k++) {
int x = A[j + k], y = w[type][i * k] * A[j + k + mid] % mod;
A[j + k] = (x + y) % mod;
A[j + k + mid] = (x - y + mod) % mod;
}
}
} if (type == 1) return;
ll inv = qp(limit, mod - 2);
for (register int i = 0; i < limit; i++) A[i] = A[i] * inv % mod;
}
inline void Inv(ll * A, ll * Iv, int deg) {
if (deg == 1) { Iv[0] = qp(A[0], mod-2); return; }
static ll tmp[N];
Inv(A, Iv, (deg+1) >> 1);
int limit = initrs(deg<<1);
memset(tmp + deg, 0, siz_ll * (limit - deg));
memcpy(tmp, A, deg * siz_ll);
NTT(tmp, limit, 1), NTT(Iv, limit, 1);
for (register int i = 0; i < limit; i++) {
Iv[i] = (2 - tmp[i] * Iv[i] % mod + mod) * Iv[i] % mod;
} NTT(Iv, limit, 0);
for (register int i = deg; i < limit; i ++) Iv[i] = 0;
}
inline void Mul(ll * F, ll * G, ll * Ret, int n, int m) { // Ret = F * G, len(F) = n, len(G) = m
static ll __F[N], __G[N];
int __n = m + n;
int limit = initrs(__n);
memcpy(__F, F, n * siz_ll);
memcpy(__G, G, m * siz_ll);
memset(__F + n, 0, siz_ll * (limit - n));
memset(__G + m, 0, siz_ll * (limit - m));
NTT(__F, limit, 1), NTT(__G, limit, 1);
for (register int i = 0; i < limit; i++) Ret[i] = __F[i] * __G[i] % mod;
NTT(Ret, limit, 0);
} inline void mul(ll * F, ll val, int n) { for (register int i = 0; i < n; i++) F[i] = F[i] * val % mod; }
inline void Div(ll * F, ll * G, ll * Quot, ll * Remd, int n, int m) { // F = G * Quot + Remd, len(F) = n, len(G) = m
static ll __F[N], __G[N], __Inv_G[N], __Q[N], __Q_G[N], __Res_G[N];
int __n = (n<<2);
memset(__F, 0, __n * siz_ll);
memset(__G, 0, __n * siz_ll);
memset(__Inv_G, 0, __n * siz_ll);
memset(__Q_G, 0, __n * siz_ll);
memset(__Q, 0, __n * siz_ll);
memset(__Res_G, 0, __n * siz_ll); // for (register int i = 0; i < __n; i++) _F[i] = _G[i] = __Inv_G[i] = __Q_G[i] = __Q[i] = __Res_G[i] = 0;
memcpy(__F, F, n * siz_ll); // for (register int i = 0; i < n; i++) _F[i] = F[i];
memcpy(__G, G, m * siz_ll);
memcpy(__Res_G, G, m * siz_ll); // for (register int i = 0; i < m; i++) _G[i] = __Res_G[i] = G[i];
reverse(__F, __F + n); reverse(__G, __G + m);
memset(__G + n - m + 1, 0, siz_ll * (__n - n + m - 1)); // for (register int i = n - m + 1; i < __n; i++) _G[i] = 0;
Inv(__G, __Inv_G, n - m + 1);
Mul(__F, __Inv_G, __Q, n, n - m + 1);
memset(__Q + n - m + 1, 0, siz_ll * (__n - n + m - 1)); // for (register int i = n - m + 1; i < __n; i++) __Q[i] = 0;
reverse(__Q, __Q + n - m + 1);
memcpy(Quot, __Q, siz_ll * (n - m + 1)); // for (register int i = 0; i < n - m + 1; i++) Quot[i] = __Q[i];
Mul(__Q, __Res_G, __Q_G, n - m + 1, m);
for (register int i = 0; i < m-1; i++) Remd[i] = (F[i] - __Q_G[i] + mod) % mod;
}
inline void Deri(ll * F, ll * Ret, int n) { // ret = F', len(F) = n
for (register int i = 1; i < n; i++) Ret[i-1] = 1ll * i * F[i] % mod;
Ret[n - 1] = 0;
}
inline void Intg(ll * F, ll * Ret, int n) { // ret' = F, len(F) = n
for (register int i = 1; i < n; i++) Ret[i] = 1ll * F[i-1] * qp(i, mod-2) % mod;
Ret[0] = 0;
}
inline void Ln(ll * F, ll * Ret, int n) { // ret = ln(F), len(F) = n
static ll __Dri[N], __Inv_F[N];
int limit = initrs(n<<1);
memset(__Dri, 0, siz_ll * limit);
memset(__Inv_F, 0, siz_ll * limit);
Deri(F, __Dri, n); Inv(F, __Inv_F, n);
memset(__Dri + n, 0, siz_ll * (limit - n));
memset(__Inv_F + n, 0, siz_ll * (limit - n));
Mul(__Dri, __Inv_F, __Dri, limit, limit);
Intg(__Dri, Ret, n);
memset(Ret + n, 0, siz_ll * (limit - n));
}
inline void Exp(ll * F, ll * Ret, int n) { // F = ln(Ret), len(F) = n
if (n == 1) { Ret[0] = 1; return; }
Exp(F, Ret, (n + 1) >> 1);
static ll __Ln_R[N], __Res_F[N];
memset(__Ln_R, 0, siz_ll * ((n<<1) + 1));
memcpy(__Res_F, F, siz_ll * ((n<<1) + 1));
Ln(Ret, __Ln_R, n);
int limit = initrs(n<<1);
memset(__Res_F + n, 0, siz_ll * (limit - n));
for (register int i = 0; i < limit; i++)
__Res_F[i] = (__Res_F[i] - __Ln_R[i] + mod) % mod;
__Res_F[0]++;
Mul(Ret, __Res_F, Ret, limit, limit);
memset(Ret + n, 0, siz_ll * (limit - n));
}
好这是封装了的
#include <bits/stdc++.h>
#define rep(i, a, b) for (register int(i) = (a); (i) <= (b); ++(i))
#define pre(i, a, b) for (register int(i) = (a); (i) >= (b); --(i))
using namespace std;
typedef long long ll;
const int N = 3e6 + 10, mod = 998244353, g = 3, invg = 332748118;
const int siz_ll = sizeof(ll);
int n, m;
ll a[N], b[N], s;
int btrs[N]; // butterfly_transform
inline int initrs(int k) {
int limit = 1;
while (limit < k) limit <<= 1;
for (register int i = 0 ; i < limit; i ++)
btrs[i] = ((btrs[i >> 1] >> 1) | ((i & 1) ? limit >> 1 : 0));
return limit;
}
inline ll qp(ll a, ll b) {
ll ret = 1;
while (b) {
if (b & 1) ret = ret * a % mod;
a = a * a % mod;
b >>= 1;
} return ret;
} inline ll inv(ll a) { return qp(a, mod-2); }
int L, w[2][1<<19];
inline int __INITILIZE__UNIT__ROOT__() {
L = 1<<19;
w[0][0] = w[1][0] = 1;
int wn = qp(g, (mod-1) / L);
for (register int i = 1; i < L; i++) w[0][L - i] = w[1][i] = 1ll * w[1][i - 1] * wn % mod;
return 1;
} int __INITIALIZER__ = __INITILIZE__UNIT__ROOT__();
struct poly {
vector <ll> f;
ll operator [] (const int & pos) const { return f[pos]; }
ll & operator [] (const int & pos) { return f[pos]; }
int deg() {return f.size(); }
int deg() const {return f.size(); }
void Set(int n) { f.resize(n); }
void Adjust() { while (f.size() > 1 and f.back() == 0) f.pop_back(); }
void Reverse() { reverse(f.begin(), f.end()); }
void deri() { for (register int i = 0; i + 1 < (int)f.size(); ++i) f[i] = 1ll * f[i+1] * (i+1) % mod; f.back() = 0; }
void intg() { for (register int i = f.size() - 1; i >= 1; --i) f[i] = 1ll * f[i-1] * inv(i) % mod; f[0] = 0; }
void scan(int n = -1) { if (n < 0) get(n); Set(n); for (register int i = 0; i < n; i++) get(f[i]); }
void print() { for (ll x : f) printf("%lld ", x); }
poly Deri() {
poly ret; ret.Set(deg());
for (int i = 0; i + 1 < deg(); ++i) {
ret[i] = 1ll * f[i+1] * (i+1) % mod;
} return ret;
} poly Intg() {
poly ret; ret.Set(deg());
for (int i = deg() - 1; i >= 1; --i) {
ret[i] = 1ll * f[i-1] * inv(i) % mod;
} return ret;
}
inline void NTT (const int lim, const int type) {
Set(lim);
for (register int i = 0; i < lim; i++) if (i < btrs[i]) swap(f[i], f[btrs[i]]);
for (register int mid = 1; mid < lim; mid <<= 1) {
for (register int i = L / (mid<<1), j = 0; j < lim; j += (mid << 1)) {
for (register int k = 0; k < mid; k++) {
int x = f[j + k], y = w[type][i * k] * f[j + k + mid] % mod;
f[j + k] = (x + y) % mod;
f[j + k + mid] = (x - y + mod) % mod;
}
}
} if (type == 1) return;
ll inv = qp(lim, mod - 2);
for (register int i = 0; i < lim; i++) f[i] = f[i] * inv % mod;
}
friend poly operator + (const poly & x, const poly & y) {
poly ret; ret.Set(max(x.deg(), y.deg()));
for (register int i = 0; i < x.deg(); ++i) ret[i] = x[i];
for (register int i = 0; i < y.deg(); ++i) ret[i] = (ret[i] + y[i]) % mod;
return ret;
} void operator += (const poly & x) {
Set(max(deg(), x.deg()));
for (register int i = 0; i < x.deg(); ++i) f[i] = (f[i] + x[i]) % mod;
}
friend poly operator - (const poly & x, const poly & y) {
poly ret; ret.Set(max(x.deg(), y.deg()));
for (register int i = 0; i < x.deg(); ++i) ret[i] = x[i];
for (register int i = 0; i < y.deg(); ++i) ret[i] = (ret[i] - y[i] + mod) % mod;
return ret;
} void operator -= (const poly & x) {
Set(max(deg(), x.deg()));
for (register int i = 0; i < x.deg(); ++i) f[i] = (f[i] - x[i] + mod) % mod;
}
friend poly operator * (const poly & x, const poly & y) {
poly ret, A = x, B = y;
int limit = initrs(A.deg() + B.deg() - 1);
A.NTT(limit, 1), B.NTT(limit, 1); ret.Set(limit);
for (register int i = 0; i < limit; i++) ret[i] = 1ll * A[i] * B[i] % mod;
ret.NTT(limit, 0); ret.Adjust();
return ret;
} void operator *= (const poly & x) {
poly A = x;
int limit = initrs(deg() + A.deg() - 1);
A.NTT(limit, 1); NTT(limit, 1);
for (register int i = 0; i < limit; i++) f[i] = 1ll * A[i] * f[i] % mod;
NTT(limit, 0); Adjust();
}
friend poly operator / (const poly & x, const poly & y) {
poly A = x, B = y;
A.Reverse(), B.Reverse(); B.Set(x.deg() - y.deg() + 1);
B = B.Inv(); B *= A;
B.Set(x.deg() - y.deg() + 1); B.Reverse();
return B;
}
friend pair<poly,poly> operator % (const poly & x, const poly & y) {
poly A = x, B = y;
A.Reverse(), B.Reverse(); B.Set(x.deg() - y.deg() + 1);
B = B.Inv(); B *= A;
B.Set(x.deg() - y.deg() + 1); B.Reverse();
poly R = x - B * y; R.Set(y.deg() - 1);
return {B, R};
}
poly Inv() {
poly ret; ret.Set(1);
if (f.empty()) return ret;
ret[0] = inv(f[0]); poly A, B;
for (register int len = 2, limit; len < (deg() << 1); len <<= 1) {
A.f.assign((*this).f.begin(), (*this).f.begin() + min(len, deg()));
B = ret; B.Set(min(len, deg()));
limit = initrs(A.deg() + B.deg() - 1);
A.NTT(limit, 1); B.NTT(limit, 1);
ret.Set(limit);
for (register int i = 0; i < limit; i++) ret[i] = (2 - A[i] * B[i] % mod + mod) * B[i] % mod;
ret.NTT(limit, 0);
ret.Set(len);
} ret.Set(deg()); return ret;
}
poly Ln() {
poly ret = Deri(), Iv = Inv();
ret *= Iv; ret.intg(); ret.Set(deg());
return ret;
}
poly Exp() {
poly ret; ret.Set(1); ret[0] = 1;
poly A, B;
for (register int len = 2; len < (deg() << 1); len <<= 1) {
ret.Set(len); A = ret.Ln();
B.f.assign((*this).f.begin(), (*this).f.begin() + min(len, deg()));
B -= A; B[0]++; ret *= B;
ret.Set(len);
} ret.Set(deg()); return ret;
}
} F, G;
以下是博客签名，与正文无关。
请按如下方式引用此页：
本文作者 joke3579，原文链接：https://www.cnblogs.com/joke3579/p/chitchat220819.html。
遵循 CC BY-NC-SA 4.0 协议。
请读者尽量不要在评论区发布与博客内文完全无关的评论，视情况可能删除。
