51Nod1220：约数之和

51Nod1220：约数之和

题意：

\(d(k)\)表示\(k\)所有约数的和。

比如说\(d(6)=1+2+3+6=12\)。

定义：\(S(N)=\sum_{i=1}^N\sum_{j=1}^Nd(i*j)\)。

给出\(N\leq 10^9\)，求\(S(N)\)。

思路：

我们知道：\(\sigma(n)\)表示约数个数和，且为积性函数。

我们需要求\(\sigma(ij)\)，但是\(\sigma\)不是完全积性函数。

这里有个结论：

\[\sigma(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]\frac{yi}{x} \]

代入原式：

\[\sum_{i=1}^N\sum_{j=1}^N\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]\frac{yi}{x} \]

反演：

\[\sum_{i=1}^N\sum_{j=1}^N\sum_{x|i}\sum_{y|j}\frac{yi}{x}\sum_{d|gcd(x,y)}\mu(d)\\\sum_{i=1}^N\sum_{j=1}^N\sum_{x|i}\sum_{y|j}\frac{yi}{x}\sum_{d|x,d|y}\mu(d) \]

改变枚举项为\(xd,yd\)，接着修改为\(id,jd\)：

\[\sum_{d=1}^N\mu(d)\sum_{i=1}^N\sum_{j=1}^N\sum_{xd|i}\sum_{yd|j}\frac{yi}{x}\\\sum_{d=1}^Nd\mu(d)\sum_{i=1}^\frac{N}{d}\sum_{j=1}^\frac{N}{d}\sum_{x|i}\sum_{y|i}\frac{yi}{x} \]

分配一下：

\[\sum_{d=1}^Nd\mu(d)(\sum_{i=1}^\frac{N}{d}\sum_{x|i}\frac{i}{x})(\sum_{j=1}^\frac{N}{d}\sum_{y|i}y) \]

那这后面不就是约数的平方嘛。

\[\sum_{d=1}^Nd\mu(d)\sum_{i=1}^\frac{N}{d}\sigma^2(i) \]

显然后面可以数论分块一下，前面可以杜教筛一下。

先看看怎么样怎样杜教筛这个式子：

设\(f=\mu id,g=id\)，那么有：

\[(f*g)(n)=\sum_{d|n}(\mu(d)d)\frac{n}{d}=n\sum_{d|n}\mu(d)=[n=1] \]

又有：

\[g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S(\frac{n}{i}) \]

可以得到：

\[S(n)=1-\sum_{i=2}^niS(\frac{n}{i}) \]

之后要求这个：

\[\sum_{i=1}^n\sigma(i) \]

怎么计算呢？

需要有点逆向思维，枚举约数不好弄就枚举倍数，所以有：

\[\sum_{i=1}^n\sigma(i)=\sum_{i=1}^ni\frac{n}{i} \]

此时理清楚就可以写代码了。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9+7;
const int maxn = 3e6;
const int inv2 = 500000004;
bool vis[maxn+10];
ll mu[maxn+10];
int primes[maxn+10], cnt;

void init(int n)
{
    mu[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!vis[i])
        {
            primes[++cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1; primes[j] <= n/i; j++)
        {
            vis[primes[j]*i] = 1;
            if(i % primes[j] == 0) break;
            else mu[i*primes[j]] = -mu[i];
        }
    }
    for(int i = 1; i <= n; i++)
    {
        mu[i] = i*mu[i]%mod;
        mu[i] = (mu[i]+mu[i-1])%mod;
    }
}

unordered_map<ll, ll> Smu;

ll getSmu(ll n)
{
    if(n <= maxn) return mu[n];
    if(Smu[n]) return Smu[n];
    ll res = 1;
    for(ll l = 2, r; l <= n; l = r+1)
    {
        r = n/(n/l);
        res -= (r-l+1)%mod*(r+l)%mod*inv2%mod*getSmu(n/l)%mod;
        res = ((res%mod)+mod)%mod;
    }
    return Smu[n] = res;
}

ll n;
unordered_map<ll, ll> G;
ll getG(ll n)
{
    if(G[n]) return G[n];
    ll res = 0;
    for(ll l = 1, r; l <= n; l = r+1)
    {
        r = n/(n/l);
        res = (res + (r-l+1)%mod*(r+l)%mod*inv2%mod*(n/l)%mod)%mod;
    }
    return G[n] = res;
}

int main()
{
    init(maxn);
    cin >> n;

    ll res = 0;
    for(ll l = 1, r; l <= n; l = r+1)
    {
        r = n/(n/l);
        res = (res + ((getSmu(r)-getSmu(l-1))%mod+mod)%mod*getG(n/l)%mod*getG(n/l)%mod)%mod;
    }
    res = ((res%mod)+mod)%mod;
    cout << res << endl;
    return 0;
}