51nod 1220 约数之和（杜教筛 + 推推推推推公式）

题意

给出\(n(1\leq n \leq 10^9)\),求\(\sum_{i=1}^n\sum_{j=1}^n\sigma(ij)\),其中\(\sigma(n)\)表示\(n\)的约数之和。

balabala

交了两道杜教筛的的板子题(51nod 1239, 1244)就看到了这题，然后不会搞，然后看题解看了一天一夜终于彻底搞明白一发A掉了。。。感觉学到了很多，写个博客整理一下，如有错请指出。

技能需求

数论函数与线性筛

莫比乌斯反演（也可以当成容斥去理解）

狄利克雷卷积

杜教筛

强大的数学推导功力（我没有）

先贴参考链接：

题目：http://www.51nod.com/Challenge/Problem.html#problemId=1220

这位这题的题解写的很好: https://www.cnblogs.com/hzoier/p/6937717.html

杜教筛：https://www.cnblogs.com/peng-ym/p/9446555.html 这位写的一看就懂

题解

推推推推推推推，我们记答案为\(Ans\)

1

我们一上来就需要用到一个公式，就是下面这个

\[\sigma(ij)=\sum_{a|i}\sum_{b|j}[gcd(a,b)=1]\frac{aj}{b} \]

感觉网上那群写博客的都是一眼就能看懂这是啥意思，但我却想了几个小时。。。所以我写详细一点

我们有\(\sigma(n)=\prod_{i=1}^k(p_i^0+p_i^1+...+p_i^{e_i})\),其中不同的\(p\)是出现在\(n\)的质因数分解式中不同的质因子，所以我们设\(i=p_1^{a_1}p_2^{a_2}...p_k^{a_k},j=p_1^{b_1}p_2^{b_2}...p_k^{b_k}(a,b\geq 0)\)

则\(\sigma(ij)=\prod_{m=1}^k(p_m^0+p_m^1+...+p_m^{a_m+b_m})\)

现在要想想上面那个公式的含义了，我们不妨考虑\(k=2\)的情况，假设我们选了\(a=p_1^{x_1},b=p_2^{x_2}\)，这时\(gcd(a,b)=1\),那么当\(x_1\)从\(1\)往\(a_1\)增大时，\(x_2=0\),则\(\frac{aj}{b}\)中\(p_1\)的幂次范围是\(1+b_1\)到\(a_1+b_1\);反过来，当\(x_2\)从\(1\)往\(b_1\)增大时，\(x_1=0\),则\(\frac{aj}{b}\)中\(p_1\)的幂次范围是\(b_1-1\)到\(0\).再加上\(x_1=x_2=0\)的情况，我们正好不重不漏的选出了\(p_1\)的所有可能选择。

回头考虑\(\sigma(ij)=\prod_{m=1}^k(p_m^0+p_m^1+...+p_m^{a_m+b_m})\),如果我们把这\(k\)个因式乘出来，变成一个多项式，从这个多项式的角度的某一项去考虑，它有\((a_1+b_1)(a_2+b_2)...(a_k+b_k)\)种选择，选取多项式\(i\)的哪一项不会对选取多项式\(j\)产生任何影响，这就是它们”互相独立“的含义。

2

下面这一段主要思路就是不停的交换求和符号，把整除条件改成枚举条件以化简式子。举个例子，我们要计算下面这个式子：

\[f(n)=\sum_{i=1}^n\sum_{j|i}i^2j \\=\sum_{j=1}^n\sum_{j|i}i^2j \\=\sum_{j=1}^n\sum_{a=1}^{\lfloor\frac{n}{j}\rfloor}(aj)^2j \\=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j^2i^3 \\=\sum_{i=1}^ni^3\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j^2 \]

回到正题：

\[Ans=\sum_{i=1}^n\sum_{j=1}^n\sum_{a|i}\sum_{b|j}[gcd(a,b)=1]\frac{aj}{b} \\=\sum_{a=1}^n\sum_{b=1}^n\sum_{p=1}^{\lfloor\frac{n}{a}\rfloor}\sum_{q=1}^{\lfloor\frac{n}{b}\rfloor}[gcd(a,b)=1]\frac{aqb}{b} \\=\sum_{a=1}^n\sum_{b=1}^n\sum_{i=1}^{\lfloor\frac{n}{a}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{b}\rfloor}[gcd(a,b)=1]aj \\=\sum_{a=1}^na\sum_{b=1}^n[gcd(a,b)=1]\sum_{i=1}^{\lfloor\frac{n}{a}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{b}\rfloor}j \\=\sum_{a=1}^na\lfloor\frac{n}{a}\rfloor\sum_{b=1}^n[gcd(a,b)=1]\frac{\lfloor\frac{n}{b}\rfloor(\lfloor\frac{n}{b}\rfloor+1)}{2} \]

