约数之和

Problem

d(k)表示k的所有约数的和。d(6) = 1 + 2 + 3 + 6 = 12。

定义S(N) = ∑1<=i<=N ∑1<=j<=N d(i*j)。

例如：S(3) = d(1) + d(2) + d(3) + d(2) + d(4) + d(6) + d(3) + d(6) + d(9) = 59，S(1000) = 563576517282。

给出正整数N，求S(N)，由于结果可能会很大，输出Mod 1000000007(10^9 + 7)的结果。

输入

输入一个数N(2 <= N <= 10^9)。

输出

输出S(N) Mod 1000000007(10^9 + 7)的结果。

输入样例

1000

输出样例

576513341

Guess

首先得会把\(d(ij)\) 写成\(\gcd\) 的形式

然后搞一波莫比乌斯变换

再用和式化一化

上个杜教筛，打个分块

好像就能解决了！

Solve

就是推公式呗。

这里这个\(d(ij)\) 与之前的不太一样，这里是约数之和的函数。

必须要知道这个式子才可以往后推：

\[d(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)==1]x\cdot\frac{j}{y}\]

（意会一下这个式子，大概和之前的除数函数相类似，枚举质因子分析的感觉）
（我不懂怎么来的）

这样的话

\[
\begin{align}
S(n) &= \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}d(ij)\ &= \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)==1]x\cdot\frac{j}{y}\ & 来一波莫比乌斯变换吧！\ &=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{x|i}\sum\limits_{y|j}\sum\limits_{d|\gcd(x,y)}\mu (d)\cdot x\cdot\frac{j}{y}\ &=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{x|i}\sum\limits_{y|j}\sum\limits_{d=1}^n\mu (d)[d|\gcd(x,y)]\cdot x\cdot\frac{j}{y}\ &我们知道d|\gcd(x,y)=d|x 且d|y \&=\sum\limits_{d=1}^n\mu (d)\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{x|i}\sum\limits_{y|j}[d|x][d|y]\cdot x\cdot\frac{j}{y}\&对于枚举1到n的数i的因子x的因子d\&其实就是枚举1到n的因子d的倍数x的倍数i\&=\sum\limits_{d=1}^n\mu (d)\sum\limits_{d|x}\sum\limits_{d|y}\frac{x}{y}\sum\limits_{x|i}\sum\limits_{y|j}j\&这里的\sum\limits_{y|j}j=y+2y+3y+\cdots+\lfloor\frac{n}{y}\rfloor y\&=\sum\limits_{d=1}^n\mu (d)\sum\limits_{d|x}\sum\limits_{d|y}\frac{x}{y}\lfloor\frac{n}{x}\rfloor\cdot y \frac{\lfloor\frac{n}{y}\rfloor(\lfloor\frac{n}{y}\rfloor+1)}{2} \&=\sum\limits_{d=1}^n\mu (d)\sum\limits_{d|x}x\lfloor\frac{n}{x}\rfloor\sum\limits_{d|y}\frac{\lfloor\frac{n}{y}\rfloor(\lfloor\frac{n}{y}\rfloor+1)}{2} \&=\sum\limits_{d=1}^n\mu (d)\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}dx\lfloor\frac{n}{dx}\rfloor\sum\limits_{y=1}^{\lfloor\frac{n}{y}\rfloor}\frac{\lfloor\frac{n}{dy}\rfloor(\lfloor\frac{n}{dy}\rfloor+1)}{2} \&=\sum\limits_{d=1}^n\mu (d)d\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}x\lfloor\frac{n}{dx}\rfloor\sum\limits_{y=1}^{\lfloor\frac{n}{y}\rfloor}\frac{\lfloor\frac{n}{dy}\rfloor(\lfloor\frac{n}{dy}\rfloor+1)}{2} \\end{align}
\]
而对于中间的和号可以看成是\(f(\lfloor\frac{n}{d}\rfloor)\)

\[f(n)=\sum\limits_{x=1}^{n}x\lfloor\frac{n}{x}\rfloor\]

对于后面的和号看成是\(g(\lfloor\frac{n}{d}\rfloor)\)

\[
\begin{align}
g(n)&=\sum\limits_{y=1}^{n}\frac{\lfloor\frac{n}{y}\rfloor(\lfloor\frac{n}{y}\rfloor+1)}{2}\&等差数列前n项和\&=\sum\limits_{y=1}^{n}\sum\limits_{i=1}^{\lfloor\frac{n}{y}\rfloor}i\&=\sum\limits_{y=1}^{n}y\lfloor\frac{n}{y}\rfloor
\end{align}
\]
所以这两个式子是一样的。

\[S(n)=\sum\limits_{d=1}^n\mu (d)\cdot d\cdot (\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}x\lfloor\frac{n}{dx}\rfloor)^2 \]

这样已经可以写了，但是我用的是这个式子

\[S(n)=\sum\limits_{d=1}^n\mu (d)\cdot d\cdot (\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{i=1}^{\lfloor\frac{n}{dx}\rfloor}i)^2 \]

1. 后面的就是分块
2. 前面的用杜教筛（其实就是构造\(h=f*g\)）
3. 外面再套一个分块

补充杜教筛的推导

要求\(A(n)=\sum\limits_{i=1}^ni\cdot \mu(i)\)

我们可以设\(f(i)=i\cdot \mu(i)\)

然后开始构造
\[
\begin{align}
h&=f*g\&=\sum\limits_{d|n}d\cdot \mu(d) \cdot g(n/d)\&可以取g(n/d)=n/d\&=n\sum\limits_{d|n}\mu(d)\&然后我们知道\mu *u=I\&=n[n==1]\&=[n==1]=I
\end{align}
\]
套一下杜教筛的公式
\[A(n)=1-\sum\limits_{i=2}^nd\cdot A(\lfloor\frac{n}{d}\rfloor)\]

Code

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

const int mod = 1e9 + 7;
const int maxn = 6e6;
const int N = 6e6 + 10;
const int inv2 = (mod+1)/2;

int prime[N];
int vis[N];
int mu[N];
ll sum[N];

void sieve()
{
    int p=0;
    mu[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!vis[i]){
            prime[++p]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=p&&i*prime[j]<=maxn;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)break;
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=maxn;i++){
        sum[i]=(sum[i-1]+1ll*mu[i]*i)%mod;
    }
}

unordered_map<int, ll> mp;

ll S(int l, int r){
    return (1ll*(1ll*(l+r)%mod*(r-l+1)%mod)%mod*inv2)%mod;
}

ll djs(int n){
    if(n<=maxn)return sum[n];
    if(mp[n])return mp[n];
    ll ans=1+mod;
    for(int l=2,r;l>=0&&l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans-1ll*S(l,r)%mod*djs(n/l)%mod+mod)%mod;
    }
    return mp[n]=ans;
}

ll F(int x){
    ll ans=0;
    for(int l=1,r;l<=x;l=r+1){
        r=x/(x/l);
        ans+=1ll*(r-l+1)*S(1,x/l)%mod;
    }
    return ans%mod;
}

int main()
{
    int n;
    scanf("%d",&n);
    sieve();
    ll ans=0;
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans%mod+((djs(r)-djs(l-1))%mod)%mod*(F(n/l)%mod*F(n/l)%mod)%mod)%mod;
        ans=(ans+mod)%mod;
    }
    cout<<ans<<endl;
    return 0;
}

原文地址：https://www.cnblogs.com/tongseli/p/11795918.html