質數篩法

參考資料：https://oi-wiki.org/math/number-theory/sieve/ 

何謂篩法？爲找出 1~n 的所有質數，而欲將其中合數除去，此爲篩也。而 1~n 中合數必有質因數處 1~sqrt(n) 中，故欲盡求 1~n  中合數，只需將 1~sqrt(n) 中質數之倍數小於等於 n 者盡求，則合數盡得矣。

由是，得埃氏篩法如下：

#include <vector>
using namespace std;

#define N 10000

vector<int> prime;    // 確定下來的質數
bool is_prime[N];    // 變動之中的質數列表，將初始化爲全 true，再將其中合數設爲 false

// Eratosthenes 篩法（埃拉託斯特尼篩法，簡稱埃氏篩法）
void Eratosthenes(int n)
{
    is_prime[0] = is_prime[1] = false;    // 0, 1 非質非合
    for (int i = 2; i <= n; ++i)
        is_prime[i] = true;
    
    for (int i = 2; i <= n; ++i)
    {
        if (is_prime[i])    // 先前循環未標記此數爲合數，則此數必是素數
        {
            prime.push_back(i);
            
            // 此處不能直接 break 因爲後面的質數還沒有放入 vector
            if ((long long)i * i > n)
                continue;
            
            // 對於 i=3 而言，它的 2 倍已經在 i=2 時篩過了；
            // 對於 i=4 而言，它的 2, 3 倍已經在 i=2,3 時篩過了；
            // 故對 i 而言，直接篩它的 i, i+1, i+2, ... 倍即可
            for (int j = i * i; j <= n; j += i)
                is_prime[j] = false;  // 是 i 的倍数的均不是素数
        }
    }
}

int main()
{
    Eratosthenes(N - 1);

    for (const auto& var : prime)
    {
        printf("%d ", var);
    }

    return 0;
}

稍優化之，寫法略美觀些：

vector<int> prime;
bool is_prime[N];

// Eratosthenes 筛法（埃拉托斯特尼筛法，简称埃氏筛法）
void Eratosthenes(int n)
{
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i <= n; ++i) is_prime[i] = true;

    // 只需遍歷 2~sqrt(n) 中質數，便可篩完 2~n 之合數
    // i * i <= n 等價於 i <= sqrt(n)
    for (int i = 2; i * i <= n; ++i)
        if (is_prime[i])
            for (int j = i * i; j <= n; j += i) is_prime[j] = false;

    // 合數篩畢，素數猶未記錄，此處補錄
    for (int i = 2; i <= n; ++i)
        if (is_prime[i]) prime.push_back(i);
}

至此，埃氏篩法已可用矣，然此法也有重篩之處，尚待優化。

例如，12 = 2 * 6 = 3 * 4，在埃氏篩之中

i = 2 時篩了：（從 i*i 開始篩）2^2 = 4, 6, 8, 10, 12, ...

i = 3 時篩了：（從 i*i 開始篩）3^3 = 9, 12, ...

發現 12 作爲一個合數被它的兩種不同分解方式篩了兩遍，浪費了時間。

對於 18, 24, 30 也是如此：

18 = 2 * 9 = 3 * 6，也會被 2 和 3 重複篩

30 = 2 * 15 = 3 * 10 = 5 * 6 會被 2, 3, 5 重複篩

數值更大時，重複現象便更爲明顯。重複之根源在於許多合數的分解方式不唯一，且可分解爲 p * i (p 爲質數，i >= p) 的形式，故而會被埃氏篩法的循環重複標記。

要解決重複，就應當採取唯一的分解方式。故此後對每個合數，只關注它的 p 最小的 p * i 分解，就可以避免重複。這就是歐拉篩法，具體原理請看代碼註釋：

#include <vector>

// 歐拉篩
std::vector<int> EulerPrimes(int n)
{
    if (n < 2)
        return {};

    std::vector<int> primes;                    // 存儲質數
    std::vector<bool> is_prime(n + 1, true);    // 臨時數組，假定所有數均爲質數，後續再標記出其中合數
    is_prime[0] = is_prime[1] = false;

    // 以 i 遍歷 2~n，既代表當前判斷到的質數上界，又代表合數 p*i 分解中的 i
    for (int i = 2; i <= n; i++)
    {
        if (is_prime[i])            // 既未被前面循環記爲合數，則必是質數。起始數 2 亦復如是。
            primes.push_back(i);

        // 此時已記錄的所有質數 p 均 <= i，因此可對每個 p 構造合數 p*i 進行標記，但須保證 p*i 是該合數的 p 最小的分解。 
        // 由於 p 在 primes 中從小到大存儲，故一旦發現 i 可被 p 整除，便說明 i 有最小質因數 p。
        // 此時的 p*i 結構依然是合數 p*i 的 p 最小的分解（因爲此式分解不出更小的質數），但對於下一個從 primes 中取出的 p'，
        // 由於 p' > p，故 p'*i 必然可以寫作 p'*p*(i/p) 即 p*(p'*i/p)，屬於 p*i 結構而非 p'*i 結構，下個 p''>p 亦是如此，故應退出。
        // 此乃減少重複之關鍵。
        for (const int& p : primes)
        {
            long long not_prime = (long long)p * i;
            if (not_prime > n)    // 超出所求質數範圍的合數無需標記。後面由於 p 越來越大，因此循環可結束矣
                break;
            is_prime[not_prime] = false;
            if (i % p == 0)
                break;
        }
    }

    return primes;
}

