洛谷 P1403 [AHOI2005]约数研究 题解 【 约数个数、约数和（普通方法+筛法）】

本文转载自： https://www.cnblogs.com/ldysy2012/p/10390857.html

最近本人脑洞大开，发现了一种线性筛约数个数和约数和的神奇方法。
目前网上的方法基本都是利用(num[i])数组记录(i)最小的质因子的个数，然后进行转移。
其实可以省去(num[i])数组，直接进行递推。

设(n)的唯一分解式为：

$n={p_1}^{r_1} {p_2}^{r_2} ... {p_k}^{r_k}$

一、(n)的约数个数公式：

$d(n) = (r_1+1) * (r_2+1) * ... * (r_k+1) $

简写为：

$d(n) = displaystyleprod_{i=1}^k (r_i+1) $

证明也很简单，以(p_1)为例，这个质数因子，可以选择(0)个，可以选择(1)个，...,最多可以选择(r_1)个，就是有(r_1+1)种选择的可能性，其它(p_2,p_3,...,p_k)都是如此，根据乘法原理，所有的可能性就是((r_1+1) * (r_2+1) * ... * (r_k+1))。

举个栗子：
(180= 2^2 * 3^2 * 5)

约数个数(=(1+2) * (1+2) * (1+1) =18)

据说小学奥数比赛中经常考查约数个数与约数和的公式，我也没有参加过奥数比赛，不知道是不是真的。

在给出具体的实现代码之前，还需要了解一个质因子分解的数学性质：

一个正整数的质因子，最多只有一个大于其平方根。

证明：用反证法，假设有超过两个质因子大于其平方根，那么二者相乘一定大于该数。得证。

下面给出朴素版本的约数个数代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;

/**
 * 根据 X =(p1^a1)* (p2^a2)*.......(pn^an)
   T = (a1+1)*(a2+1)*(a3+1)....*(an+1)(a1,a2,a3...an位 p1,p2,p3....pn的系数）
 */

/**
 * 功能：计算约数个数
 * @param n
 * @return
 */
LL getDivisorCount(LL x) {
    unordered_map<int, int> primes; //key:质数 value:个数
    //求质数因子
    for (int i = 2; i <= x / i; i++)
        while (x % i == 0) x /= i, primes[i]++; //primes[i]表示质数i因子的个数+1

    //如果还有质数，那就加上
    if (x > 1) primes[x]++;
    //公式大法
    LL res = 1;
    for (auto p : primes) res = res * (p.second + 1);
    return res;
}

/**
 感悟：
 这是单个处理约数的办法，可以找出所有质数因子，和每个质数因子的个数。再通过公式(1+r1)*(1+r2)*...*(1+rk)，就能知道约数的个数了。
 优点：
 直白，公式，暴力

 缺点：
 不能处理在某个区间内筛出约数的场景，挨个计算，性能差，需要再有一个面对区间的约数筛法。
 */
LL res;

int main() {
    int n;
    cin >> n;
    for (int i = 1; i <= n; i++) cout << i << " " << getDivisorCount(i) << endl;
    return 0;
}

二、(n)的约数和公式：

$sigma(n)=({p_1}^{0}+{p_1}^{1}+{p_1}^{2}+...+{p_1}^{r_1})({p_2}^{0}+{p_2}^{1}+{p_2}^{2}+...+{p_2}^{r_2})...({p_k}^{0}+{p_k}^{1}+{p_k}^{2}+...+{p_k}^{r_k}) $

简写为：

$sigma(n)=prodlimits_{i=1}^{k}(1+p_{i}+p_{i}^{2}+cdots+p_{i}^{r_{i}}) $

Sigma（大写Σ，小写σ），是第十八个希腊字母。

仔细观察一下就知道了，(p_1^0+p_1^1+p_1^2+...+p_1^{r_1})是一个等比数列的求和，有公式滴，不用白不用：

$frac{{p_1}^{a_1+1}-1}{p_1-1} * frac{{p_2}^{a_2+1}-1}{p_2-1} * ... * frac{{p_k}^{a_k+1}-1}{p_k-1}$

简写为：

$prodlimits_{i=1}^{k} frac{p_i^{a_i+1}-1}{p_i-1}$

举个栗子：
(180= 2^2 * 3^2 * 5)

约数和(=(1+2+4) * (1+3+9 ) * (1+5)=546)

下面给出朴素版本求约数和的代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;

/**
 * 功能：计算约数之和
 * @param n
 * @return
 */
