天天看點
傳回

[軟工課程部落格] 求解第N個素數

任務

求解第 10,0000、100,0000、1000,0000 ... 個素數（要求精确解）。

想法

Sieve of Eratosthenes

學習初等數論的時候曾經學過埃拉托斯特尼篩法(Sieve of Eratosthenes)，這是一種非常古老但是非常有效的求解\(p_n\)的方法，其原理非常簡單：從2開始，将每個素數的各個倍數都标記成合數。

其原理如下圖所示：

[軟工課程部落格] 求解第N個素數

圖引自維基百科

埃拉托斯特尼篩法相比于傳統試除法最大的優勢在于：篩法是将素數的各個倍數标記成合數，而非判定每個素數是否是素數的倍數，使用了加法代替了除法，降低了時間複雜度。

舉個簡單的例子，已知2、3、5、7是素數，求小于49的其餘素數。

試除法

對于每一個數\(m ( 7 < m < 49 )\)而言，都要進行如下判定：

  • 求出 \({\lceil}\sqrt{m}{\rceil} + 1 = n\) ，試除 $ n $以内的素數即可判定 $ m $是否為素數。
  • m / 2 ... 不整除。
  • m / 3 ... 不整除。
  • ...

對于m而言，隻有兩種情況可以結束試除：

  1. m 是素數，要把所有候選素數都試除一遍。
  2. m 是合數，可以被某個素數整除。

對于試除法而言，假設新找到 x 個素數， y 個合數 ，需要試除的素數因子有 m 個，則至少需要做 $ x * m + y $次除法。

埃拉托斯特尼篩法

  • 2，2 + 2 = 4, 4 + 2 = 6 ... 48 + 2 = 50 > 49，下一個。
  • 3，3 + 3 = 6, 6 + 3 = 9 ... 48 + 3 = 51 > 49，下一個。
  • 5，5 + 5 = 10, 10 + 5 = 15 ... 45 + 5 = 50 > 49，下一個。
  • 7, 7 + 7 = 14, 14 + 7 = 21 ... 42 + 7 = 49 = 49，下一個。
  • 沒有小于 8 的素數了，停止标記，沒有标記到的都是素數。

實際上，假設在 \(\sqrt{x}\) 以内有 $ m $個素數，在 \((\sqrt{x},x]\)範圍内又找到了 $n $個素數，可以明顯看到試除法和篩法的差異：

  1. 試除法至少要做 $ n * m + x - n - m (1)$次除法
  2. 篩法需要做 $ x - n - m (2)$次加法

做一次除法運算的時間要大于加法，且有 \((1)\) > \((2)\)，是以試除法開銷遠比篩法大。

現在回到問題本身，問題本身是求解第 n 個素數，是以我們一開始并不知道我們需要在多大的範圍内求解素數，但是許多數學家給了我們很多定理，比如這個第\(n\)個素數\(p_n\)的不等式。

\[p_n < n \ln n + n \ln \ln n\! ( n ≥ 6 )

\]

有了這個函數，我們可以在輸入 $ n $之後用它估算第 $ n $個素數的上界，繼而求解。

優化一：樸素篩法

輸入的 \(n\) 值為素數的序号，假設求得的素數上界為

limit

。為了讓篩法跑得更快，在記憶體允許的範圍内，我們可以直接用一個

limit

大小的數組

sieve[limit]

求解：下标是

int

，值為

bool

，數組初始值均為

true

，表示都未标記。假設下标為

index

。

  • index = 2，sieve[index] = true，開始标記，逐次将 sieve[index+2]标記為false，直到下标超過limit。
  • index = 3，sieve[index] = true，開始标記，逐次将 sieve[index+3]标記為false，直到下标超過limit。
  • index = 4, sieve[index] = false，跳過。

按照上述方式标記完成後，最終在數組中過濾一遍，求得第 n 個未被标記的下标值，即為第 n 個素數。

附代碼如下：

public static void Normal_Sieve(int nth)
{
    int startTime = System.Environment.TickCount;
    int limit = nth < 6 ? 25 : (int)(nth * (Math.Log(nth) + Math.Log(Math.Log(nth))));
    int count = 0;

    List<bool> is_prime = new List<bool>(limit+1);

    for (int i = 0; i < limit+1; i++)
        is_prime.Add(true);

    for (int i = 2; i * i <= limit; i++)
        if (is_prime[i])
            for (int j = i * i; j <= limit; j += i)
                is_prime[j] = false;

    for(int i=2;i<is_prime.Count();i++)
    {
        if (is_prime[i])
            count++;
        if(count == nth)
        {
            Console.WriteLine("The nth_prime is:{0} SpentTime:{1}ms",i,Environment.TickCount- startTime);
            break;
        }
    }
}

