[BZOJ3561]DZY Loves Math VI(莫比乌斯反演)

$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)$

#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std;

const int N=500010,mod=1e9+7;
int n,m,ans,tot,sm[N],pw[N],mu[N],p[N],b[N];

int ksm(int a,int b){
    int res=1;
    for (; b; a=1ll*a*a%mod,b>>=1)
        if (b & 1) res=1ll*res*a%mod;
    return res;
}

void init(int n){
    mu[1]=1;
    rep(i,2,n){
        if (!b[i]) p[++tot]=i,mu[i]=-1;
        for (int j=1; j<=tot && p[j]*i<=n; j++){
            b[p[j]*i]=1;
            if (i%p[j]==0) { mu[p[j]*i]=0; break; }
                else mu[p[j]*i]=-mu[i];
        }
    }
}

int main(){
    freopen("bzoj3561.in","r",stdin);
    freopen("bzoj3561.out","w",stdout);
    scanf("%d%d",&n,&m);
    if (n>m) swap(n,m); init(n);
    rep(i,1,m) pw[i]=1;
    rep(i,1,n){
        int res=0; sm[0]=0;
        rep(d,1,m/i) pw[d]=1ll*pw[d]*d%mod,sm[d]=(sm[d-1]+pw[d])%mod;
        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;
        ans=(ans+1ll*ksm(i,i)*res)%mod;
    }
    printf("%d
",(ans%mod+mod)%mod);
    return 0;
}