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

最新推荐文章于 2020-08-17 16:35:44 发布

原创最新推荐文章于 2020-08-17 16:35:44 发布 · 241 阅读

0 ·

CC 4.0 BY-SA版权

数论专栏收录该内容

48 篇文章

订阅专栏

本文针对求∑i=1n∑j=1nlcm(i,j)的问题，通过数学推导给出了详细的解答过程，并运用杜教筛算法进行高效求解。

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

题目链接：哆啦A梦传送门

题意：

求 $\begin{aligned}\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i,j) \end{aligned}$

我们设:
$\begin{aligned}ans&=\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i,j)\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[(i,j)=1]i*j \end{aligned}$
这步公式推导见：详情

参考博客

我们知道：

$\begin{aligned}\sum_{i=1}^{n}[(n,i)=1]*i=\frac{n*\varphi(n)+[n==1]}{2}\end{aligned}$

那么上式就可化为：
$\begin{aligned}&\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}i\sum_{j=1}^{\frac{n}{d}}[(i,j)=1]j\\ =&\sum_{d=1}^{n}d((2\sum_{i}^{\frac{n}{d}}i\sum_{j=1}^{i}[(j,i)=1]*j)-1)......解释1\\ =&\sum_{d=1}^{n}d((2\sum_{i}^{\frac{n}{d}}i\frac{i*\varphi(i)+[i==1]}{2})-1)\\ =&\sum_{d=1}^{n}d\sum_{i}^{\frac{n}{d}}i^2\varphi(i)\end{aligned}$

解释1：因为两倍的话，我们多算了(1,1)，故要减去1。两倍的时候，其实(x,x)是算了两次，但因为它们不互质，其实最后也没算进去的。

此时我们只需分块就能算出上式的结果，但n太大，没法预处理 $i^2\phi(i)$ 的前n项和。
故我们知道杜教筛就是求这个的了，我们只需把 $i^2\phi(i)$ 的前n项和给筛出来，题目就迎刃而解了。

我们设
$h(x)=x^2\varphi(x)$ ， $S(n)=\sum_{i=1}^{n}h(i)$

我们知道
$\begin{aligned}n=\sum_{d|n}\varphi(d)\end{aligned}$
跟上面式子联想一下，有个 $x^3$ 的关系在里面。

$\begin{aligned}\sum_{d=1}^{n}d^3=\frac{n^2*(n+1)^2}{4}\end{aligned}$
那么我们就有文章可做了。
$\begin{aligned}\sum_{d=1}^{n}d^3&=\sum_{d=1}^{n}\sum_{i|d}\varphi(i)*d^2\\ &=\sum_{d=1}^{n}\sum_{i|d}h(i)(\frac{d}{i})^2\\ &=\sum_{d=1}^{n}\sum_{i=1}^{\frac{n}{d}}h(i)(\frac{d*i}{i})^2.......解释2\\ &=\sum_{d=1}^{n}d^2\sum_{i=1}^{\frac{n}{d}}h(i) \end{aligned}$

解释2：我们之前是枚举每个数d的因子，也就是h(i)中i表示的是d的因子， $\frac{d}{i}$ 表示这个数d的另一个因子。
那么我们就可以等价于：枚举因子t为1到n，那么说明有 $\frac{n}{t}$ 个数a，它们的因子都为t，即根据之前的表达式，h(i)中i还是为因子，此时的另外一个因子就为 $\frac{i*t}{i}$ 。也就是枚举每个数，它对它倍数的贡献。

那么上式就很明显了：

$\begin{aligned}&\frac{n^2*(n+1)^2}{4}=\sum_{i=1}^{n}h(i)+\sum_{d=2}^{n}d^2\sum_{i=1}^{\frac{n}{d}}h(i)\\ 恒等于&S(n)=\frac{n^2*(n+1)^2}{4}-\sum_{d=2}^{n}d^2S(\frac{n}{d})\end{aligned}$

故最终结果就为：

$\begin{aligned}ans=\sum_{d=1}^{n}d*S(\frac{n}{d})\end{aligned}$

我们只需把S()的前6000000预处理出来，再杜教筛一下就求出来了。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<tr1/unordered_map>

using namespace std;
typedef long long LL;
const LL mod=1000000007;
const int N=6e6;

LL sum[N];///存放S的前N项
bool vis[N+10];
LL phi[N+10];
int pri[N+10],tot;

tr1::unordered_map<LL,LL> w;

LL inv2,inv4,inv6;

void init()
{


    phi[1]=1;

    for(int i=2;i<=N;i++)
    {
        if(!vis[i]){
            phi[i]=i-1;
            pri[++tot]=i;
        }
        for(int j=1;j<=tot&&i*pri[j]<=N;j++)
        {
            int t=pri[j];
            vis[i*t]=1;
            if(i%t==0){
                phi[i*t]=phi[i]*t;break;
            }
            phi[i*t]=phi[i]*(t-1);
        }
    }
    
    ///预处理 \sum_{i=1}^{N}i*i*phi[i]
    for(int i=1;i<=N;i++)
    {
        sum[i]=(sum[i-1]+1LL*phi[i]%mod*i%mod*i%mod)%mod;
    }

}

LL fast_pow(LL a1,LL n)
{
    LL ans=1;
    while(n)
    {
        if(n&1) ans=ans*a1%mod;
        a1=a1*a1%mod;
        n>>=1;
    }
    return ans;
}

LL s2(LL x) ///\sum_{i=1}^{x}i*i
{
    LL t=x%mod;
    return t*(t+1)%mod*(2*t+1)%mod*inv6%mod;
}

LL s1(LL x) ///\sum_{i=1}^{x}i
{
    x=x%mod;
    return x*(x+1)%mod*inv2%mod;
}


LL djs(LL x)
{
    if(x<=N) return sum[x];
    if(w[x]) return w[x];

    LL t=x%mod;
    LL ans=t*t%mod*(t+1)%mod*(t+1)%mod*inv4%mod;
    for(LL l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ans=(ans-(s2(r)-s2(l-1)+mod)%mod*djs(x/l)+mod)%mod;
    }

    return w[x]=ans;


}



int main()
{

    LL n;

     inv2=fast_pow((LL)2,mod-2);
     inv4=fast_pow((LL)4,mod-2);
     inv6=fast_pow((LL)6,mod-2);

    init();
    scanf("%lld",&n);


    LL ans=0;///分块求ans
    for(LL l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
     ans=(ans+(s1(r)-s1(l-1)+mod)%mod*djs(n/l)%mod)%mod;
        }

    printf("%lld\n",(ans));
}