在数论中,素数计数函数 π(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