【51nod 1227】 平均最小公倍数(杜教筛)

题目来源： Project Euler
基准时间限制：1.5 秒 空间限制：131072 KB 分值: 640 难度：8级算法题 重点内容
Lcm(a,b)表示a和b的最小公倍数，A(n)表示Lcm(n,i)的平均数(1 <= i <= n)，
例如：A(4) = (Lcm(1,4) + Lcm(2,4) + Lcm(3,4) + Lcm(4,4)) / 4 = (4 + 4 + 12 + 4) / 4 = 6。
F(a, b) = A(a) + A(a + 1) + …… A(b)。(F(a,b) = ∑A(k), a <= k <= b)
例如：F(2, 4) = A(2) + A(3) + A(4) = 2 + 4 + 6 = 12。
给出a,b，计算F(a, b)，由于结果可能很大，输出F(a, b) % 1000000007的结果即可。
Input
输入2个数a,b，中间用空格分隔(1 <= a <= b <= 10^9)。
Output
输出F(a, b) % 1000000007的结果。
Input示例
1 100
Output示例
122726

推吧。我们先看$A(n)$.

$$A(n)=\frac{1}{n}\sum_{i=1}^nlcm(i,n)=\sum_{i=1}\frac{i}{gcd(i,n)}$$
我们枚举gcd: $$A(n)=\sum_{d|n}\sum_{d|i}[gcd(i,n)==d]\frac{i}{d}\\ =\sum_{d|n}\sum_{d|i}[gcd(\frac{i}{d},\frac{n}{d})==1]\frac{i}{d}\\ =\sum_{d|n}\sum_{i}^{\frac{n}{d}}[gcd(i,\frac{n}{d})==1]i$$
这里可以引入欧拉函数了： $$A(n)=\frac{1}{2}+\sum_{d|n}\frac{d\phi(d)}{2}$$

设：

$$F(n)=\frac{n}{2}+\frac{1}{2}\sum_{i=1}^n\sum_{j|i}j\phi(j)$$
设： $$F'(n)=\sum_{i=1}^n\sum_{j|i}j\phi(j)\\=\sum_{j=1}^nj\phi(j)\lfloor\frac{n}{j}\rfloor$$

这里就很显然了，用杜教筛搞$j\phi(j)$的前缀和，然后分块求答案。
设$f(x)=x\phi(x),g(x)=x$,这两个函数搞一下卷积。

$$f*g(n)=\sum_{d|n}d\phi(d)\frac{n}{d}=n^2$$
设 $$sum(n)=\sum_{i=1}^ni\phi(i)$$

推得：

$$sum(n)=\frac{(n+1)n(2n+1)}{6}-\sum_{i=2}^ni*sum(\lfloor\frac{n}{i}\rfloor)$$

代码：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define maxx 5000000
#define mod 1000000007
#define ll long long
using namespace std;
bool isP[maxx];
int prime[400000];
int cnt;
ll phi[maxx];
ll sum[maxx];
ll inv2=500000004,inv6=166666668;
void init()
{
    phi[1]=1;
    for(int i=2;i<maxx;i++)
    {
        if(!isP[i])
        {
            prime[cnt++]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<cnt&&(ll)i*prime[j]<maxx;j++)
        {
            isP[i*prime[j]]=true;
            if(i%prime[j])
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            else
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
        }
    }
    for(int i=1;i<maxx;i++)
        sum[i]=i*phi[i]%mod;
    for(int i=1;i<maxx;i++)
        sum[i]=(sum[i]+sum[i-1])%mod;
}
map<ll,ll>M;
ll get(ll a,ll b)
{
    return (b-a+1)*(b+a)%mod*inv2%mod;
}
ll work(ll x)
{
    if(x<maxx)return sum[x];
    if(M[x])return M[x];
    ll t=x%mod;
    ll ans=t*(t+1)%mod*(2*t%mod+1)%mod*inv6%mod;
    for(ll i=2,last;i<=x;i=last+1)
    {
        last=x/(x/i);
        ans=(ans-get(i,last)*work(x/i)%mod)%mod;
    }
    if(ans<0)ans=(ans+mod)%mod;
    M[x]=ans;
    return ans;
}
ll F(ll x)
{
    ll ans=0;
    for(ll i=1,last;i<=x;i=last+1)
    {
       last=x/(x/i);
       ll now=(work(last)-work(i-1)+mod)%mod;
       ll base=x/i%mod;
       ans=(ans+base*now%mod)%mod;
    }
    ans=ans*inv2%mod;
    ans=(ans+x*inv2%mod)%mod;
    return ans;
}
int main()
{
    init();
    ll a,b;
    cin>>a>>b;
    cout<<(F(b)-F(a-1)+mod)%mod<<endl;
    return 0;
}