// Eratosthenes 篩法（埃拉託斯特尼篩法，簡稱埃氏篩法）
std::vector<int> EratosthenesPrimes(int n)
{
    if (n < 2)
        return {};

    std::vector<int> primes;                    // 存儲質數
    std::vector<bool> is_prime(n + 1, true);    // 臨時數組
    is_prime[0] = is_prime[1] = false;

    // 只需遍歷 2~sqrt(n) 中質數，便可篩完 2~n 之合數
    // i * i <= n 等價於 i <= sqrt(n)
    for (int i = 2; i * i <= n; ++i)
        if (is_prime[i])
            for (int j = i * i; j <= n; j += i) is_prime[j] = false;

    // 合數篩畢，素數猶未記錄，此處補錄
    for (int i = 2; i <= n; ++i)
        if (is_prime[i]) primes.push_back(i);

    return primes;
}

int main()
{
    std::vector<int> p1 = EulerPrimes(100);
    std::vector<int> p2 = EratosthenesPrimes(100);

    for (const auto& var : p1)
    {
        printf("%d ", var);
    }
    printf("\n");
    for (const auto& var : p2)
    {
        printf("%d ", var);
    }

    return 0;
}

經粗略測定知，輸入 n = 10^7 時，歐拉篩比埃氏篩法快三倍。但是數據量小的時候優勢不是非常明顯。

歐拉函數

定義：在數論，對正整數 n，歐拉函數是小於 n 的正整數中與 n 互質的數的數目。

我的描述：euler_phi(n) = sum( gcd(n, k) == 1 , k : 0 -> n-1)

歐拉函數 phi(0) ~ phi(100) 的值：

0 1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20 12 18 12 28 8 30 16 20 16 24 12 36 18 24 16 40 12 42 20 24 22 46 16 42 20 32 24 52 18 40 24 36 28 58 16 60 30 36 32 48 20 66 32 44 24 70 24 72 36 40 36 60 24 78 32 54 40 82 24 64 42 56 40 88 24 72 44 60 46 72 32 96 42 60 40

只要會了歐拉篩，並且瞭解一些歐拉函數的性質，就可以通過稍微修改上述歐拉篩代碼實現求解歐拉函數。

我通過一篇博客學習了歐拉函數的幾條重要性質，由此完成了下面的歐拉公式求解代碼。博客連接：https://www.cnblogs.com/JustinRochester/p/12230413.html 

需要用到的歐拉函數性質

打字不好備述，若需要查看請看圖片：

https://s2.loli.net/2024/12/14/wF3jBpJiPOIl9Kr.jpg 
https://s2.loli.net/2024/12/14/MIVag7EtsOBn3j1.jpg 
https://s2.loli.net/2024/12/14/aKNMl8sDFOcBfIQ.png 

計算 0~n 的歐拉函數值的代碼（基本上就是歐拉篩的代碼）：

#include <vector>

// 求 0~n 的歐拉函數值
std::vector<int> EulerPhi(int n)
{
    if (n < 2)
        return {};

    std::vector<int> primes;                    // 存儲質數
    std::vector<int> euler_phi(n + 1, 0);        // 歐拉函數值
    std::vector<bool> is_prime(n + 1, true);    // 臨時數組，假定所有數均爲質數，後續再標記出其中合數
    is_prime[0] = is_prime[1] = false;
    euler_phi[1] = 1;

    for (int i = 2; i <= n; i++)
    {
        if (is_prime[i])
        {
            primes.push_back(i);
            euler_phi[i] = i - 1;    // 質數 p 的歐拉函數值爲 p-1
        }

        for (const int& p : primes)
        {
            long long not_prime = (long long)p * i;
            if (not_prime > n)
                break;

            is_prime[not_prime] = false;
            if (i % p == 0)
            {
                euler_phi[not_prime] = p * euler_phi[i];        // p, i 不互質的情況
                break;
            }
            else
            {
                euler_phi[not_prime] = (p - 1) * euler_phi[i];    // p, i 互質的情況
            }
        }
    }

    return euler_phi;
}

int main()
{
    std::vector<int> p1 = EulerPhi(100);

    for (const auto& var : p1)
    {
        printf("%d ", var);
    }

    return 0;
}