讲下一个思路比较自然的基于自然数幂和的 \(O(n\log n)\) 且复杂度与 \(m\) 几乎无关的做法。
不难发现让我们计数的问题是保序回归 \(L_1\) 中一条链的情况。这个情况有一个简单的 slope-trick 做法：用堆维护斜率，每次 push 进去两个当前的数，然后 pop 出一个最大值。最终所有数的和减去堆中的数的和就是答案。
有一个来自 ARC128F 经典思维技巧：对于这类弹堆压堆还要求堆中元素和的计数问题，考虑转化成值域为 01 的问题。即活用 \(i=\sum_{x=1}^{\infin} [i\ge x]\) 的等式。将 \(\ge x\) 的数标成 \(1\)，剩余的数标成 \(0\)。那么如果原先堆中有 \(s\) 个 \(1\)，遇到一个 \(1\) 会变成一个 \(s+1\)，遇到一个 \(0\) 会变成 \(\max(s-1,0)\)。这就是一个格路游走问题，容易验证最终 \(s\) 的值等于路程中 \(1\) 的个数减去路程中 \(0\) 的个数再加上如果没有对 \(0\) 取 \(\max\) 的情况下，\(s\) 在整个过程中的最小值。
发现唯一难算的就是最后一部分的“最小值之和”。考虑继续运用上面的思维技巧，继续活用等式 \(i=\sum_{x=1}^{\infin} [i\ge x]\)，将 \(s\) 的游走过程看作 \((0,0)\) 到 \((n,\times)\) 的游走过程，那么最小值之和（的相反数）可以拆成触碰 \(y=-t\ (t>0)\) 这条线的方案数之和。钦定触碰一条线的格路游走就是我们熟悉的反射容斥。我们设原来有 \(p\) 个 \(1\)，那么原先终点在 \((n,2p-n)\)，如果 \(2p-n>t\) 则经过那条线的方案数等于到终点 \((n,-2t-2p+n)\) 的方案数。
所以对于一个固定的 \(p\)，我们需要对以下东西求和：
\[S_p=\sum_{t=1}^{\infin} {n\choose \min(n-t-p,p)} \]这个可以预处理 \({n\choose \times}\) 这一行的组合数前缀和简单算出来。
现在考虑对于 \(x\) 你需要依次带入 \(x=1,2,\dots,m\)，不妨把结果看成一个关于 \(x\) 的多项式，则这个多项式实际上是：
\[F(x)=\sum_{p=0}^n x^p(m-x)^{n-p}S_p \]展开 \((m-x)^{n-p}\)，容易发现可以用一遍卷积求出 \(F(x)\)。
最后我们只需要解决 \(\sum_{i=1}^m F(i)\)，相当于要对一个固定的 \(n\) 求出 \(k=1,2,\dots,n\) 的自然数幂和 \(S_k(n)=\sum_{i=1}^n i^k\)，可以用多项式求逆求出伯努利数，然后卷一次得到自然数幂和。于是我们就做到复杂度 \(O(n\log n)\)。
写给自己看的备忘笔记：伯努利数感觉还是直接生成
#include <cstdio>
#include <vector>
#include <cassert>
#include <algorithm>
// headers
struct poly{/** my poly template **/};
int n,m;
const int N=100103;
int fac[N],fiv[N];
int arr[N],pre[N];
int coe[N],pw[N];
inline int C(int a,int b){return (ll)fac[a]*fiv[b]%P*fiv[a-b]%P;}
void calc(int lim){
pw[0]=1;
for(int i=1;i<=lim;++i) pw[i]=(ll)pw[i-1]*m%P;
fac[0]=1;
for(int i=1;i<=lim;++i) fac[i]=(ll)fac[i-1]*i%P;
fiv[lim]=qp(fac[lim]);
for(int i=lim;i;--i) fiv[i-1]=(ll)fiv[i]*i%P;
}
int main(){
n=read();m=read();
calc(n+3);int res=(ll)n*(m-1)%P*pw[n]%P;
if(res&1) res+=P;
res>>=1;
pre[0]=1;
for(int i=1;i<=n;++i){
pre[i]=pre[i-1]+C(n,i);
if(pre[i]>=P) pre[i]-=P;
}
poly F(n+1),G(n+1);
for(int i=0;i<=n;++i){
if(2*i<=n) F[i]=(pre[i-1]+(ll)C(n,i)*(n-i*2))%P;
else F[i]=pre[n-1-i];
F[i]=(ll)F[i]*fac[n-i]%P;
}
for(int i=0;i<=n;++i) if(i&1) G[i]=fiv[i];else G[i]=P-fiv[i];
G=F*G;
for(int i=0;i<=n;++i) coe[i]=(ll)G[i]*pw[n-i]%P*fiv[n-i]%P;
F.f.resize(n+2);
G.f.resize(n+2);
for(int i=0;i<=n+1;++i) F[i]=fiv[i+1];
F=F.inv(n+2);
for(int i=0,tt=1;i<=n+1;++i){
G[i]=(ll)tt*fiv[i]%P;
tt=(ll)tt*(m+1)%P;
}
G=F*G;
for(int i=0;i<=n;++i) res=(res+(ll)coe[i]*(G[i+1]-F[i+1]+P)%P*fac[i])%P;
res-=coe[0];
if(res<0) res+=P;
printf("%d\n",res);
return 0;
}