[机房测试]太阳神

\(2021.9.13\)

Description

给定 \(n\leq 10^{10}\)，求

\[\sum_{i=1}^n \sum_{j=1}^n \Big[\frac{ij}{(i,j)}> n\Big] \]

Solution

补集转换，然后进一步化简。

\[\begin{align} &n^2-\sum_{i=1}^n \sum_{j=1}^n \Big[\frac{ij}{(i,j)} \leq n\Big] \\ =&n^2-\sum_{d=1}^n \sum_{k=1}^{\frac{n}{d}} \mu(k) \sum_{i=1}^{\frac{n}{kd}} \sum_{j=1}^{\frac{n}{kd}} \Big[ijk^2d\leq n\Big] \end{align} \]

特别的点在于很多上界都是无效的，因为往大了取对答案没有影响。然后 \(k\) 比较特殊，带了一个 \(\mu\) 的项，其他项都是类似的，只是上界不同。考虑把 \(k^2\) 除过去，那么 \(k\) 就可以枚举了，因为有 \(k^2\leq n\)。然后更新一下上界，就可以得到

\[\sum_{k=1}^{\sqrt{n}} \mu(k) \sum_{d=1}\sum_{i=1}\sum_{j=1} [ijd\leq \frac{n}{k^2}] \]

所以可以枚举 \(k\)，然后计算三元组 \((i,j,d)\) 的个数，令 \(rg=\frac{n}{k^2}\)。因为这三者完全对称，所以可以钦定一种大小关系然后乘上一个系数。例如，一种是 \(a<b<c\) ，就有 \(6\) 种可以将 \(i,j,d\) 代入的方式。考虑如何计算。因为 \(a\) 是最小的，所以有 \(a^3\leq rg\)，那么 \(a\) 就可以枚举了，再枚举一个 \(b\) 满足 \(b^2\leq \frac{rg}{a}\)，然后 \(d\) 的范围就可以快速计算。剩下还有其中两个数相同和三个都相同的情况，不再赘述。复杂度 \(O(能过)\)。

#include<stdio.h>
#include<math.h>

typedef long long ll;

const int N=1e5+7;
const int Mod=1e9+7;

bool mk[N];
int mu[N],p[N],cnt=0;

ll calc(ll rg){
    ll ret=0;
    for(ll i=1;i*i*i<=rg;i++)
        for(ll j=i+1;j+1<=(rg/(i*j));j++)
            ret=(ret+6ll*(rg/(i*j)-j)%Mod)%Mod;
    for(ll i=1;i*i<=rg;i++)
        ret=(ret+3ll*((rg/(i*i))%Mod-(i<=(rg/(i*i)))))%Mod;
    for(ll i=1;i*i*i<=rg;i++) ret++;
    return ret%Mod;
}

int main(){
    freopen("ra.in","r",stdin);
    freopen("ra.out","w",stdout);
    ll n; scanf("%lld",&n);
    int rg=(int)(sqrt(n)+0.5); mu[1]=1;
    for(int i=2;i<=rg;i++){
        if(!mk[i]) p[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&p[j]*i<=rg;j++){
            mk[p[j]*i]=1;
            if(i%p[j]==0) break;
            mu[p[j]*i]=-mu[i];
        }
    }
    ll ans=0;
    for(ll d=1;d<=rg;d++)
        ans=(ans+mu[d]*calc(n/(d*d))+Mod)%Mod;
    printf("%lld",((n%Mod)*(n%Mod)-ans+Mod)%Mod);
}