BZOJ 4176 Lucas的数论 莫比乌斯反演

题目大意：给定$n(n\leq10^9)$，求$\sum_{i=1}^n\sum_{j=1}^nd(ij)$

推错式子害死人。。。
由$d|ij$等价于$\frac d{\gcd(i,d)}|j$可得
$\sum_{i=1}^n\sum_{j=1}^nd(ij)$
$=\sum_{i=1}^n\sum_{d=1}^{n^2}\lfloor\frac{n*\gcd(i,d)}d\rfloor$
$=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac{n^2}d\rfloor}\lfloor\frac nj\rfloor[\gcd(i,j)=1]$
$=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^n\lfloor\frac nj\rfloor[\gcd(i,j)=1]$
$=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^n\lfloor\frac nj\rfloor\sum_{k|i,k|j}\mu(k)$
$=\sum_{k=1}^n\mu(k)(\sum_{d=1}^{\lfloor\frac nk\rfloor}\lfloor\frac n{kd}\rfloor)^2$
$O(\sqrt n)$枚举$\lfloor\frac nk\rfloor$
$\mu(k)$的部分同BZOJ3944
后面那个平方里面的东西每次$O(\sqrt{\frac nk})$求 时间复杂度$O(\sqrt\frac n1)+O(\sqrt\frac n2)+...+O(\sqrt\frac n{\sqrt n})=O(n^{\frac34})$
总时间复杂度$O(n^{\frac 34})$

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 5700000
#define MOD 1000000007
using namespace std;

int n,S;

int Get_Sum1(int n)
{
    int i,last,re=0;
    for(i=1;i<=n;i=last+1)
    {
        last=n/(n/i);
        (re+=(long long)(last-i+1)*(n/i)%MOD)%=MOD;
    }
    return re;
}

int _sum2[M];
int mu[M],prime[M],tot;
bool not_prime[M];
void Linear_Shaker()
{
    int i,j;
    mu[1]=1;
    for(i=2;i<=S;i++)
    {
        if(!not_prime[i])
        {
            mu[i]=-1;
            prime[++tot]=i;
        }
        for(j=1;prime[j]*i<=S;j++)
        {
            not_prime[prime[j]*i]=true;
            if(i%prime[j]==0)
            {
                mu[prime[j]*i]=0;
                break;
            }
            mu[prime[j]*i]=-mu[i];
        }
    }
}
int Get_Sum2(int x)
{
    if(x<=S) return mu[x];
    return _sum2[n/x];
}
void Pretreatment()
{
    int i,j,last;
    Linear_Shaker();
    for(i=1;i<=S;i++)
        mu[i]+=mu[i-1];
    for(i=1;n/i>S;i++);
    for(j=i;j;j--)
    {
        int n=::n/j;
        _sum2[j]=1;
        for(i=2;i<=n;i=last+1)
        {
            last=n/(n/i);
            (_sum2[j]-=(long long)(last-i+1)*Get_Sum2(n/i)%MOD)%=MOD;
        }
    }
}

int main()
{
    int i,last,ans=0;
    cin>>n;S=ceil(pow(n,0.75)-1e-7)+1e-7;
    Pretreatment();
    for(i=1;i<=n;i=last+1)
    {
        last=n/(n/i);
        long long temp=Get_Sum1(n/i);
        (ans+=(Get_Sum2(last)-Get_Sum2(i-1))*temp%MOD*temp%MOD)%=MOD;
    }
    cout<<(ans+MOD)%MOD<<endl;
    return 0;
}