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;
}