JZOJ 4872 太阳神

首先取补集。
即枚举 i=agcd(a,b),j=bgcd(a,b),d=gcd(a,b)，条件即 gcd(i,j)=1,dij≤n。
固定了 i,j 后，d 的取值只有 ⌊nij⌋ 种。

然后有 n∑i=1n∑j=1[gcd(i,j)=1]⌊nij⌋=n∑i=1n∑j=1⌊nij⌋∑d|gcd(i,j)μ(d)=n∑d=1μ(d)⌊nd⌋∑i=1⌊nd⌋∑j=1⌊nd2ij⌋

修正上界，原式变为 ⌊√n⌋∑d=1μ(d)⌊nd2⌋∑i=1⌊nd2i⌋∑j=1⌊nd2ij⌋

设 f(n)=n∑i=1⌊ni⌋∑j=1⌊nij⌋。
发现后面的和式很有趣，如果能求出 g(n)=n∑i=1⌊ni⌋ 就可以数论分块了。

观察发现 g 的定义，发现 g(n)=n∑i=1σ0(i)。
于是可以类似杜教筛，先筛出 n2/3 以内的 g 值，然后在 n2/3 以外的 g 值暴力数论分块做。
可以发现这一部分很像杜教筛所以复杂度也不会超过 O(n2/3)（人家还没递归）
其实可以令 σ0∗μ 得到 1 这样就可以了两遍杜教筛带走了（但是常数过大）

然后对 f 数论分块，求出 ⌊√n⌋∑d=1μ(d)f(⌊nd2⌋) 即可。
别忘了取补集回去。

代码：

#include <cstdio>
#include <cstring>
using namespace std;
const long long N = 1e10;
const int MX = 1e7;
const int mod = 1e9 + 7;
long long n;
int vis[MX + 5],cnt,prime[MX + 5],mu[MX + 5],d[MX + 5],mn[MX + 5],ans;
int mem[MX + 5];
inline int id(long long x)
{
    return n / x;
}
int g(long long n)
{
    if(n <= MX)
        return d[n];
    if(~mem[id(n)])
        return mem[id(n)];
    int ret = 0;
    for(register long long l = 1,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ret = (ret + (long long)(r - l + 1) * (n / l) % mod) % mod;
    }
    return mem[id(n)] = ret;
}
int f(long long n)
{
    int ret = 0;
    for(register long long l = 1,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ret = (ret + (long long)(r - l + 1) * g(n / l) % mod) % mod;
    }
    return ret;
}
int main()
{
    freopen("ra.in","r",stdin),freopen("ra.out","w",stdout);
    memset(mem,-1,sizeof mem),mu[1] = d[1] = 1;
    for(register int i = 2;i <= MX;++i)
    {
        if(!vis[i])
            mu[prime[++cnt] = i] = mod - 1,d[i] = mn[i] = 2;
        for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))
            {
                d[i * prime[j]] = d[i] / mn[i] * (mn[i] + 1),mn[i * prime[j]] = mn[i] + 1;
                break;
            }
            else
                d[i * prime[j]] = d[i] * 2,mn[i * prime[j]] = 2,mu[i * prime[j]] = mod - mu[i];
        }
        d[i] = (d[i - 1] + d[i]) % mod;
    }
    scanf("%lld",&n);
    for(register int i = 1;(long long)i * i <= n;++i)
        ans = (ans + (long long)mu[i] * f(n / i / i) % mod) % mod;
    ans = ((long long)(n % mod) * (n % mod) % mod - ans + mod) % mod;
    printf("%d\n",ans);
}