LL getSumOfDivisors(LL x) {
    //拆出所有质数因子及质数因子个数
    unordered_map<int, int> primes;
    for (int i = 2; i <= x / i; i++)
        while (x % i == 0) {
            x /= i;
            primes[i]++;
        }
    if (x > 1) primes[x]++;

    //计算约数个数
    LL res = 1;
    for (auto p : primes) {
        LL a = p.first, b = p.second;
        LL t = 1;
        while (b--) t = (t * a + 1);
        res = res * t;
    }
    return res;
}

LL res;

int main() {
    //测试用例：180
    //参考答案：546
    int n;
    cin >> n;
    cout<<getSumOfDivisors(n) << endl;
    return 0;
}

小结：
上面给出了求约数个数、约数和的公式，适用于求指定数字的约数个数、约数和。对于从 (1)~(n)这类问题，效率低下，可以考虑采用类似于筛法的思路解决区间内约数个数与约数和的问题。

三、线性筛法求约数个数

要想线性筛一个积性函数(f(i))，只需要知道4个东西。

扩展知识：
1、积性函数的定义：
积性函数指对于所有互质的整数(a)和(b)有性质(f(a⋅b)=f(a)f(b))的数论函数。

2、积性函数的证明及相关性质
https://blog.csdn.net/WuBaizhe/article/details/76711158

1: (f(1)=?)
　　很显然，我们这里的两个函数：(d(1)=1,σ(1)=1;)
(d(1)=1)意味着数字(1)的约数个数是(1)个。(σ(1)=1)意味着数字(1) 的约数和是(1),只有(1)个(1)，加在一起也是(1)。

2: (f(p)=?)，其中(p)为质数。

　　同样很显然，(d(p)=2,σ(p)=p+1)
(d(p)=2)意味着数字(p)的约数有(2)个，一个是(1)，另一个是就是(p),不可能再有别的约数了，质数嘛。
(σ(p)=p+1) 约数是两个，(1)和(p)，那么约数和是确定的，就是(1+p)。

3、4:(f(i⋅p_j)=?)，(i)与(p_j)互质或不互质。
先说线性筛约数个数的方法:

(1)若(i)与(p_j)互质
我们可以利用积性函数的性质直接推出：(d(i⋅p_j)=d(i)⋅d(p_j)=2d(i))
这里我有两个疑问：
疑问1：为啥(d(p))是积性函数，怎么知道的？
疑问2：为啥(d(p_j)==2)?现在理解就是(p_j)在本文的筛法中是一个质数，(i)是它的倍数，既然(p_j)是质数，那么(d(p_j)=2),这是上面讨论过的内容。

(2)若(i)与(p_j)不互质
考虑(i)，(i * p_j)，(frac{i}{p_j})的关系(由线性筛过程可知,(i)与(p_j)不互质即(p_j|i))

这里不妨设(p_j=p_1,r_j=r_1)：

$i={p_1}^{r_1} {p_2}^{r_2} .... {p_k}^{r_k}$ $i p_1 = {p_1}^{r_1+1} {p_2}^{r_2} .... {p_k}^{r_k}$ $frac{i}{p_1}={p_1}^{r_1-1} {p_2}^{r_2} ... {p_k}^{r_k}$ 则有： $d(i)=(r_1+1) (r_2+1) ... (r_k+1)$ $d(i*p_1)=(r_1+2) (r_2+1) ... (r_k+1)$ $d(frac{i}{p_1})=r_1 (r_2+1) ... (r_k+1)$

设(d(i))中(r_1+1)后面的一大坨为(T),即可表示为：

即: $T=(r_2+1) (r_3+1) ... (r_k+1) $

则：(d(i)=(r_1+1)T) ①

(d(i⋅p_1)=(r_1+2)T=(r_1+1+1)T=(r_1+1)T+T= d(i)+T) ②

(d(frac{i}{p_1})=r_1T=(r_1+1-1)T=(r_1+1)T-T=d(i)−T) ③

将(2、3)式相加，整理得
$d(i⋅p_1) + d(frac{i}{p_1}) = d(i)+T +d(i)−T = 2d(i) $

(d(i⋅p_1)=2d(i)−d(frac{i}{p_1}))

因为最初设(p_j=p_1,r_j=r_1)：所以，就是

(d(i⋅p_j)=2d(i)−d(frac{i}{p_j}))

说白了，费了这么半天劲，就是为了找出这个递推式。

有了递推式，我们来看一下线性筛约数个数的代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;

