BZOJ3561 DZY Loves Math VI

BZOJ3561 DZY Loves Math VI

题目在这里呀~

这是我做的第一道莫比乌斯反演题啊！ 其实这题算是很基础的吧，只要学过莫比乌斯反演估计就会做（显然我一开始是不会做滴）。
这其实也算个比较难的数论问题了啦.

推的过程我就发照片吧（字写得不好w）
用片段和维护sum i^p * sum j^p
剩下的时间复杂度和普通欧拉筛的时间复杂度一样。

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define N 500000+5000
#define ll long long
using namespace std;
const int mod=1e9+7;
int n,m,tot,prime[N],vis[N],mu[N];
ll s[N],sum[N],ans,res;
inline ll power(ll a,int b){
    ll mul=1;
    while(b){
        if(b & 1) mul=mul*a % mod;
        a=a*a % mod;
        b>>=1;
    }
    return mul;
}
int main()
{
    scanf("%d%d",&n,&m);
    if(n>m) swap(n,m);
    mu[1]=1;
    for(int i=2;i<=m;i++){
        if(!vis[i]){
            prime[++tot]=i;
            mu[i]=mod-1;
        }
        for(int j=1;j<=tot && i*prime[j]<=m;j++){
            vis[i*prime[j]]=1;
            if(i % prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }else mu[i*prime[j]]=(-mu[i]+mod) % mod;
        }
    }

    for(int i=1;i<=m;i++) s[i]=1;
    for(int p=1;p<=n;p++){
        for(int j=1;j<=m / p;j++){
            s[j]=(ll)s[j]*j % mod;
            sum[j]=(sum[j-1]+s[j]) % mod;
        }    //维护前缀和
        res=0;
        for(int d=1;d<=n/p;d++)
            res=(res+(ll)mu[d] * s[d] % mod * s[d] % mod * sum[n/d/p] % mod * sum[m/d/p] % mod) % mod;
        ans=(ans+res*power(p,p) % mod) % mod;  //求值
    }
    printf("%lld\n",ans);
    return 0;
}

ps：哈哈这样就看不出我字写得差啦！