優化二：位篩法

位篩法相比于簡單篩法的改進就是：用1個比特位來标記某個下标是否是素數。1個bool類型要占 8 位，用1個比特位可以使程式在極限記憶體容量的情況下，比普通篩法能多計算一些素數。

舉個例子，按照上面的普通篩法，因為記憶體大小的限制最多隻能求解第 1000,000,000 個素數，那麼位篩法就可以求解到第 7000,000,000 個素數。

下面附上位篩法的代碼

public static void Bit_Sieve(int nth)
{
    int startTime = Environment.TickCount;
    int limit = nth < 6 ? 25 : (int)(nth * (Math.Log(nth) + Math.Log(Math.Log(nth))));
    int count = 0;

    int total = limit + 1;
    int sqrt = (int)Math.Sqrt(limit) + 1;

    //[31 30 29 ... 0] every number maps to a bit in uint.
    List<uint> is_prime = new List<uint>((total >> 5) + 1);

    for (int i = 0; i < (total >> 5) + 1; i++)
        is_prime.Add(0x0);


    for (int i = 2; i <= sqrt; i++)
        // is_prime[i>>5] bit i % 32 == 0 means it is a prime
        if ((is_prime[i >> 5] & (1 << (i & 31))) == 0)
        {
            for (int j = i * i; j <= total; j += i)
                // is_prime[j>>5] bit j % 32 = 1;
                is_prime[j >> 5] |= (uint)1 << (j & 31);
        }

    for (int i = 2; i < total; i++)
    {
        if ((is_prime[i >> 5] & (1 << (i & 31))) == 0)
        {
            count++;
            if (count == nth)
            {
                Console.WriteLine("The {0}th_prime is:{1} SpentTime:{2}ms",nth , i, Environment.TickCount - startTime);
                break;
            }
        }
    }
}

在位篩法中大部分運算都是移位、與和或的運算，是以在測試時發現要比簡單的篩法更快一些。

優化三：局部篩法

假設計算得出的第\(n\)個素數上界為\(x\)，在位篩法中，我們申請了一個大小為\(x bit\)的數組用來标記合數。随之而來一個問題，難道在篩法中，我們必須使用 \(x bit\)才能求出第 n 個素數嗎？

當然不是，對于素數上界\(x\)，其實不需要這麼大的空間，我們隻需把\(\sqrt{x}\)以前的素數儲存下來即可。

為什麼一定是 \(\sqrt{x}\)呢？想象一下最暴力的試除法：它在判定一個數是否為素數時，需要周遊試除\(\sqrt{x}\)及以内的素因子。如果\(x\)是一個合數，則必然存在一個素因子整除\(x\)，且其值小于等于\(\sqrt{x}\)。那麼反過來想一下，在篩法中，如果\(\sqrt{x}\)及以内的所有素因子都沒有标記\(x\)為合數，那麼它一定是一個素數了。

基于這樣的思想，我們将\([0,\sqrt{x}]\)内的素數儲存下來，然後對\([\sqrt{x},x]\)分段去掃，每掃一段就記錄一下本段掃到了多少個素數。這樣每次需要載入記憶體的就隻有用于掃描的小素數數組和被掃描的段數組，相比位篩法可以節省更多記憶體空間。

下面舉個簡單的例子來說明這種算法：

假設要求解第11個素數(31)，我們估計出上界為 36，然後下面就是局部篩法的求解過程。

  • \(\sqrt{36}=6\)，然後利用埃拉托斯特尼篩法求出[2,6]内的素數，即數組 [2,3,5]，此時素數有 3 個。
  • 申請一個大小為10的數組 **A **存儲被掃描段。
  • 将[7,36]内的元素按照10個元素一組的方式分成三組（Ps:最後一組不夠10個元素）
  • 初始化數組A，此時數組A的下标[0,9]分别映射到[7,16]。
    • 用素數2開始标記，因為8、10、12...是2的倍數，是以标記A[1]、A[3]、A[5]...為合數。
    • 用素數3開始标記，因為9、12、15...是3的倍數，标記A[2]、A[5]、A[8]...為合數。
    • 用素數5開始标記，因為10、15是5的倍數，标記A[3]、A[8]為合數。
    • 掃描A數組内未被标記元素下标為：0、4、6、10，是以這一段有 4 個素數，把它們對應的實際數字存入素數數組。
    • 現在已經發現了 7 個素數，還未達到預期的11個，繼續掃描。
  • 重新初始化數組A，此時數組A的下标[0,9]分别映射到[17,26]。
    • 與上述過程一樣，用素數2、3、5，分别掃描一遍。
    • 記下未被标記的數字，這一段掃描到了 2個素數，到現在已經發現了 9個素數。
  • 重新初始化數組A，繼續掃描尋找。

