太阳神

题目描述

太阳神拉很喜欢最小公倍数,有一天他想到了一个关于最小公倍数的题目。
求满足如下条件的数对(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);
}
阅读更多
版权声明:本文是蒟蒻写出来的,神犇转载也要说一声哦! https://blog.csdn.net/WerKeyTom_FTD/article/details/53106093
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

关闭
关闭
关闭