51nod 1227 平均最小公倍数杜教筛

最新推荐文章于 2020-08-20 19:53:49 发布

olahiuj

最新推荐文章于 2020-08-20 19:53:49 发布

阅读量229

点赞数

CC 4.0 BY-SA版权

分类专栏： c++ 杜教筛文章标签： 51nod

本文链接：https://blog.youkuaiyun.com/jpwang8/article/details/81122620

c++ 同时被 2 个专栏收录

1072 篇文章

订阅专栏

杜教筛

9 篇文章

订阅专栏

本文详细解析了一个涉及数论和组合数学的复杂数学问题，通过引入杜教筛算法优化了求解过程，提供了一段完整的C++代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Description

求

\sum i = a b 1 n \sum j = 1 i l c m (i, j)

$\sum_{i=a}^b\frac{1}{n}\sum_{j=1}^i lcm(i,j)$

Solution

明天回家，现在有点浑浑噩噩
终于推出了与题解一致的柿子，so moved
照例推柿子，区间可以变成前缀和之差，枚举gcd

a n s = \sum d = 1 b \sum i = 1 ⌊ b d ⌋ \sum j = 1 i j \cdot [gcd (i, j) = 1]

$ans=\sum_{d=1}^b\sum_{i=1}^{\lfloor\frac{b}{d}\rfloor}\sum_{j=1}^{i}j\cdot [\gcd(i,j)=1]$
然后后面这一坨实际上就是

∑ni=1i⋅[gcd(n,i)=1]=φ(n)⋅n2∑i=1ni⋅[gcd(n,i)=1]=φ(n)⋅n2 $\sum_{i=1}^ni\cdot[gcd(n,i)=1]=\frac{\varphi(n)\cdot n}{2}$ 的应用，带进去就得到

a n s = \sum d = 1 b \sum i = 1 ⌊ b d ⌋ φ ( i ) \cdot i 2

$ans=\sum_{d=1}^b\sum_{i=1}^{\lfloor\frac{b}{d}\rfloor}\frac{\varphi(i)\cdot i}{2}$
需要注意的是当i=1的时候后面的贡献应当是1但是柿子体现不出来，需要我们单独加上去。这个和除二一起先不管
令

S(n)=sumni=1φ(i)⋅iS(n)=sumi=1nφ(i)⋅i $S(n)=sum_{i=1}^n\varphi(i)\cdot i$ ，那么

a n s = b + 1 2 \sum d = 1 b S (⌊ b d ⌋)

$ans=b+\frac{1}{2}\sum_{d=1}^bS\left(\lfloor\frac{b}{d}\rfloor\right)$
然后问题就变成了怎么求

S(n)S(n) $S(n)$ ，这个是非常经典的杜教筛问题我就不写了

Code

#include <stdio.h>
#include <string.h>
#include <map>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)

typedef long long LL;
const int MOD=1000000007;
const int N=10000005;

std:: map <LL,LL> map;

int prime[N/10];
LL f[N+5],ny6,ny2;
bool not_prime[N+5];

void pre_work(int n) {
    f[1]=1;
    rep(i,2,n) {
        if (!not_prime[i]) {
            prime[++prime[0]]=i;
            f[i]=(i-1);
        }
        for (int j=1;i*prime[j]<=n&&j<=prime[0];j++) {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) {
                f[i*prime[j]]=f[i]*prime[j]%MOD;
                break;
            }
            f[i*prime[j]]=f[i]*(prime[j]-1)%MOD;
        }
    }// f[0]=1;
    rep(i,1,n) f[i]=(f[i-1]+f[i]*(LL)i%MOD)%MOD;
}

LL ksm(LL x,LL dep) {
    LL ret=1;
    while (dep) {
        if (dep&1) ret=ret*x%MOD;
        x=x*x%MOD; dep/=2;
    }
    return ret;
}

LL get_f(LL n) {
    if (n<=N) return f[n];
    if (map[n]) return map[n];
    LL ret=n*(n+1)%MOD*(2*n%MOD+1)%MOD*ny6%MOD;
    for (LL i=2,j;i<=n;i=j+1) {
        j=n/(n/i);
        LL tmp=(j*(j+1)%MOD-(i-1)*i%MOD+MOD)%MOD*ny2%MOD;
        ret=(ret+MOD-tmp*get_f(n/i)%MOD)%MOD;
    }
    ret=(ret%MOD+MOD)%MOD;
    return map[n]=ret;
}

LL solve(LL n) {
    if (n==0) return 0;
    LL ret=0;
    for (LL i=1,j;i<=n;i=j+1) {
        j=n/(n/i);
        ret=(ret+(j-i+1)*get_f(n/i)%MOD)%MOD;
    }
    ret=(ret+n)%MOD;
    ret=ret*ny2%MOD;
    return ret;
}

int main(void) {
    ny2=ksm(2,MOD-2);
    ny6=ksm(6,MOD-2);
    pre_work(N);
    LL a,b; scanf("%lld%lld",&a,&b);
    LL ans=solve(b);
    ans=ans-solve(a-1);
    ans=(ans%MOD+MOD)%MOD;
    printf("%lld\n", ans);
    return 0;
}