【51nod1238】 最小公倍数之和 V3

Description

出个数N，输出小于等于N的所有数，两两之间的最小公倍数之和。

相当于计算这段程序（程序中的lcm(i,j)表示i与j的最小公倍数）：
由于结果很大，输出Mod 1000000007的结果。

G=0;
for(i=1;i< N;i++)
for(j=1;j<=N;j++)
{
G = (G + lcm(i,j)) % 1000000007;
}

Solution

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{n}{i*j\over gcd(i,j)}=\sum_{d=1}^{n}d*\sum_{i=1}^{n/i}i*\sum_{j=1}^{n/i}j[gcd(i,j)==1]=2*(\sum_{d=1}^{n}d(1+\sum_{i=1}^{n/i}i*i*\phi(i)/2)))-(1+n)*n/2=\sum_{d=1}^{n}d\sum_{i=1}^{n/i}i*i*\phi(i)$$
然后这就是杜教筛的标准形式了，直接筛就好了。

Code

#include<iostream>
#include<math.h>
#include<string.h>
#include<stdio.h>
#include<algorithm>
#define ll long long
using namespace std;
const ll maxn=1e6,mo=1e9+7,mo4=250000002,mo2=mo4*2,mo6=166666668;
ll d[maxn],bz[maxn+5],p[maxn+5],h[maxn],f[maxn];
ll n,i,t,j,k,l,x,y,z,ans;
ll hash(ll x){
    ll t=x%maxn;
    while (h[t] && h[t]!=x) t=(t+1)%maxn;
    return t;
}
ll sqr(ll n){
    n%=mo;return n*n%mo;
}
ll make(ll n){
    return n*(n+1)%mo*(2*n+1)%mo*mo6%mo;
}
ll dg(ll n){
    if (n<=maxn) return p[n];
    ll x=hash(n),y,t,k=sqr(1+n)*sqr(n)%mo*mo4%mo,i=2;
    if (h[x]) return f[x];
    while (i<=n){
        t=n/(n/i);
        k-=(make(t%mo)-make((i-1)%mo)+mo)%mo*dg(n/i)%mo;i=t+1;
    }
    h[x]=n;f[x]=(k%mo+mo)%mo;return f[x];
}
int main(){
    freopen("data.in","r",stdin);freopen("data.out","w",stdout);
    scanf("%lld",&n);p[1]=1;
    for (i=2;i<=maxn;i++){
        if (!bz[i]) d[++d[0]]=i,p[i]=i-1;
        for (j=1;j<=d[0];j++){
            if (i*d[j]>maxn) break;
            bz[i*d[j]]=1;
            if (i%d[j]==0) {
                p[i*d[j]]=p[i]*d[j];break;
            }else p[i*d[j]]=p[i]*p[d[j]];
        }
    }
    for (i=2;i<=maxn;i++)
        p[i]=(p[i]*i%mo*i%mo+p[i-1])%mo;
    i=1;
    while (i<=n){
        t=n/(n/i);
        k=dg(t)-dg(i-1);x=n/i%mo;
        l=(x+1)*x%mo*mo2%mo;
        ans+=k*l%mo;i=t+1;
    }
    ans=(ans%mo+mo)%mo;
    printf("%lld\n",ans);
}