51nod 1190 最小公倍数之和 V2(莫比乌斯反演)

转自:http://blog.csdn.net/u014610830/article/details/49493279

没有看懂,膜拜一下

题意: 

给定1<=a<=b<=109,求∑i=ablcm(i,b)?

 

分析: 

∑i=ablcm(i,b)=∑i=abi∗bgcd(i,b)=b∗∑d∑a<=i<=bid∗(gcd(i,b)=d)=b∗∑d|b∑⌈ad⌉<=i<=bdi∗(gcd(i,bd)=1)=b∗∑d|b∑⌈ad⌉<=i<=bdi∗∑d′|gcd(i,bd)μ(d′)=b∗∑d|b∑d′|bdμ(d′)∑⌈ad⌉<=i<=bd,d′|ii=b∗∑d|b∑d′|bdμ(d′)∗d′∗(⌈⌈ad⌉d′⌉+⌊⌊bd⌋d′⌋)∗(⌊⌊bd⌋d′⌋−⌈⌈ad⌉d′⌉+1)2=b∗∑d|b∑d′|bdμ(d′)∗d′∗(⌈add′⌉+⌊bdd′⌋)∗(⌊bdd′⌋−⌈add′⌉+1)2=b∗∑d|b(⌈ad⌉+⌊bd⌋)∗(⌊bd⌋−⌈ad⌉+1)2∑d′|dμ(d′)∗d′

 

最后一步转化是另d=dd′。

做到最后发现需要枚举b的所有约数,此处我用dfs,上一次写的状压现在t掉了,因为枚举了好多重复的约数。b的约数做多不超过200个。还有一个问题,如何求f(d)=∑d′|dμ(d′)∗d′。

易知:f(d)是一个积性函数,f(p)=1−p,f(pk)=f(p)此处p为素数。 
所以枚举约数d是,暴力判断b的所有不同的质因子是不是能够整除d,如果能的话,res=res∗f(p)=res∗(1−p)。

算一下复杂度:n分解质因数复杂度不超过o(n√),dfs找n所有的约数,最大复杂度不超过o(256),枚举约数,暴力判断是不是能够整除,o(256∗8),n的不同约数个数不超过9个。 
所以总的最大复杂度为o(T∗(n√+256∗8)),注意所有数据都达到不了这个复杂度。

下面是代码:

 

#include <bits/stdc++.h>
#define LL long long
#define FOR(i,x,y)  for(int i = x;i < y;++ i)
#define IFOR(i,x,y) for(int i = x;i > y;-- i)

using namespace std;

const LL Mod = 1000000007;
const int maxn = 100010;
const LL inv = 500000004;

LL a,b;
int prime[maxn];
bool check[maxn];
int pri[35],pri_cnt,len[35];
vector <int> fac;

void Mobius(){
    memset(check,false,sizeof(check));
    prime[0] = 0;
    FOR(i,2,maxn){
        if(!check[i]){
            prime[++prime[0]] = i;
        }
        FOR(j,1,prime[0]+1){
            if(i*prime[j] >= maxn)  break;
            check[i*prime[j]] = true;
            if(i%prime[j] == 0) break;
        }
    }
}

void dfs(int res,int l){
    if(l >= pri_cnt)    {fac.push_back(res);return;}
    int tem = 1;
    dfs(res,l+1);
    FOR(i,1,len[l]+1){
        tem *= pri[l];
        dfs(res*tem,l+1);
    }
}

void Get_Fac(){
    pri_cnt = 0;
    int b_c = b;
    FOR(i,1,prime[0]+1){
        if(prime[i]*prime[i] > b_c) break;
        if(b_c % prime[i] == 0) {pri[pri_cnt++] = prime[i]; len[pri_cnt-1] = 0;}
        while(b_c % prime[i] == 0)  {b_c /= prime[i];len[pri_cnt-1] ++;}
    }
    if(b_c > 1) {pri[pri_cnt++] = b_c; len[pri_cnt-1] = 1;}
    fac.clear();
    dfs(1,0);
}

void work(){
    Get_Fac();
    LL ans = 0;
    FOR(i,0,(int)fac.size()){
        int v = fac[i];
        LL tt = a+v-1;
        LL tem = (((tt/v+b/v)%Mod)*((b/v-tt/v+1+Mod)%Mod)%Mod)*inv%Mod;
        LL res = 1;
        FOR(j,0,pri_cnt){
            if(v % pri[j])  continue;
            LL t = (1-pri[j]+Mod)%Mod;
            res = (res*t)%Mod;
        }
        ans += (res*tem)%Mod;
        ans %= Mod;
    }
    ans = (ans*b)%Mod;
    printf("%lld\n",ans);
}

int main()
{
    //freopen("test.in","r",stdin);
    //freopen("out.txt","w",stdout);
    Mobius();
    int T;  scanf("%d",&T);
    while(T--){
        scanf("%lld%lld",&a,&b);
        work();
    }
    return 0;
}

 

 

 

 

 

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值