讲下一个思路比较自然的基于自然数幂和的 $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$，我们需要对以下东西求和：

这个可以预处理 ${n\choose \times}$ 这一行的组合数前缀和简单算出来。

现在考虑对于 $x$ 你需要依次带入 $x=1,2,\dots,m$，不妨把结果看成一个关于 $x$ 的多项式，则这个多项式实际上是：

展开 $(m-x)^{n-p}$，容易发现可以用一遍卷积求出 $F$。

最后我们只需要解决 $\sum_{i=1}^m F(i)$，相当于要对一个固定的 $n$ 求出 $k=0,1,2,\dots,n$ 的自然数幂和 $S_k(n)=\sum_{i=0}^n i^k$，可以用多项式求逆求出伯努利数，然后卷一次得到自然数幂和。于是我们就做到复杂度 $O(n\log n)$。

关于伯努利数求自然数幂和：。

写给自己看的备忘笔记：伯努利数感觉还是直接生成函数定义更好理解。其 EGF 为 $B=\frac{x}{e^x-1}$。可以对 $\sum_{i=0}^{\infin} \frac{1}{(i+1)!}x^i$ 多项式求逆单 $\log$ 求得。

对于自然数幂和，我们考虑研究其关于 $k$ 的 EGF：

则：

#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;
}