[20191003机房测试] 太阳神

Posted 2022-12-17 tqr06

太阳神拉很喜欢最小公倍数，有一天他想到了一个关于最小公倍数的题目
求满足如下条件的数对(a,b)对数：
a,b 均为正整数且 a,b<=n 而lcm(a,b)>n
其中的 lcm 当然表示最小公倍数
答案对 1,000,000,007取模

数据范围是1e10的，打表找了半天规律发现没用……
那就莫比乌斯反演呗

题目求：

\[\sum_i=1^n\sum_j=1^n[lcm(i,j)> n]\]

但是大于的太多了，那我们就反过来求，最后用总的来减

也就是求：

\[\sum_i=1^n\sum_j=1^n[lcm(i,j)\leq n]\]

先把\(lcm\)拆开

\[\sum_i=1^n\sum_j=1^n[\dfraci·jgcd(i,j)\leq n]\]

看到\(gcd\)，枚举 \(i,j\) 的约数 \(d\)

\[\sum_d=1^n\sum_i=1^n\sum_j=1^n[\dfraci·jd\leq n]\]

拉出去反演：

\[\sum_d=1^n\sum_i'=1^\lfloor\fracnd\rfloor\sum_j'=1^\lfloor\fracni'·d\rfloor[gcd(i',j')=1]\]

引入莫比乌斯函数：

\[\sum_d=1^n\sum_i'=1^\lfloor\fracnd\rfloor\sum_j'=1^\lfloor\fracni'·d\rfloor\mu(gcd(i',j'))\]

再枚举 \(i',j'\) 的约数 \(d'\)，继续反演

\[\sum_d=1^n\sum_d'=1^\sqrt\fracnd\mu(d')\sum_i''=1\sum_j''=1\lfloor\fracnd·d'^2i''\rfloor\]

把\(d'\)拿出来：

\[\sum_d'=1^\sqrtn\mu(d')[d·i''·j''\leq \fracnd'^2\text的组数]\]

然后就是求满足\(d·i''·j''\leq \fracnd'^2\)的三元组(d,i‘‘,j‘‘)的个数

这个可以\(\Theta(n^\frac23)\)算出来

就可以了

代码：

#include<bits/stdc++.h>
#define ll long long
#define N 100005
#define M 100000
#define mod 1000000007
using namespace std;

ll n,ans,siz,sum,maxn,x;
ll mu[N],prime[N],primenum;
bool isprime[N];

template<class T>inline void read(T &res)

    char c;T flag=1;
    while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;res=c-'0';
    while((c=getchar())>='0'&&c<='9')res=res*10+c-'0';res*=flag;


void init()

    mu[1]=1;
    for(register int i=2;i<=M;++i)
    
        if(!isprime[i])
        
            prime[++primenum]=i;
            mu[i]=-1;
        
        for(register int j=1;j<=primenum;++j)
        
            if(i*prime[j]>M) break;
            isprime[i*prime[j]]=1;
            if(!(i%prime[j]))
            
                mu[i*prime[j]]=0;
                break;
            
            mu[i*prime[j]]=-mu[i];
        
    


int main()

    freopen("ra.in","r",stdin);
    freopen("ra.out","w",stdout);
    read(n);
    init();
    siz=(ll)(sqrt(n)+0.5);
    for(register ll i=1;i<=siz;++i)
    
        sum=0;
        maxn=n/i/i;
        for(register ll j=1;j*j*j<=maxn;++j)
         
            for(register ll k=j;k*k<=maxn/j;++k)
            
                x=(maxn/j/k-k+1);
                if(j==k) sum=(sum+1+(x-1)*3)%mod;
                else sum=(sum+3+(x-1)*6)%mod;
            
        
        sum%=mod;
        ans=(ans+mu[i]*sum)%mod;
    
    ans=((n%mod)*(n%mod)-ans)%mod;
    while(ans<0) ans+=mod;
    printf("%lld\n",ans);
    return 0;

/*
10000000000
210705255
*/