//用途：线性筛约数个数的代码
const int N = 1e6 + 10;
int d[N];      //约数个数数组,记录每个数字的约数个数
int primes[N]; //质数数组
bool st[N];     //是不是已经被筛掉了
int idx;       //质数数组的下标游标，从1开始，它的最终结果值，就是质数的数量

void getDivisorCount(int n) {
    //1的约数个数是1
    d[1] = 1;
    for (int i = 2; i <= n; i++) { //从2开始到n,进行筛
        if (!st[i])primes[++idx] = i, d[i] = 2; //如果没有被筛掉，那么是质数。记录到质数数组中，同时将d[i]设置为2.因为上文论证了，i为质数时，d[i]=2
        for (int j = 1; i * primes[j] <= n && j <= idx; j++) {//遍历每一个已经出现的质数（欧拉筛大法好！），对primes[j]的整数倍进行筛除,i * primes[j] <= n表示是范围内，超过边界就不用了
            st[i * primes[j]] = true;//标识已筛出
            //i与primes[j]不互质
            if (i % primes[j] == 0) {
                d[i * primes[j]] = d[i] * 2 - d[i / primes[j]];//利用费了牛劲才得出的递推公式进行计算，因为是从小到大过来的，依赖项提前算好，所以可以算出来
                break;  //都说这个break值钱，只为第一个质数因子筛掉，欧拉筛的精髓，妙！
            }
            //i与primes[j]互质
            d[i * primes[j]] = d[i] * 2; //积性函数的性质直接推出,在上文中可以找到理由
        }
    }
}

LL res;

int main() {
    int n = 100;
    //线性筛约数个数
    getDivisorCount(n);
    for (int i = 1; i <= n; i++) res += d[i];
    cout << res << endl;
    return 0;
}

四、线性筛法求约数和
1、(i) 与(p_j)互质时

$sigma(icdot p_{j})=sigma(i)cdot sigma(p_{j})=sigma(i)(p_{j}+1)$ 解读：因为在本文的上下文中，$p_j$是个质数，而质数的$sigma(p)=p+1$，因为它只有两个约数，$1$和$p$,和当然是$p+1$了。

2、(i) 与(p_j)不互质时
这里对于上下文中，如果不互质，那么肯定是(p_j|i),不存在其它有公约数的情况。
这时考虑(sigma(i)),(sigma(frac{i}{p_{1}})),(sigma(icdot p_{1}))三者之间的关联性。

疑问：
你要说考虑(sigma(i)),(sigma(icdot p_{1}))的关系，我能理解，因为递推嘛，肯定要找两者之间的关系，怎么想到与sigma(frac{i}{p_{1}})也有关联呢？这是大神的思路啊？还是有其它不为人知的推导过程，我没看过，所以孤陋寡闻了~

设：

$i=p_{1}^{r_{1}}p_{2}^{r_{2}}cdots p_{k}^{r_{k}}$ $icdot p_{1}=p_{1}^{r_{1}+1}p_{2}^{r_{2}}cdots p_{k}^{r_{k}}$ $frac{i}{p_{1}}=p_{1}^{r_{1}-1}p_{2}^{r_{2}}cdots p_{k}^{r_{k}}$

则有：
(sigma(i)=(1+p_{1}+cdots+p_{1}^{r_{1}})(1+p_{2}+cdots+p_{2}^{r_{2}})cdots(1+p_{k}+cdots+p_{k}^{r_{k}}))

(sigma(icdot p_{1})=(1+p_{1}+cdots+p_{1}^{r_{1}+1})(1+p_{2}+cdots+p_{2}^{r_{2}})cdots(1+p_{k}+cdots+p_{k}^{r_{k}}))

(sigma(frac{i}{p_{1}})=(1+p_{1}+cdots+p_{1}^{r_{1}-1})(1+p_{2}+cdots+p_{2}^{r_{2}})cdots(1+p_{k}+cdots+p_{k}^{r_{k}}))

同理：设后面的那一大串为(T).
即：(T=(1+p_{2}+cdots+p_{2}^{r_{2}})cdots(1+p_{k}+cdots+p_{k}^{r_{k}}))
则：
(sigma(i)=(1+p_{1}+cdots+p_{1}^{r_{1}})T)
(sigma(icdot p_{1})=(1+p_{1}+cdots+p_{1}^{r_{1}+1})T=(1+p_{1}+cdots+p_{1}^{r_{1}})T + p_{1}^{r_{1}+1}T = sigma(i)+p_{1}^{r_{1}+1}T)
(sigma(frac{i}{p_{1}})=(1+p_{1}+cdots+p_{1}^{r_{1}-1})T=(1+p_{1}+cdots+p_{1}^{r_{1}-1} +p_{1}^{r_{1}} -p_{1}^{r_{1}})T =(1+p_{1}+cdots+p_{1}^{r_{1}-1} +p_{1}^{r_{1}})T-p_{1}^{r_{1}}T =sigma(i)-p_{1}^{r_{1}}T)

