51 Nod 1238 最小公倍数之和 V3 杜教筛

题目链接：http://www.51nod.com/Challenge/Problem.html#!#problemId=1238

题意：求$sum_{i=1}^{n}sum_{j=1}^{n}lcm(i,j)=sum_{i=1}^{n}sum_{j=1}^{n}{frac{i*j}{gcd(i,j)}}$,$1leq{n}leq10^{10}$.

知识提要：小于等于n中与n互质的数总和为$sum_{i=1}^{n}[(n,i)=1]i=frac{varphi(n)*n+[n=1]}{2}$

解析：

枚举最大公约数d，

$$Ans=sum_{d=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor}[(i,j)=1]i*j$$

我们先考虑 j<=i 的情况，

$$quadsum_{d=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{i}[(i,j)=1]i*j\$$

$$=sum_{d=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}frac{varphi(i)*i+[i=1]}{2}*i$$

还有i<=j的情况没考虑,其实两者是对称的 ,上面的式子乘2就好了,然后(1,1)这一对多算了一次了,所以-1就好了，

$$Ans=sum_{d=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}varphi(i)*i^2$$

令$F(n)=sum_{i=1}^{n}varphi(i)*i^2$ 

$$Ans=sum_{d=1}^{n}d*F(lfloorfrac{n}{d} floor)$$

欧拉函数的前缀和$phi(n)$之前博客里写过 按照类似的方法可以推出来

$$F(n)=frac{n^2*(n+1)^2}{4}-sum_{i=2}^{n}sum_{j=1}^{lfloorfrac{n}{i} floor}varphi(j)*i^2*j^2\$$
$$=frac{n^2*(n+1)^2}{4}-sum_{i=2}^{n}i^2sum_{j=1}^{lfloorfrac{n}{i} floor}varphi(j)*j^2\$$
$$=frac{n^2*(n+1)^2}{4}-sum_{i=2}^{n}i^2F(lfloorfrac{n}{i} floor)$$

到此为止可以$O(n^{frac{2}{3}})$求出Ans

AC代码

#include <bits/stdc++.h>
#define pb push_back
#define mp make_pair
#define fi first
#define se second
#define all(a) (a).begin(), (a).end()
#define fillchar(a, x) memset(a, x, sizeof(a))
#define huan printf("
");
#define debug(a,b) cout<<a<<" "<<b<<" ";
using namespace std;
const int maxn=1e6+10,inf=0x3f3f3f3f;
typedef long long ll;
const ll mod = 1000000007;
typedef pair<int,int> pii;
int check[maxn],prime[maxn],phi[maxn],sum[maxn];
void Phi(int N)//线性筛
{
    int pos=0;sum[0]=0;
    sum[1]=phi[1]=1;
    for(ll i = 2 ; i <= N ; i++)
    {
        if (!check[i])
            prime[pos++] = i,phi[i]=i-1;
        for (int j = 0 ; j < pos && i*prime[j] <= N ; j++)
        {
            check[i*prime[j]] = 1;
            if (i % prime[j] == 0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            else
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
        sum[i]=(sum[i-1]+(phi[i]*i%mod)*i%mod)%mod;
    }
}
unordered_map<ll,ll> ma;
ll inv2=500000004;
ll inv4=250000002;
ll inv6=166666668;
ll solve(ll n)
{
    if(n<=1e6)
        return sum[n];
    else if(ma.count(n))
        return ma[n];
    ll temp = ((n%mod)*((n+1)%mod)%mod)*inv2%mod;
    temp=temp*temp%mod;
    for(ll i=2,j;i<=n;i=j+1)
    {
        j=n/(n/i);
        ll r=(((j%mod)*((j+1)%mod)%mod)*((2*j+1)%mod)%mod)*inv6%mod;
        ll l=(((i%mod)*((i-1)%mod)%mod)*((2*i-1)%mod)%mod)*inv6%mod;
        r=(r-l+mod)%mod;
        temp = (temp-solve(n/i)*r%mod+mod)%mod;
    }
    return ma[n]=temp;
}
int main()
{
    ll n;
    Phi(1e6);
    scanf("%lld",&n);
    ll ans=0;
    for(ll i=1,j;i<=n;i=j+1)
    {
        j=n/(n/i);
        ll r=((j%mod)*((j+1)%mod)%mod)*inv2%mod;
        ll l=((i%mod)*((i-1)%mod)%mod)*inv2%mod;
        r=(r-l+mod)%mod;
        ans=(ans+solve(n/i)*r%mod)%mod;
    }
    printf("%lld
",ans);
}