接下来还要搞的话只能在\([gcd(a,b) = 1]\)上做文章了

3

公式:\([gcd(a,b)=1]=\sum_{d=1}^{min(a,b)}\mu(d)[d|gcd(a,b)]\)

原理你可以当成莫比乌斯反演去搞，也可以当成容斥原理来想，（因为这两个本质上是一个东西----huchi）当\(gcd=1\)时，右边等于\(1\)，当\(gcd>1\)时，以\(gcd=12=2^2*3\)为例，右边\(=\mu(1)+\mu(2)+\mu(3)+\mu(4)+\mu(6)+\mu(12)\)，去掉为\(0\)的值，按质因数的个数分个组,\(=(\mu(1))+(\mu(2)+\mu(3))+(\mu(6))=C_2^01^0(-1)^2+C_2^11^1(-1)^1+C_2^21^2(-1)^0=(1-1)^2=0\)

以此类推，如果\(gcd\)有\(x\)个质因子那么就是\((1-1)^x\)，这个公式是这么来的，如果需要严谨的证明的话我不会或许应该有人写过

4

继续。。

\[Ans=\sum_{a=1}^na\lfloor\frac{n}{a}\rfloor\sum_{b=1}^n\frac{\lfloor\frac{n}{b}\rfloor(\lfloor\frac{n}{b}\rfloor+1)}{2}\sum_{d=1}^n\mu(d)[d|gcd(a,b)] \]

老套路，继续交换求和符号，这次要连着\([d|gcd(a,b)]\)一起考虑。\(d|gcd(a,b)\)等价于\(d|a,d|b\).

\[Ans=\sum_{d=1}^n\sum_{d|a}\sum_{d|b}a\lfloor\frac{n}{a}\rfloor\frac{\lfloor\frac{n}{b}\rfloor(\lfloor\frac{n}{b}\rfloor+1)}{2}\mu(d) \\=\sum_{d=1}^n\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{q=1}^{\lfloor\frac{n}{d}\rfloor}pd\lfloor\frac{n}{pd}\rfloor\frac{\lfloor\frac{n}{pd}\rfloor(\lfloor\frac{n}{pd}\rfloor+1)}{2}\mu(d) \\=\sum_{d=1}^nd\mu(d)\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}p\lfloor\frac{n}{pd}\rfloor\sum_{q=1}^{\lfloor\frac{n}{d}\rfloor}\frac{\lfloor\frac{n}{qd}\rfloor(\lfloor\frac{n}{qd}\rfloor+1)}{2} \]

咦，后两个求和符号那两个式子好像是独立的，那就不妨设

\(f(n)=\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor\)

\(=\sum_{i=1}^n\sum_{j=1}^ni[i|j]\) 因为\(\lfloor\frac{n}{i}\rfloor\)意思就是\(1\)到\(n\)中\(i\)的倍数的个数。

\(=\sum_{j=1}^n\sum_{i=1}^ni[i|j]\)

\(=\sum_{j=1}^n\sigma(j)\).

\(g(n)=\sum_{i=1}^n\frac{\lfloor\frac{n}{i}\rfloor(\lfloor\frac{n}{i}\rfloor+1)}{2}\) 之前好像推出过这个东西，我们再把它变回去试试

\(=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j\)

\(=\sum_{j=1}^n\sum_{i=1}^{\lfloor\frac{n}{j}\rfloor}j\)

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

\(=f(n)\).

我们把\(f(n)\)记成\(S_{\sigma}(n)\)吧.

因此有

\[Ans=\sum_{d=1}^nd\mu(d)(S_{\sigma}(\lfloor\frac{n}{d}\rfloor))^2 \]

5

套杜教筛公式的时间到啦

\(\sigma\)是积性函数，可以线筛可以杜教筛，做法很多，我是卷上了个\(\mu\)得到\(id\)来筛的。

