$sumlimits_{p=1}^{n}p^psumlimits_{d=1}^{lfloorfrac{n}{p} floor}mu(d)d^{2p}sumlimits_{i=1}^{lfloorfrac{n}{pd} floor}i^psumlimits_{j=1}^{lfloorfrac{n}{pd} floor}j^p$
枚举p，发现预处理自然数幂前缀和的复杂度是$O(n/p)$的，后面整个式子也就可以在$O(n/p)$内求出。
复杂度为调和级数$O(nlog n)$
1 #include<cstdio> 2 #include<algorithm> 3 #define rep(i,l,r) for (int i=(l); i<=(r); i++) 4 using namespace std; 5 6 const int N=500010,mod=1e9+7; 7 int n,m,ans,tot,sm[N],pw[N],mu[N],p[N],b[N]; 8 9 int ksm(int a,int b){ 10 int res=1; 11 for (; b; a=1ll*a*a%mod,b>>=1) 12 if (b & 1) res=1ll*res*a%mod; 13 return res; 14 } 15 16 void init(int n){ 17 mu[1]=1; 18 rep(i,2,n){ 19 if (!b[i]) p[++tot]=i,mu[i]=-1; 20 for (int j=1; j<=tot && p[j]*i<=n; j++){ 21 b[p[j]*i]=1; 22 if (i%p[j]==0) { mu[p[j]*i]=0; break; } 23 else mu[p[j]*i]=-mu[i]; 24 } 25 } 26 } 27 28 int main(){ 29 freopen("bzoj3561.in","r",stdin); 30 freopen("bzoj3561.out","w",stdout); 31 scanf("%d%d",&n,&m); 32 if (n>m) swap(n,m); init(n); 33 rep(i,1,m) pw[i]=1; 34 rep(i,1,n){ 35 int res=0; sm[0]=0; 36 rep(d,1,m/i) pw[d]=1ll*pw[d]*d%mod,sm[d]=(sm[d-1]+pw[d])%mod; 37 rep(d,1,n/i) res=(res+1ll*mu[d]*pw[d]%mod*pw[d]%mod*sm[n/i/d]%mod*sm[m/i/d])%mod; 38 ans=(ans+1ll*ksm(i,i)*res)%mod; 39 } 40 printf("%d ",(ans%mod+mod)%mod); 41 return 0; 42 }