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

Description

---

求

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

Solution

---

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

$$ans=\sum_{d=1}^b\sum_{i=1}^{\lfloor\frac{b}{d}\rfloor}\sum_{j=1}^{i}j\cdot [\gcd(i,j)=1]$$
然后后面这一坨实际上就是 $\sum_{i=1}^ni\cdot[gcd(n,i)=1]=\frac{\varphi(n)\cdot n}{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)=sum_{i=1}^n\varphi(i)\cdot i$，那么 $$ans=b+\frac{1}{2}\sum_{d=1}^bS\left(\lfloor\frac{b}{d}\rfloor\right)$$
然后问题就变成了怎么求 $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;
}