素數計數函式

黃志斌發表於2016-11-24

在數論中,素數計數函式 π(x) 返回不超過 x 的素數的個數。

可以使用以下 C# 語言程式來計算 π(x):

 1: using System;
 2: 
 3: namespace Skyiv.Utils
 4: {
 5:   sealed class PrimeHelper
 6:   {
 7:     long m; int z; int[] a; long[] b;
 8: 
 9:     public PrimeHelper(long n)
10:     {
11:       a = new int[(z = (int)Math.Sqrt(m = n)) + 1];
12:       b = new long[z + 1];
13:       for (int i = 1; i <= z; a[i] = i-1, ++i) b[i] = n/i-1;
14:       for (int p = 2; p <= z; ++p) {
15:         if (a[p] <= a[p - 1]) continue;
16:         long d = p, q = (long)p * p;
17:         int c = a[p - 1], m = (int)Math.Min(z, n / q);
18:         for (int i = 1; i <= m; ++i, d += p)
19:           b[i] += c - ((d <= z) ? b[d] : a[n / d]);
20:         for (int i = z; i >= q; --i) a[i] += c - a[i/p];
21:       }
22:     }
23: 
24:     public long π(long n)
25:     { // 當 n > sqrt(m) 時,僅保證對 n == floor(m/i) 正確
26:       return (n > z) ? b[m / n] : a[n];
27:     }
28: 
29:     public int[] GetPrimes(bool hasSentinel = false)
30:     { // 返回不超過 sqrt(m) 的所有素數
31:       var primes = new int[a[z] + (hasSentinel ? 1 : 0)];
32:       for (int j = 0, i = 2; i <= z; ++i)
33:         if (a[i] > a[i - 1]) primes[j++] = i;
34:       if (hasSentinel) primes[a[z]] = z + 1;
35:       return primes;
36:     }
37:   }
38: }

簡要說明:

  • 第 24 至 27 行的 π(n) 方法返回不超過 n 的素數的個數。當 n > 時,僅保證對 正確,其中 i = 1,2,...。
  • 第 29 至 36 行的 GetPrimes 函式返回不超過 的所有素數。當 hasSentinel 為真時,尾部附加有一個哨兵(該哨兵的值大於 ,但不一定是素數)。這個方法對於計算 π(n) 來說不是必需的,僅是順手提供一個可能有用的方法。也就是說,刪除這個方法並不影響計算 π(n)。

參考資料

  1. Wikipedia: Prime-counting function

相關文章