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

Posted 2022-06-16 zhugezy

题意

给出\\(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]\\fracajb \\]
感觉网上那群写博客的都是一眼就能看懂这是啥意思，但我却想了几个小时。。。所以我写详细一点

我们有\\(\\sigma(n)=\\prod_i=1^k(p_i^0+p_i^1+...+p_i^e_i)\\),其中不同的\\(p\\)是出现在\\(n\\)的质因数分解式中不同的质因子，所以我们设\\(i=p_1^a_1p_2^a_2...p_k^a_k,j=p_1^b_1p_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\\),则\\(\\fracajb\\)中\\(p_1\\)的幂次范围是\\(1+b_1\\)到\\(a_1+b_1\\);反过来，当\\(x_2\\)从\\(1\\)往\\(b_1\\)增大时，\\(x_1=0\\),则\\(\\fracajb\\)中\\(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|ii^2j \\\\=\\sum_j=1^n\\sum_j|ii^2j \\\\=\\sum_j=1^n\\sum_a=1^\\lfloor\\fracnj\\rfloor(aj)^2j \\\\=\\sum_i=1^n\\sum_j=1^\\lfloor\\fracni\\rfloorj^2i^3 \\\\=\\sum_i=1^ni^3\\sum_j=1^\\lfloor\\fracni\\rfloorj^2 \\]
回到正题：
\\[ Ans=\\sum_i=1^n\\sum_j=1^n\\sum_a|i\\sum_b|j[gcd(a,b)=1]\\fracajb \\\\=\\sum_a=1^n\\sum_b=1^n\\sum_p=1^\\lfloor\\fracna\\rfloor\\sum_q=1^\\lfloor\\fracnb\\rfloor[gcd(a,b)=1]\\fracaqbb \\\\=\\sum_a=1^n\\sum_b=1^n\\sum_i=1^\\lfloor\\fracna\\rfloor\\sum_j=1^\\lfloor\\fracnb\\rfloor[gcd(a,b)=1]aj \\\\=\\sum_a=1^na\\sum_b=1^n[gcd(a,b)=1]\\sum_i=1^\\lfloor\\fracna\\rfloor\\sum_j=1^\\lfloor\\fracnb\\rfloorj \\\\=\\sum_a=1^na\\lfloor\\fracna\\rfloor\\sum_b=1^n[gcd(a,b)=1]\\frac\\lfloor\\fracnb\\rfloor(\\lfloor\\fracnb\\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\\fracna\\rfloor\\sum_b=1^n\\frac\\lfloor\\fracnb\\rfloor(\\lfloor\\fracnb\\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|ba\\lfloor\\fracna\\rfloor\\frac\\lfloor\\fracnb\\rfloor(\\lfloor\\fracnb\\rfloor+1)2\\mu(d) \\\\=\\sum_d=1^n\\sum_p=1^\\lfloor\\fracnd\\rfloor\\sum_q=1^\\lfloor\\fracnd\\rfloorpd\\lfloor\\fracnpd\\rfloor\\frac\\lfloor\\fracnpd\\rfloor(\\lfloor\\fracnpd\\rfloor+1)2\\mu(d) \\\\=\\sum_d=1^nd\\mu(d)\\sum_p=1^\\lfloor\\fracnd\\rfloorp\\lfloor\\fracnpd\\rfloor\\sum_q=1^\\lfloor\\fracnd\\rfloor\\frac\\lfloor\\fracnqd\\rfloor(\\lfloor\\fracnqd\\rfloor+1)2 \\]
咦，后两个求和符号那两个式子好像是独立的，那就不妨设

\\(f(n)=\\sum_i=1^ni\\lfloor\\fracni\\rfloor\\)

\\(=\\sum_i=1^n\\sum_j=1^ni[i|j]\\) 因为\\(\\lfloor\\fracni\\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\\fracni\\rfloor(\\lfloor\\fracni\\rfloor+1)2\\) 之前好像推出过这个东西，我们再把它变回去试试

\\(=\\sum_i=1^n\\sum_j=1^\\lfloor\\fracni\\rfloorj\\)

\\(=\\sum_j=1^n\\sum_i=1^\\lfloor\\fracnj\\rfloorj\\)

\\(=\\sum_j=1^n\\sum_i=1^nj[j|i]\\)

\\(=f(n)\\).

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

因此有
\\[ Ans=\\sum_d=1^nd\\mu(d)(S_\\sigma(\\lfloor\\fracnd\\rfloor))^2 \\]

5

套杜教筛公式的时间到啦

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

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

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

代码

\\(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;