洛谷 P1403 [AHOI2005]約數研究題解【約數個數、約數和（普通方法+篩法）】

時間 2021-08-14

標籤 c++ 數組 ide 函數性能測試 spa get 數學欄目 C&C++ 简体版

原文原文鏈接

最近本人腦洞大開，發現了一種線性篩約數個數和約數和的神奇方法。
目前網上的方法基本都是利用$num[i]$數組記錄$i$最小的質因子的個數，而後進行轉移。
其實能夠省去$num[i]$數組，直接進行遞推。c++

設$n$的惟一分解式爲：數組

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

1、$n$的約數個數公式：ide

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

簡寫爲：函數

$d(n) = \displaystyle\prod_{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$spa

聽說小學奧數比賽中常常考查約數個數與約數和的公式，我也沒有參加過奧數比賽，不知道是否是真的。ci

在給出具體的實現代碼以前，還須要瞭解一個質因子分解的數學性質：get

一個正整數的質因子，最多隻有一個大於其平方根。

證實：用反證法，假設有超過兩個質因子大於其平方根，那麼兩者相乘必定大於該數。得證。數學

下面給出樸素版本的約數個數代碼：

#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;
}

2、$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)=\prod\limits_{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}$

簡寫爲：

$\prod\limits_{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$這類問題，效率低下，能夠考慮採用相似於篩法的思路解決區間內約數個數與約數和的問題。

3、線性篩法求約數個數

要想線性篩一個積性函數$f(i)$，只須要知道4個東西。

擴展知識：
一、積性函數的定義：
積性函數指對於全部互質的整數$a$和$b$有性質$f(a⋅b)=f(a)f(b)$的數論函數。

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$。

三、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$ ③

將$二、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;
}

4、線性篩法求約數和
一、$i$ 與$p_j$互質時

$\sigma(i\cdot p_{j})=\sigma(i)\cdot \sigma(p_{j})=\sigma(i)(p_{j}+1)$ 解讀：由於在本文的上下文中，$p_j$是個質數，而質數的$sigma(p)=p+1$，由於它只有兩個約數，$1$和$p$,和固然是$p+1$了。

二、$i$ 與$p_j$不互質時
這裏對於上下文中，若是不互質，那麼確定是$p_j|i$,不存在其它有公約數的狀況。
這時考慮$\sigma(i)$,$\sigma(\frac{i}{p_{1}})$,$\sigma(i\cdot p_{1})$三者之間的關聯性。

疑問：
你要說考慮$\sigma(i)$,$\sigma(i\cdot p_{1})$的關係，我能理解，由於遞推嘛，確定要找二者之間的關係，怎麼想到與\sigma(\frac{i}{p_{1}})也有關聯呢？這是大神的思路啊？仍是有其它鮮爲人知的推導過程，我沒看過，因此孤陋寡聞了~

設：

$i=p_{1}^{r_{1}}p_{2}^{r_{2}}\cdots p_{k}^{r_{k}}$ $i\cdot 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(i\cdot 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(i\cdot 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(i\cdot 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(i\cdot p_{1})+p_{1}\sigma(\frac{i}{p_{1}})=(p_{1}+1)\sigma(i)$
整理，即：
$\large{\sigma(i\cdot 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:模板，約數個數，約數
 */
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;
}