51nod [1220 约数之和]

题意

题意：已知 $ \sigma(i)=\sum\limits_{d|i}d $ ，求 $ \sum\limits_ {i=1}^n\sum\limits_ {j=1}^n\sigma(ij) $ ，其中 $ 1\le n\le10^9 $。

题目链接

题解

首先我们先来看一个结论： $ \sigma(ij)=\sum\limits_ {x|i}\sum\limits_ {y|j} xy\times\left[\gcd\left(\frac{i}{x},y\right)=1\right] $ 。其中 $ [x] $ 表示在 $ x $ 成立的时候这个值为 $ 1 $ ，不成立时这个值为 $ 0 $ 。

为啥这个结论是对的呢？前面的 $ xy $ 实际上是在枚举这个约数的值，但是为什么在 $ xy $ 后面还要乘上一个数字呢？

因为如果不乘那个数字会重复计算，例如在 $ i=6,j=3 $ 时 $ x=1,y=3 $ 和 $ x=3,y=1 $ 这两组取值都可以得到一个 $ ij $ 的约数 $ 3 $ ，这时候这个约数就被计算了两遍。加上这个约束就不会重复计算了。

下来我们来推式子。

欢乐的颓柿子时间

\[\begin{aligned} &\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sigma(ij)\\ =&\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_ {x|i}^i\sum\limits_ {y|j}^j xy\times\left[\gcd\left(\frac{i}{x},y\right)=1\right]&\textbf{展开上面的式子}\\ =&\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_ {x|i}^i\sum\limits_ {y|j}^j \frac{iy}{x}\times\left[\gcd\left(x,y\right)=1\right]&\textbf{这里换了一下枚举的方式}\\ =&\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_ {x|i}^i\sum\limits_ {y|j}^j \frac{iy}{x}\times\varepsilon\left(\gcd\left(x,y\right)\right)&\varepsilon(n)\textbf{和}[n=1]\textbf{等价}\\ =&\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_ {x|i}^i\sum\limits_ {y|j}^j \frac{iy}{x}\times\sum\limits_{d|\gcd(x,y)}\mu(d)&\textbf{由迪利克雷卷积知}\mu\ast I=\varepsilon\\ =&\sum_{d=1}^n\mu(d)\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_ {x|i\wedge d|x}^i\kern{5pt}\sum\limits_ {y|j\wedge d|y}^j \frac{iy}{x}&\textbf{换一下枚举顺序}\\ =&\sum_{d=1}^nd\mu(d)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_ {x|i}\sum\limits_ {y|j} \frac{iy}{x}& i,j,x,y\textbf{都是 $ d $ 的倍数，我们枚举他们是 $ d $ 的几倍}\\ =&\sum_{d=1}^nd\mu(d)\left(\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_ {x|i}\frac{i}{x}\right)\left(\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_ {y|j}y\right)&\textbf{换一下枚举顺序}\\ =&\sum\limits_{d=1}^nd\mu(d)\left(\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sigma(i)\right)^2&\textbf{然后会发现一些神奇的东西}\\ \end{aligned} \]

令 $ h(n)=\left(\sum\limits_ {i=1}^ n\sigma(i)\right)^2 $

那么上面的式子可以变成 $ \sum\limits_{d=1}^nd\mu(d)h\left(\frac{n}{d}\right)$ 。

可以发现，如果可以比较快速求出 $ h(i) $ 的值以及 $ d\mu(d) $ 的前缀和就可以对整个式子进行数论分块。

Part1

我们先来看一下 $ h $ 的值如何求。

\[\begin{aligned} &h(n)\\ =&\left(\sum\limits_{i=1}^n\sigma(i)\right)^2\\ =&\left(\sum\limits_{i=1}^n\sum_{j|i}j\right)^2\\ =&\left(\sum_{j=1}^nj\left\lfloor\frac{n}{j}\right\rfloor\right)^2\\ \end{aligned} \]

这个玩意也可以数论分块（瞬间窒息

Part2

我们设 $ f(i)=i\mu(i) $ 。

那么我们要求的 $ S(i)=\sum\limits_{j=1}^if(j) $ 。

既然要在 $ 10^9 $ 范围内求一个函数（就是 $ f $ ）的前缀和，并且这还是一个积性函数，我们首先考虑杜教筛。

不难发现 $ f=id\cdot \mu $ ，那我们设 $ g=id $ 。

让我们把杜教筛的式子写出来：

\[S(n)g(1)=\sum\limits_{i=1}^n(f\ast g)(i)-\sum\limits_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)g(i) \]

$ g $ 的前缀和挺好求，不说了。

我们来推一下 $ f\ast g $ 的。

\[\begin{aligned} &\sum_{i=1}^n(f\ast g)(i)\\ =&\sum_{i=1}^n\sum_{d|i}d\times\mu(d)\times \frac{i}{d}\\ =&\sum_{i=1}^ni\sum_{d|i}\mu(d)\\ =&\sum_{i=1}^ni\times(\mu\ast I)(i)\\ =&\sum_{i=1}^ni\times\varepsilon(i)\\ =&1 \end{aligned} \]

Game over!

代码

#include<cstdio>
#include<map>
using namespace std;
const int mod=1000000007;
using ll=long long;
inline int read()
{
    int s=0,w=1;char ch;
    while((ch=getchar())>'9'||ch<'0') if(ch=='-') w=-1;
    while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
    return s*w;
}
int pri[1000001];
bool is[1000001];
ll f[1000001];
map<int,ll>mf;
ll ans;
int n;
ll run(int l,int r)
{
    return (ll)(r-l+1)*(l+r)/2%mod;
}
ll get(int num)
{
    if(num<=1000000) return f[num];
    if(mf.count(num)) return mf[num];
    // printf("get %d\n",num);
    ll val=1;
    for(int l=2,r;l<=num;l=r+1)
    {
        r=num/(num/l);
        val=(val-run(l,r)*get(num/l))%mod;
    }
    return mf[num]=(val+mod)%mod;
}
ll getpre(int num)
{
    ll ans=0;
    for(int l=1,r;l<=num;l=r+1)
    {
        r=num/(num/l);
        ans=(ans+run(l,r)*(num/l))%mod;
    }
    return ans*ans%mod;
}
int main()
{
    n=1000000,f[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!is[i]) pri[++pri[0]]=i,f[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=n;j++)
        {
            is[i*pri[j]]=1;
            if(i%pri[j]==0) break;
            f[i*pri[j]]=-f[i];
        }
    }
    for(int i=1;i<=n;i++) f[i]=(f[i]*i+f[i-1]+mod)%mod;
    n=read();
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=(ans+(get(r)-get(l-1))*getpre(n/l))%mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}