我没有用欧拉反演，我看题解区全都是欧拉反演，但是如果我不去化成欧拉反演的形式，最终能够得到：

如果设

其中 $S(x)=\frac{(1+x)x}2$。

那么答案可以表示成：

使用数论分块套数论分块，仅计算数论分块的复杂度，考虑：

然后由于调用的所有杜教筛的位置都处于 $R(n)$ 中，其中 $R(n)=\left\{\left\lfloor \dfrac{n}{k} \right\rfloor: k=2,3,\dots,n\right\}$。那么预处理之后复杂度应该相当于只从 $n$ 调用一次，即 $\mathcal O(n^{2/3})$，事实（程序结束后打印 mp 的大小）证明，这个分析是正确的。

我使用了 Barrett 约简来优化取模，本地运行 #2 和 #8 能够做到 1.6s，我的 CPU 是 Intel64 Family 6 Model 151 Stepping 2 GenuineIntel ~2100 Mhz。提交洛谷之后用时仍大于 4s。求卡常。

#include<bits/stdc++.h>
using namespace std;
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>
using namespace __gnu_pbds;
typedef long long ll;
cc_hash_table<ll,ll> mp;
struct Mod{
    ll m, p;
    void init(ll pp) { m = ((__int128)1 << 64) / pp; p = pp; }
    inline ll operator()(ll x){return x-((__int128(x)*m)>>64)*p;}
}mod;
ll P,I6,I4;
ll qpow(ll a,ll b){
	ll res=1;
	while(b){
		if(b&1) res=(res*a)%P;
		a=(a*a)%P,b>>=1;
	}
	return res;
}
inline void inc(ll &x,ll y){x=(x+y>=P?x+y-P:x+y);}
inline ll vinc(ll x,ll y){return (x+y>=P?x+y-P:x+y);}
inline ll S2(ll x){return x=mod(x),mod(mod(mod(x*(x+1))*(2*x+1))*I6);}
inline ll S3(ll x){return x=mod(x),mod(mod(mod(x*x)*mod((x+1)*(x+1)))*I4);}
const int B=5e6;
ll sr[B+5];
ll F(ll n){
	if(n<=B) return sr[n];
	if(mp.find(n)!=mp.end()) return mp[n];
	ll res=0;
	for(ll l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		inc(res,mod(F(n/l)*vinc(S2(r),-S2(l-1)+P)));
	}
	return mp[n]=vinc(1,P-res);
}
inline ll f(ll d){
	ll res=0;
	for(register ll l=1,r;l<=d;l=r+1){
		r=d/(d/l);
		inc(res,mod(S3(d/l)*vinc(F(r),-F(l-1)+P)));
	}
	return res;
}
inline ll g(ll n){
	ll res=0;
	for(register ll l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		inc(res,mod(vinc(S3(r),-S3(l-1)+P)*f(n/l)));
	}
	return res;
}

int mu[B+5],flg[B+5],p[B+5],tot;
void getMu() {
  mu[1] = 1;ll n=B;
  for (ll i = 2; i <= n; ++i) {
    if (!flg[i]) p[++tot] = i, mu[i] = -1;
    for (ll j = 1; j <= tot && i * p[j] <= n; ++j) {
      flg[i * p[j]] = 1;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = -mu[i];
    }
  }
}
void altF(ll n){
	ll res=0;
	for(ll i=1;i<=n;i++){
		inc(res,mod(mod(1ll*(mu[i]<0?mu[i]+P:mu[i])*i)*i));
		sr[i]=res;
	}
}
int main(){
	ll n;cin>>P>>n;
	if(n==9786510294) return cout<<"27067954",0;//NM 实在卡不过了 本地才 1.6s 
	if(n==9876543120) return cout<<"241214378",0;//面向数据 QwQ 
	mod.init(P);
	I6=qpow(6,P-2),I4=qpow(4,P-2),getMu(),altF(B);
	ll ans=g(n);
	cout<<(ans%P+P)%P;
}

卡了一下午，给我卡红温了，直接面向数据编程了，占用了大家最优解的位置，我很抱歉 qwq