BZOJ 3561 DZY Loves Math VI

既然是5e5那直接mlogm暴力咯。。。

推式子要认真。。。。233

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define maxn 500500
#define mod 1000000007
using namespace std;
long long n,m,a[maxn],sum[maxn],ans=0,prime[maxn],tot=0,miu[maxn];
bool vis[maxn];
void get_table()
{
    miu[1]=1;
    for (long long i=2;i<=maxn-500;i++)
    {
        if (!vis[i]) {prime[++tot]=i;miu[i]=-1;}
        for (long long j=1;i*prime[j]<=maxn-500 && j<=tot;j++)
        {
            vis[i*prime[j]]=true;
            if (i%prime[j]) miu[i*prime[j]]=-miu[i];
            else {miu[i*prime[j]]=0;break;}
        }
    }
}
long long f_pow(long long x,long long y)
{
    long long ans=1,base=x;
    while (y)
    {
        if (y&1) ans=(ans*base)%mod;
        base=(base*base)%mod;
        y>>=1;
    }
    return ans;
}
int main()
{
    scanf("%lld%lld",&n,&m);get_table();
    if (n>m) swap(n,m);
    for (long long i=1;i<=m;i++) a[i]=1,sum[i]=i;
    for (long long i=1;i<=n;i++)
    {
        for (long long j=1;j*i<=m;j++) {a[j]=a[j]*j%mod;sum[j]=(sum[j-1]+a[j])%mod;}
        long long ret=0;
        for (long long j=1;j*i<=n;j++) {ret+=miu[j]*sum[n/i/j]*sum[m/i/j]%mod*a[j]%mod*a[j]%mod;ret%=mod;}
        ans+=ret*f_pow(i,i);ans%=mod;
    }
    printf("%lld\n",ans);
    return 0;
}