局部篩法的求解過程如上所述，因為每個時刻記憶體中隻需要一個段的記憶體來存放需要掃描的段，而不需要一次性把所有的段都加載到記憶體中進行篩選，是以其對記憶體的要求更低。

下面附上局部篩法的代碼

public static void Local_Bit_Sieve(int nth)
{
    int startTime = Environment.TickCount;
    int limit = nth < 6 ? 25 :(int)(nth * (Math.Log(nth) + Math.Log(Math.Log(nth))));
    int sqrt = (int) Math.Sqrt(limit) + 1;
    
    //Get all primes which are less than \sqrt{limit}
    List<uint> isPrime = new List<uint>((sqrt >> 5) +1);
    for (int i = 0; i < (sqrt >> 5) + 1; i++)
        isPrime.Add(0x0);
    for (int i = 2; i * i <= sqrt; i++)
        // is_prime[i>>5] bit i % 32 == 0 means it is a prime
        if ((isPrime[i >> 5] & (1 << (i & 31))) == 0)
        {
            for (int j = i * i; j <= sqrt; j += i)
                // is_prime[j>>5] bit j % 32 = 1;
                isPrime[j >> 5] |= (uint)1 << (j & 31);
        }

    //smallPrimes store the primes
    List<int> smallPrimes = new List<int>();
    for (int i = 2; i < sqrt; i++)
    {
        if ((isPrime[i >> 5] & (1 << (i & 31))) == 0)
        {
            smallPrimes.Add(i);
        }
    }

    int segSize = Math.Max(sqrt,256 * 256);

    uint[] primeSeg = new uint[segSize];
    
    //allPrimes store all primes which are found.
    List<int> allPrimes = new List<int>();
    allPrimes.AddRange(smallPrimes);

    int high = segSize << 5;
    //chunk [2,limit] into different segments
    for (int low = sqrt; low <= limit; low += segSize << 5)
    {
        Array.Clear(primeSeg,0,segSize);
        //for each prime, use them to mark the [low,low + segSize]
        foreach (var sPrime in smallPrimes)
        {
            int initValue = low % sPrime;
            for (int i = (initValue == 0 ? initValue : sPrime - initValue); i < high; i += sPrime)
            {
                primeSeg[i >> 5] |= (uint) 1 << (i & 31);
            }
        }

        for (int i = 0; i < high; i++)
        {
            if ((primeSeg[i >> 5] & (1 << (i & 31))) == 0)
                allPrimes.Add(i + low);
        }

        if (allPrimes.Count() > nth)
        {
            Console.WriteLine("The {0}th_prime is:{1} SpentTime:{2}ms",nth, allPrimes[nth-1],
                Environment.TickCount - startTime);
            break;
        }
    }
}

優化四：局部篩法段優化

我們之前之是以說使用樸素篩法可以跑得飛快，是因為從直覺上說，樸素篩法一次性将所有的資料都加載到記憶體中，一次篩完。而局部篩法卻需要篩取多次，篩取多次好像時間開銷就要比一次加載篩取要多。

看起來局部篩法好像需要篩多次，是以它所耗費的時間就要比一次篩的要慢，但實際上卻不是這樣。實際上，計算機系統中Cache的存在對局部篩法更加有利，所耗費的時間也更小。我們分段篩選，反而更加符合空間上的局部性，是以程式可以更高效地利用Cache，Cache的存取速度要遠大于記憶體，是以局部篩法耗費的時間要比樸素篩更少。這個問題可以這麼形容：是用一個振動頻率比較慢，但是很大的篩子篩選比較快，還是用連續用多個振動頻率較快，但是比較小的篩子篩快？

我們可以通過選擇一個合理的段大小來減少 Cache miss的機率。我電腦上的L1Cache容量是 256KB，是以最後我選擇了 256 * 256 作為段的大小。

Benchmark

n 簡單篩法 位篩法 局部位篩法
1,000,000 1.75s 0.938s 0.5s
2,000,000 4.562s 4.328s 2.156s
3,000,000 7.063s 7.140s 3.140s
5,000,000 13.328s 11.407s 4.750s
8,000,000 22.484s 17.110s 9.140s
10,000,000 23.563s 22.347s 11.078s
20,000,000 57.219s 53.313s 22.734s
目前的測試隻是大緻估計，這個數字不是很準，僅供參考。
N素數問題 nth_prime
上一篇: MySQL中concat函數（連接配接字元串）
下一篇: mySql的UDF是什麼