HZOJ 太阳神

所以我刚学反演还没学反演就要做这么一道神仙题……

首先大于n不好求，补集转化。

$ans=n*n-sum limits _{i=1}^{n} sum limits _{j=1}^{n} left [  lcm(i,j)leqslant n ight ] $

那么我们要求：

$sum limits _{i=1}^{n} sum limits _{j=1}^{n} left [  frac{i*j}{gcd(i,j) } leqslant n ight ]$

枚举d=gcd(i，j),

原式=$sum limits _{d=1}^{n} sum limits _{i=1}^{n} sum limits _{j=1}^{n} left [ i*j*dleqslant n ,gcd(i,j)=1 ight ]$

=$sum limits _{d=1}^{n} sum limits _{i=1}^{left lfloor frac{n}{d} ight floor} sum limits _{j=1}^{left lfloor frac{n}{d*i} ight floor} left [ gcd(i,j)=1 ight ]$

根据莫比乌斯函数的性质：$sum limits _{dmid n}u(d) =left [ n=1 ight ]$

于是原式=$sum limits _{d=1}^{n} sum limits _{i=1}^{n} sum limits _{j=1}^{n} sum limits _{gmid gcd(i,j)} u(g)$

所以就要反演了？其实就是交换求和的顺序。

个人这步稍难理解（因为我没学过反演），将g提前后相当于求u（g）出现的次数，那么修改g的定义，令${i}'=frac{i}{g},{j}'=frac{j}{g}$.

原式=$sum limits _{d=1}^{n} sum limits _{g=1} u(g) sum limits _{{i}'=1}^{left lfloor frac{n}{d*g} ight floor} sum limits _{{j}'=1}^{left lfloor frac{n}{d*i*g} ight floor} 1$

将g提前，原式=$sum limits _{g=1}^{sqrt{n}}u(g) sum limits _{{i}'=1} sum limits _{{j}'=1} sum limits _{d=1} left [ {i}'*{j}'*dleqslant frac{n}{g*g} ight ]$

到此式子就推完了，可是看起来还是不是很可做……但是可以发现g是根号n范围内的，u线筛即可，同时枚举g。

不妨设${i}'leqslant {j}'leq d$,那么设$m=frac{n}{g*g}$可以以$Oleft ( m^{frac{1}{3}} ight )$的复杂度枚举i，以$sqrt{frac{m}{i}}$的复杂度枚举j，O1算出d的个数，之后乘$A_3^3$。

但是要考虑算重的情况，手动讨论一下就行了。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define int LL
#define LL long long
using namespace std;
const int mod=1e9+7;
LL n;
bool isprime[100010];
int prime[100010],cnt,mu[100010];
void shai(int n)
{    
    mu[1]=1;
    for(int i=2;i<=n;i++)isprime[i]=1;
    for(int i=2;i<=n;i++)
    {
        if(isprime[i])mu[i]=-1,prime[++cnt]=i;
        for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
        {
            isprime[prime[j]*i]=0;
            if(i%prime[j]==0)break;
            else mu[i*prime[j]]=-mu[i];
        }
    }
}
signed main()
{
    cin>>n;int maxn=sqrt(n);shai(maxn+1);
    LL ans=0;
    for(int i=1;i<=maxn;i++)
    {
        LL res=0;
        int m=n/(i*i);
        for(int a=1;a*a*a<=m;a++)
        {
            int maxb=sqrt((1.0*m)/a);
            for(int b=a;b<=maxb;b++)
            if(m/(a*b)>=b)
            {
                if(a==b)res=(res+(m/(a*b)-b)*3+1)%mod;
                else    res=(res+(m/(a*b)-b)*6+3)%mod;
            }
        }
        ans=(ans+mu[i]*res%mod)%mod;
    }
    printf("%lld
",(n%mod*(n%mod)%mod-ans%mod+mod)%mod);
}