太阳神

题目描述

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

推式子

我们先做个补集转化。
然后就是统计这样的三元组(a,b,d)个数了
1、a*b*d<=n
2、(a,b)=1
a和b定了,d的取值就有 nab
ans=na=1nb=1[(a,b)=1]nab
我们尝试莫比乌斯反演一下,也就是
f[k]=na=1nb=1[(a,b)=k]nab
g[k]=nki=1f[ki]
反演一下变成
f[k]=ni=1μ(i)g[ki]
g是很容易算的,所以
f[k]=ni=1μ(i)nkia=1nkib=1(nk2)ab
现在我们需要f[1],代入k=1得
ans=nd=1μ(d)nd2a=1nd2b=1(nd2)ab
观察后面,可以改变d的枚举上界
ans=nd=1μ(d)nd2a=1nd2b=1(nd2)ab
后面那个东西只与 nd2 有关,设其为 f(nd2)

算f

f(x)=xi=1xj=1xij
可以看出可以根据 xi 的取值分块!每一块答案显然相同,块数为 x
假如i确定了,那么答案为 xij=1(xi)j
后面那个东西只与 xi 有关,设其为 g(xi)

算g

g(x)=xi=1xi
我们应该怎么算g呢?
有这样一个奇妙的方法。
首先,先知道一个结论。
g(x)=xi=1d(i)
d函数表示约数个数。
为什么呢?你考虑 xi ,那么说明 1i,2i,3ixii 都在x内,而且i是上述提到的 xi 个数的约数。
因此可以表示为约数个数前缀和。
那么也就是说g可以线筛?
但显然不能全都线筛啊?
dalao说,以 n2/3 为界,记为N。我们需要筛出g(1~N),然后对于x<=N便直接得解。
如果x>N?观察g的算式,我们显然可以按照取值分块,得到一个根号算法。那我们就这样做最快咯。不过这样会T,加个记忆化即可。
怎么记忆化呢?dalao说,x可以映射到n/x,这样一定不会重复。这样就可以开小记忆化数组。
总复杂度?dalao说,是 O(n2/3)

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const int N=3000000,mo=1000000007;
int mu[N+10],sum[N+10],mem[N+10],d[N+10],pri[N+10];
bool bz[N+10];
int i,j,k,l,t,m,top,ans;
ll n;
void prepare(){
    int i,j;
    mu[1]=1;
    fo(i,2,N){
        if (!bz[i]){
            pri[++top]=i;
            mu[i]=-1;
        }
        fo(j,1,top){
            if ((ll)i*pri[j]>N) break;
            bz[i*pri[j]]=1;
            if (i%pri[j]==0){
                mu[i*pri[j]]=0;
                break;
            }
            mu[i*pri[j]]=-mu[i];
        }
    }
    fo(i,1,N)
        fo(j,1,N/i)
            d[i*j]++;
    fo(i,1,N) sum[i]=sum[i-1]+d[i],mem[i]=-1;
}
int calcg(ll x){
    if (x<=N) return sum[x];
    if (mem[n/x]>=0) return mem[n/x];
    ll i=1,j;
    int t=0;
    while (i<=x){
        j=x/(x/i);
        (t+=(ll)(j-i+1)*(x/i)%mo)%=mo;
        i=j+1;
    }
    return mem[n/x]=t;
}
int calcf(ll x){
    ll i=1,j;
    int t=0;
    while (i<=x){
        j=x/(x/i);
        (t+=(ll)calcg(x/i)*(j-i+1)%mo)%=mo;
        i=j+1;
    }
    return t;
}
int main(){
    freopen("ra.in","r",stdin);freopen("ra.out","w",stdout);
    prepare();
    scanf("%lld",&n);
    fo(i,1,floor(sqrt(n)))
        if (mu[i]) (ans+=mu[i]*calcf(n/i/i))%=mo;
    ans=((ll)(n%mo)*(n%mo)%mo-ans)%mo;
    (ans+=mo)%=mo;
    printf("%d\n",ans);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值