整理一下：
(sigma(i)=(1+p_{1}+cdots+p_{1}^{r_{1}})T) ①
(sigma(icdot p_{1})=(1+p_{1}+cdots+p_{1}^{r_{1}+1})T=sigma(i)+p_{1}^{r_{1}+1}T) ②
(sigma(frac{i}{p_{1}})=(1+p_{1}+cdots+p_{1}^{r_{1}-1})T=sigma(i)-p_{1}^{r_{1}}T) ③

为了消元，去掉(T),所以③*(p_1)+②
得到：
(sigma(icdot p_{1})+p_{1}sigma(frac{i}{p_{1}})=(p_{1}+1)sigma(i))
整理，即：
(large{sigma(icdot p_{1})=(p_{1}+1)sigma(i)-p_{1}sigma(frac{i}{p_{1}})})

再次费了牛劲，得到了第二个递推关系式！

下面给出代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;

//用途：线性筛约数和的代码
const int N = 1e6 + 10;
int primes[N];  //质数数组
int idx;        //质数数组下标游标
bool st[N];     //是否被筛出
int sigma[N];   //约数和数组

void getDivisorSum(int n) {
    sigma[1] = 1; //σ(1)=1,因为1的约数只有1，约数和就是1，这是递推的起点
    for (int i = 2; i <= n; i++) { //倍数
        if (!st[i])primes[++idx] = i, sigma[i] = i + 1; //如果是质数，那么放入到质数数组中，并且质数的约数和是i+1
        for (int j = 1; i * primes[j] <= n && j <= idx; j++) {//遍历每一个已经出现的质数（欧拉筛大法好！），对primes[j]的整数倍进行筛除,i * primes[j] <= n表示是范围内，超过边界就不用了
            st[i * primes[j]] = true;   //标识为质数
            if (i % primes[j] == 0) {    //i与primes[j]不互质
                //利用费了牛劲才得出的递推公式进行计算，因为是从小到大过来的，依赖项提前算好，所以可以算出来
                sigma[i * primes[j]] = sigma[i] * (primes[j] + 1) - primes[j] * sigma[i / primes[j]];
                break;//都说这个break值钱，只为第一个质数因子筛掉，欧拉筛的精髓，妙！
            }
            //i与primes[j]互质
            sigma[i * primes[j]] = sigma[i] * (primes[j] + 1);//约数和是当前质数+1，因为一共两个约数，一个是1一个是自己，和当然是p+1
        }
    }
}

LL res;

int main() {
    int n = 100;
    //线性筛约数和
    getDivisorSum(n);
    for (int i = 1; i <= n; i++) res += sigma[i];
    cout << res << endl;
    return 0;
}

事实上它还可以再短一点(附上约数个数和约数和放在一起的版本)：

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;

/**
 * 功能：线性筛出约数个数与约数和
 * Tag:模板，约数个数，约数和
 https://www.cnblogs.com/littlehb/p/14944596.html
 */
const int N = 1e6 + 10;
int n;
int primes[N];      //质数数组
int idx;            //质数数组下标游标
bool st[N];          //是否已被筛出
int d[N];           //约数个数数组
int sigma[N];       //约数和数组

void get_divisor(int n) {
    //积性函数的出发值
    d[1] = sigma[1] = 1;

    for (int i = 2; i <= n; i++) {  //倍数
        if (!st[i])primes[++idx] = i, d[i] = 2, sigma[i] = i + 1;
        for (int j = 1; i * primes[j] <= n & j <= idx; j++) {
            st[i * primes[j]] = true;
            d[i * primes[j]] = d[i] << 1;
            sigma[i * primes[j]] = sigma[i] * (primes[j] + 1);
            if (i % primes[j] == 0) {
                d[i * primes[j]] -= d[i / primes[j]];
                sigma[i * primes[j]] -= primes[j] * sigma[i / primes[j]];
                break;
            }
        }
    }
}

LL res;

int main() {
    cin >> n;
    //开始筛约数个数，约数和
    get_divisor(n);
    //输出约数个数和
    for (int i = 1; i <= n; i++) res += d[i];
    cout << res << endl;

    //输出约数和
    cout << sigma[n] << endl;

    return 0;
}