BZOJ3561 DZY Loves Math VI【莫比乌斯反演】

Time Limit: 10 Sec
Memory Limit: 256 MB

Description

给定正整数n,m。求

Input

一行两个整数n,m

Output

一个整数，为答案模1000000007后的值

HINT

1<=n,m<=500000，共有3组数据

---

题目分析

和 洛谷P1829 [国家集训队]Crash的数字表格的推导过程很像

$$\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j{)}^{gcd(i,j)}=\sum _{i=1}^n\sum _{j=1}^m(\frac{ij}{gcd(i,j)}{)}^{gcd(i,j)}$$

把gcd提出来枚举

$$\sum _{d=1}^{min(n,m)}\sum _{i=1}^n\sum _{j=1}^m(\frac{ij}{d}{)}^d[gcd(i,j)=d]$$

若$gcd(i,j)==d$，则有$gcd(i/d,j/d)=1$

$$\sum _{d=1}^{min(n,m)}\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }(\frac{xd\ast yd}{d}{)}^d[gcd(x,y)=1]=\sum _{d=1}^{min(n,m)}{d}^d\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }(xy{)}^d[gcd(x,y)=1]$$

莫比乌斯反演一下

$$\sum _{d=1}^{min(n,m)}{d}^d\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }(xy{)}^d\sum _{z|gcd(x,y)}\mu (z)=\sum _{d=1}^{min(n,m)}{d}^d\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }(xy{)}^d\sum _{z=1}^{gcd(x,y)}\mu (z)[z|gcd(x,y)]$$

变换一下枚举顺序
$$\sum _{d=1}^{min(n,m)}{d}^d\sum _{z=1}^{min(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )}\mu (z)\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }(xy{)}^d[z|gcd(x,y)]$$
$x,y$的枚举替换成枚举$z$的倍数
$$\sum _{d=1}^{min(n,m)}{d}^d\sum _{z=1}^{min(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )}\mu (z)\sum _{t_1=1}^{\lfloor \frac{n}{dz}\rfloor }\sum _{t_2=1}^{\lfloor \frac{m}{dz}\rfloor }(z^2{t}_1{t}_2{)}^d=\sum _{d=1}^{min(n,m)}{d}^d\sum _{z=1}^{min(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )}{z}^{2d}\mu (z)\sum _{t_1=1}^{\lfloor \frac{n}{dz}\rfloor }{t}_1^d\sum _{t_2=1}^{\lfloor \frac{m}{dz}\rfloor }{t}_2^d$$

虽然由于每一项几乎都有d次幂不方便维护前缀和
但认真观察发现最后那两个和式每次从上一次d-1次幂继承，只用$O(m/d)$处理
$z^{2d}\mu (z)$需要$O(n/d\ast logd)$处理
总复杂度是$O(nlogn)$级的

由于每个测试点只有一组数据，完全可以接受

---

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<queue>
using namespace std;
typedef long long lt;
#define lowbit(x) ((x)&(-x))
 
int read()
{
    int f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return f*x;
}
 
const int mod=1000000007;
const int maxn=500010;
int T;
int miu[maxn];
int prim[maxn],vis[maxn],cnt;
lt sum[maxn],a[maxn];
 
lt qpow(lt a,int k)
{
    lt res=1;
    while(k>0){
        if(k&1) res=(res*a)%mod;
        a=(a*a)%mod; k>>=1;
    }
    return res;
}
 
void Miu(int n)
{
    miu[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!vis[i]) prim[++cnt]=i,miu[i]=-1;
        for(int j=1;j<=cnt;++j)
        {
            if(i*prim[j]>n) break;
            vis[i*prim[j]]=1;
            if(i%prim[j]==0) break;
            else miu[i*prim[j]]=-miu[i];
        }
    }
}
 
lt query(int n,int m)
{
    lt res=0; if(n>m) swap(n,m);
    for(int i=1;i<=m;++i) a[i]=1;
     
    for(int d=1;d<=n;++d)
    {
        for(lt i=1;i<=m/d;++i)
        {
            a[i]=a[i]*i%mod;
            sum[i]=(sum[i-1]+a[i])%mod;
        }
         
        lt tt=0;
        for(int j=1;j<=n/d;++j)
        {
            tt+=(qpow(j,2*d)*miu[j]%mod+mod)%mod * sum[n/(j*d)]%mod * sum[m/(j*d)]%mod;
            tt%=mod;
        }
         
        res+=qpow(d,d)*tt%mod; 
        res%=mod;
    }
    return res;
}
 
int main() 
{
    int n=read(),m=read();
    Miu(min(n,m));
    printf("%lld",query(n,m));
    return 0;
}