记\(f(d)=d\mu(d)\),\(f\)也是个积性函数，可以线筛可以杜教筛，卷上\(id\)得到\(e\),然后杜教筛起来就完事了。

最后求\(Ans\)的时候分块求，需要啥前缀和了就\(Calc\)啥，反正都是\(O(n^{\frac{2}{3}})\)我也不会分析

代码

\(main\)里是主分块，\(Calcmu\)是\(\mu\)的筛，\(Calcf\)是\(f\)的筛，\(Calcs\)是\(\sigma\)的筛

\(359ms\)

#include <bits/stdc++.h>
#define MAXN 1000000
#define LL long long
using namespace std;
const int Mx = 10000;
const int inv2 = 500000004;
const int mod = 1000000007;
LL gcd(LL a, LL b)
{
    return b?gcd(b,a%b):a;
}

int prime[MAXN + 8];
bool notprime[MAXN + 8];
int mu[MAXN + 8];
LL fsum[MAXN + 8], ssum[MAXN + 8], dsum[MAXN + 8];
int ind = 0;
void getprime()
{
    mu[1] = 1;
    ssum[1] = dsum[1] = 1;
    for (int i = 2; i <= MAXN; ++i)
    {
        if(!notprime[i])
        {
            prime[++ind] = i;
            mu[i] = -1;
            ssum[i] = dsum[i] = i + 1;
        }
        for (int j = 1; j <= ind && i * prime[j] <= MAXN; ++j)
        {
            notprime[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                ssum[i * prime[j]] = ssum[i] / dsum[i] * (dsum[i] * prime[j] + 1); dsum[i * prime[j]] = dsum[i] * prime[j] + 1;
                break;
            }
            mu[i * prime[j]] = -mu[i];
            ssum[i * prime[j]] = ssum[i] * (1 + prime[j]); dsum[i * prime[j]] = 1 + prime[j];
        }
    }
    for (int i = 1; i <= MAXN; ++i)
    {
        fsum[i] = (fsum[i - 1] + i * mu[i] % mod + mod) % mod;
        ssum[i] = (ssum[i - 1] + ssum[i] % mod) % mod;
        mu[i] = (mu[i - 1] + mu[i] + mod) % mod;
    }
}
unordered_map<int,LL> _musum, _fsum, _ssum;

LL Calcmu(int n)
{
    if (n <= MAXN)
        return mu[n];
    auto it = _musum.find(n);
    if (it != _musum.end())
        return it->second;
    int last;
    LL ret = 1;
    for (int i = 2; i <= n; i = last + 1)
    {
        last = n / (n / i);
        ret = (ret - (1LL * last - i + 1) * Calcmu(n / i) % mod + mod) % mod;
    }
    return _musum[n] = ret;
}
LL Calcf(int n)
{
    if (n <= MAXN)
        return fsum[n];
    auto it = _fsum.find(n);
    if (it != _fsum.end())
        return it->second;
    int last;
    LL ret = 1;
    for (int i = 2; i <= n; i = last + 1)
    {
        last = n / (n / i);
        ret = (ret - ((1LL * last + i) % mod * (last - i + 1) % mod * inv2 % mod * Calcf(n / i) % mod) + mod) % mod;
    }
    return _fsum[n] = ret;
}

LL Calcs(int n)
{
    if (n <= MAXN)
        return ssum[n];
    auto it = _ssum.find(n);
    if (it != _ssum.end())
        return it->second;
    int last;
    LL ret = 1LL * n * (n + 1) % mod * inv2 % mod;
    for (int i = 2; i <= n; i = last + 1)
    {
        last = n / (n / i);
        ret = (ret - ((Calcmu(last) - Calcmu(i - 1) + mod) % mod * Calcs(n / i) % mod) + mod) % mod;
    }
    return _ssum[n] = ret;
}

int main()
{
    getprime();
    LL n, last;
    cin >> n;
    LL ans = 0;
    for (int i = 1; i <= n; i = last + 1)
    {
        last = n / (n / i);
        LL cc = Calcs(n / i);
        ans = (ans + (Calcf(last) - Calcf(i - 1) + mod) % mod * cc % mod * cc % mod) % mod;
    }
    cout << ans << endl;
    return 0;
}