Abnormal Permutation Pairs (hard version)
两个限制：字典序小、逆序对大，一个显然的想法就是确保一对关系，统计另一对关系。
确保哪一对呢？我们想了想，决定确保字典序小，因为字典序是可以贪心的。
具体而言，我们考虑两个排列自第 \(i\) 位开始出现了不同。这样子，我们便将两个排列各自划成两段，即 \([1,i-1]\) 与 \([i,n]\)。则两排列的第一段是相同的。
考虑借鉴 CDQ 分治的思想，将逆序对分为三类，即首段中、末段中、两段间。因为两排列首段相同，所以首段中的逆序对数相等；因为两排列中每一段的组成元素相同，所以两段间的逆序对数亦相等；有区别的只有末端中的逆序对数。
于是，我们考虑枚举这个 \(i\)，则前 \(i-1\) 位就有 \(n^{\underline{i-1}}\) 种取值方案。
我们还剩下 \(n-i+1\) 个互不相同的数。因为逆序对仅与大小关系相关，所以我们完全可以将其映射作一个长度为 \(n-i+1\) 的排列而不改变逆序对数。
现在考虑枚举第 \(i\) 位两个排列分别填了什么。设排列 \(p\) 填了 \(j\)（实际上是剩余数中第 \(j\) 小的数），\(q\) 填了 \(k\)（同理，实质是第 \(k\) 小的数），则要满足字典序限制则应有 \(j<k\)。\(j\) 对后面元素贡献了 \(j-1\) 个逆序对，\(k\) 对后面元素贡献了 \(k-1\) 个逆序对，则现在 \(p\) 比 \(q\) 少了 \(k-j\) 个逆序对。
因为后 \(n-i\) 位以后就可以随便填了，所以我们就设一个 \(f_{i,j}\) 表示长度为 \(i\) 的排列出现 \(j\) 个逆序对的方案数。枚举数字 \(i\) 填哪里就可以做到 \(n^4\)，用前缀和/差分优化就能做到 \(n^3\)，过于基础不再赘述。
现在回到 \(p\) 填 \(j\)、\(q\) 填 \(k\) 的情形。则，\(p\) 要保证后 \(n-i\) 位中逆序对数至少比 \(q\) 多 \(k-j+1\) 个才能满足逆序对的限制。
考虑枚举 \(i,j,k\) 以及 \(p,q\) 在后 \(n-i\) 位中逆序对数分别是 \(J,K\)，可以做到垃圾的 \(n^7\) 复杂度；
但注意到我们对于同样的 \(J-K\) 只关心 \(k-j+1\)，于是我们仅枚举 \(J,K\)，再预处理出此时合法的 \(j,k\) 对数，就能做到 \(n^5\)；
假如再观察到 \(j-k+1\) 最大只到 \(n-i\)，这样 \(J-K\) 在大于 \(n-i\) 后，对于再之前的所有 \(K\) 来说合法的 \(j,k\) 对数便会一直相等，这样枚举 \(K\) 时就只用枚举 \(J\) 之前的 \(n-i\) 个位置，之前的用一个前缀和就能解决，这样就能做到 \(n^4\)；
优化到正解 \(n^3\) 的做法是考虑 \(J\) 增加时，对于不同的 \(K\) 来说 \(j-k\) 数对的增量是什么；然后发现这一增量在 \(J-K=2\) 时是 \(n-i,J-K=3\) 时是 \(n-i-1,J-K=k\) 时是 \(n-i-k+2\)。放到平面直角坐标系上就是一个三角形的形式，可以通过预处理 \(f_{i,j}\) 数组的前缀和以及 \(j\times f_{i,j}\) 的前缀和来 \(O(1)\) 求出。
细节很多，详见代码。
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=510;
const int M=125100;
int n,mod,res;
int com[N],f[M],fac[N],s[M],t[M],C[N][N],d[N];
signed main(){
cin>>n>>mod;
for(int i=1;i<=n;i++) com[i]=i*(i-1)/2;
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%mod;
for(int i=0;i<=n;i++) C[i][0]=1;
for(int i=1;i<=n;i++) for(int j=1;j<=i;j++) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
f[0]=1;
for(int i=0;i<n;i++){
for(int j=1;j<=com[i];j++) (f[j]+=f[j-1])%=mod;
f[com[i]+1]=0;
for(int j=0;j<=i;j++) d[j]=0;
for(int J=0;J<=i;J++) for(int K=J+1;K<=i;K++) d[K-J]++;
for(int j=1;j<=i;j++) d[j]+=d[j-1];
for(int j=0;j<=com[i];j++){
s[j]=f[j];
t[j]=1ll*f[j]*j%mod;
}
for(int j=1;j<=com[i];j++){
(s[j]+=s[j-1])%=mod;
(t[j]+=t[j-1])%=mod;
}
for(int j=2,k=0;j<=com[i];j++){
if(j>i){
(res+=1ll*d[i]*fac[n-i-1]%mod*C[n][n-i-1]%mod*f[j]%mod*s[j-i-1]%mod)%=mod;
(k+=mod-1ll*f[j-i-1]*d[i]%mod)%=mod;
}
(k+=(0ll+t[j-2]+mod-1ll*(j-i-2)*s[j-2]%mod)%mod)%=mod;
if(j-1>i) (k+=(2ll*mod-t[j-i-2]+1ll*(j-i-2)*s[j-i-2]%mod)%mod)%=mod;
(res+=1ll*C[n][n-i-1]*fac[n-i-1]%mod*f[j]%mod*k%mod)%=mod;
}
for(int j=com[i];j>=0;j--) (f[j+i+1]+=mod-f[j])%=mod;
}
cout<<res;
return 0;
}