51Nod1222:最小公倍数计数 (莫比乌斯反演)

题目传送门:http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1222


题目分析:这题又是skywalkert博客里的练习题,我一开始自己手推原式子推出了一个 O(nlog2(n)) 的东西(没有用 μ ),完全没办法化简,还是要跑去看题解,后来发现别人的思路和我一开始就不是一个方向的……也许我数论还是太弱了吧。


题目要我们求lcm(x,y)在[a,b]内的数对(x,y)的数量,其中x<=y。我们先用差分,设F(n)表示lcm(x,y)在[1,n]内的数对(x,y)的数量(x<=y),然后用F(b)-F(a-1)即可。
然后我们发现x<=y这个限制很难做,于是我们考虑G(n)=lcm(x,y)在[1,n]内的数对(x,y)的数量(不要求x<=y),则F(n)=(G(n)+n)/2。
现在变成要求:

G(n)=x=1ny=1n[xy(x,y)<=n]

接下来我们枚举d=(x,y),令 i=xd,j=yd 得到:
G(n)=d=1ni=1ndj=1nd[ijd<=n][(i,j)=1]

接下来我们套用莫比乌斯反演的公式 d|nμ(d)=[n=1] ,将最右边的那个条件限制(i,j)=1改为枚举:
G(n)=d=1ni=1ndj=1nd[ijd<=n]k|i,k|jμ(k)

然后我们把对k的枚举提到第二层(即d之后),重新令 x=ik,y=jk ,得到:
G(n)=d=1nk=1ndμ(k)x=1ndky=1ndk[dk2xy<=n]

我们可以把最外层的两层枚举看成是 先枚举d,然后再枚举d的倍数k,那么我们可以变成 先枚举倍数k,再枚举有多少个符合条件的约数d,也就是说 nd=1ndk=1μ(k) nk=1μ(k)nkd=1 是一样的。
同时,我们把 [dk2xy<=n] 变形为 [dxy<=nk2] ,得到:
G(n)=k=1nμ(k)d=1nkx=1nkdy=1nkd[dxy<=nk2]


推到这里之后,我们就可以看到这道题的巧妙之处了。关键在于最右边的 [dxy<=nk2] ,这意味着我们的k只需要枚举到 n 就可以了,再枚举下去也不会对答案产生贡献。同时我们的d只需要枚举到 nk2 ,x只需枚举到 nk2d ,y只需枚举到 nk2dx 即可。再枚举下去也不会对答案产生贡献,这样原式就变成了:

G(n)=k=1nμ(k)H(nk2)

其中:
H(n)=a=1nb=1nac=1nab1

那么我们怎么快速地求H(n)呢?我们不妨令a<=b<=c,然后将a从1枚举到 n13 ,b从a枚举到 na ,直接算出c的个数为num。根据组合数学,当 a<b<c 的时候,将6*num加进答案;当 a=b<c a<b=c 的时候,将3*num加进答案;当a=b=c的时候,将num加进答案。
那么这样的时间复杂度是多少呢?网上神犇说可以用积分证得约为 O(n23) (然而我并不会证,我们学校搞数学竞赛的人好像也不会证……),总之6组数据的总时间为33s(然而出题人用10s,我的慢好多)。


CODE:

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=400100;
typedef long long LL;

int miu[maxn];
bool vis[maxn];
int prime[maxn];
int cur=0;

LL a,b;

void Preparation()
{
    miu[1]=1;
    for (int i=2; i<maxn; i++)
    {
        if (!vis[i]) prime[++cur]=i,miu[i]=-1;
        for (int j=1; j<=cur && i*prime[j]<maxn; j++)
        {
            int k=i*prime[j];
            vis[k]=true;
            if (i%prime[j]) miu[k]=-miu[i];
            else
            {
                miu[k]=0;
                break;
            }
        }
    }
}

LL Get(LL x)
{
    LL temp=0;
    int A=(int)pow(x,1.0/3.0)+1;
    for (int a=1; a<=A; a++)
    {
        LL B=x/(long long)a;
        B=(long long)floor( sqrt( (long double)B )+0.0000001 );
        for (LL b=a+1; b<=B; b++)
        {
            LL c=x/( (long long)a*b );
            c-=b;
            if (c>0) temp+=( c*6LL );
            if ( (long long)a*b*b<=x ) temp+=3LL;
        }
    }
    for (int a=1; a<=A; a++)
    {
        LL c=x/( (long long)a*(long long)a );
        c-=a;
        if (c>0) temp+=( c*3LL );
        if ( (long long)a*(long long)a*(long long)a<=x ) temp++;
    }
    return temp;
}

LL Work(LL n)
{
    int sn=(int)floor( sqrt( (long double)n )+0.0000001 );
    LL temp=0;
    for (int i=1; i<=sn; i++)
    {
        LL x=n;
        x/=(long long)i;
        x/=(long long)i;
        if (!x) break;
        temp+=( (long long)miu[i]*Get(x) );
    }
    temp+=n;
    temp/=2LL;
    return temp;
}

int main()
{
    freopen("c.in","r",stdin);
    freopen("c.out","w",stdout);

    Preparation();
    cin>>a>>b;
    LL ans=Work(b)-Work(a-1);
    cout<<ans<<endl;

    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值