模拟赛t3 太阳神(ra) 题解

太阳神

求满足如下条件的数对$(a,b)$对数:$a,b$均为正整数且$a,b \leq n$而$lcm(a,b)>n$。

答案对$10^9+7$取模 $n\leq 10^{10}$.

原题题解写的看不懂

题意即为求

$\sum _{a=1} ^{N} \sum _{b=1} ^{N} [lcm(a,b)>N]$.

转化为

$N^2-\sum _{a=1} ^{N} \sum _{b=1} ^{N} [a*b/gcd(a,b) \leq N]$

把$a,b$最大公约数$d$提出来，化为

$N^2-\sum _{d=1}^{N}\sum _{a=1} ^{\lfloor \frac{N}{d} \rfloor} \sum _{b=1} ^{\lfloor \frac{N}{d} \rfloor} [a*b \leq \lfloor \frac{N}{d} \rfloor][gcd(a,b)==1]$

之后令

$f(x)=\sum _{a=1} ^{x} \sum _{b=1} ^{x} [a*b \leq x][gcd(a,b)==1]$

于是我们可以先枚举$\lfloor \frac{N}{d} \rfloor$，算出$\sum _{d=1} ^{N} f(\lfloor \frac{N}{d}\rfloor )$即可，由于$\lfloor \frac{N}{d} \rfloor$只有$\sqrt N$个，所以只需要利用数论分块枚举$\sqrt N$次即可！

接下来考虑如何快速求出$f(d)$

我们发现：

$f(d)=\sum _{a=1} ^{d} \sum _{b=1} ^{d} [a*b \leq d][gcd(a,b)==1]$

$=\sum _{a=1} ^{d-1} \sum _{b=1} ^{d-1} [a*b \leq d-1][gcd(a,b)==1]+\sum _{a=1} ^{d} \sum _{b=1} ^{d} [a*b==d][gcd(a,b)==1]$

$=f(d-1)+\sum _{a=1} ^{d} \sum _{b=1} ^{d} [a*b==d][gcd(a,b)==1]$

相当于把$d$质因数分解为$\prod p_i ^ {a_i}$,对于每个$p_i ^ {a_i}$,要么给$a$，要么给$b$，所以有$2^{sum[d]}$种方案，$sum[d]$即$d$的质数约数个数。

于是

$f(d)=f(d-1)+2^{sum[d]}$

$sum[d]$是可以线性筛的，于是我们就可以$O(N)$求$f(d)$了!

但是考虑到过于大的数据范围，我们考虑优化：

令$S(x)=\sum _{a=1} ^{x} \sum _{b=1} ^{x} [a*b \leq x]$

则有$f(x)=\sum _{a=1} ^{x} \sum _{b=1} ^{x} [a*b \leq x] - \sum _{t=2} ^{\sqrt x} \sum _{a=1} ^{x} \sum _{b=1} ^{x} [a*b \leq x][gcd(a,b)==t]$

$=S(x)-\sum _{t=2} ^{\sqrt x} \sum _{a=1} ^{\frac{x}{t^2}} \sum _{b=1} ^{\frac{x}{t^2}} [a*b \leq \frac{x}{t^2}][gcd(a,b)==1]$

$=S(x)-\sum _{t=2} ^{\sqrt x}f(\lfloor \frac{x}{t^2} \rfloor)$

这里面$S(x)$是可以数论分块$O(\sqrt x)$求得，于是我们可以递归来求$f(d)$了！别忘了加上记忆化，这里可以用$N/x$来记忆化。

我们在不爆空的情况下尽可能多筛一些$f(d)$，我这里就筛到了$10^7$.对于每一个$\lfloor \frac{N}{d} \rfloor$,时间复杂度就是$O(\sqrt {\frac{N}{d}})$.

计算一下总时间复杂度($log _{10^{10}} {10^7} = 0.7$)：

$n^{0.7} \leq \frac{n}{d} \Rightarrow d \leq n^{0.3}$

于是有

$\displaystyle \int _1 ^{n^{0.3}} \sqrt {\frac{n}{x}}\,dx$ 

$=\sqrt n \displaystyle \int _1 ^{n^{0.3}} x^{-\frac{1}{2}}\,dx$

$=n^{\frac{1}{2}} 2x^{\frac{1}{2}}| _{1} ^{n^{0.3}}$

$= 2n^{0.5}(n^{0.15}-1)$

$\approx n^{0.65}$

所以时间复杂度为$O(n^{0.65})$

#define ll long long
#define chr 10000000
const int maxn=1e7+5;
int f[maxn],num[maxn],p[maxn/10],cnt=0;
bool vis[maxn];
ll N;
const int mod=1e9+7;
void gt(int n){
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            p[++cnt]=i;
            num[i]=1;
        }
        for(int j=1;j<=cnt&&p[j]*i<=n;j++){
            vis[p[j]*i]=1;
            if(i%p[j]==0){
                num[p[j]*i]=num[i];
                break;
            }else num[p[j]*i]=num[i]+1;
        }
    }
    for(int i=1;i<=n;i++)
        f[i]=(f[i-1]+(1ll<<num[i]))%mod;
}
int F[maxn];
int f_(ll x){
    if(x<=chr) return f[x];
    if(~F[N/x]) return F[N/x];
    ll ans=0;
    for(ll i=1;i<=x;){
        ll j=x/(x/i);
        ans+=x/i*(1+j-i)%mod;
        ans%=mod;
        i=j+1;
    }
    for(ll i=2;i*i<=x;i++) (ans-=f_(x/(i*i)))%=mod;
    return F[N/x]=ans; 
}
int MAIN(){
    // file();
    cin>>N;gt(chr);memset(F,-1,sizeof F);
    ll ans=0;
    for(ll i=1;i<=N;"chr"){
        ll j=N/(N/i);
        (ans+=(ll)f_(N/i)*(j-i+1)%mod)%=mod;
        i=j+1;
    }
    N%=mod;
    cout<<((N*N-ans)%mod+mod)%mod<<endl;
    return 0;
}