题目

Problem 609. π sequences

For every n ≥ 1 the prime-counting function π(n) is equal to the number of primes not exceeding n.
E.g. π(6) = 3 and π(100) = 25.

We say that a sequence of integers u = (u0, …, um) is a π sequence if

  • un ≥ 1 for every n
  • un + 1 = π(un)
  • u has two or more elements

For u0 = 10 there are three distinct π sequences: (10,4), (10,4,2) and (10,4,2,1).

Let c(u) be the number of elements of u that are not prime.
Let p(n,k) be the number of π sequences u for which u0 ≤ n and c(u) = k.
Let P(n) be the product of all p(n,k) that are larger than 0.
You are given: P(10) = 3×8×9×3 = 648 and P(100) = 31038676032.

Find P(108). Give your answer modulo 1000000007.

分析

计算 P(10) 的过程如下图所示:

  • 黄色的行对应 u0 为素数的情况
  • 青色的行对应 u0 为合数的情况
  • 青色的行的 c(u) 值刚好比其上方的对应的黄色的行的 c(u) 大 1 。

可以看出,p(10, 0) = 3, p(10, 1) = 8, p(10, 2) = 9, p(10, 3) = 3 。

解答

根据以上分析,我们有以下 C# 程序:

 1: using System;
 2: using Skyiv.Utils;
 3: 
 4: static class E609
 5: {
 6:   static readonly int n = (int)1e8, m = 1000000007;
 7:   static readonly int[] a = Primes.GetPrimes(0, n);
 8: 
 9:   static bool IsPrime(int n) { return Array.BinarySearch(a, n) >= 0; }
10: 
11:   static int PrimeCount(int max)
12:   {
13:     var n = Array.BinarySearch(a, max);
14:     return (n < 0) ? ~n : (n + 1);
15:   }
16: 
17:   static void Main()
18:   {
19:     var ps = new int[64];
20:     for (var i = 0; i < a.Length; i++)
21:       for (int c = 0, d = ((i < a.Length-1) ? (a[i+1]-1) : n) - a[i],
22:         j = a[i]; (j = PrimeCount(j)) > 0; ps[c]++, ps[c + 1] += d)
23:         if (!IsPrime(j)) c++;
24:     long z = 1;
25:     foreach (var p in ps) if (p != 0) z = (z * p) % m;
26:     Console.WriteLine(z);
27:   }
28: }

简要说明:

  • 第 7 行的数组 a 存放不超过 n 的素数,Primes 类见参考资料1。
  • 第 19 行的数组 ps 存放 p(n, k) 的值,其中 n 在第 6 行定义,k 就是数组 ps 的下标。实际上,当 u0 = 108 时,u 序列有 12 个元素,其中 2 个是素数,c(u) = 10。也就是说,数组 ps 的长度设为 12 就够了。
  • 第 20 至 23 行的 for 循环遍历素数数组 a 。
  • 第 21 至 23 行的 for 循环对每个素数 a[i] 生成 π 序列。
  • 第 21 行的 d 就是每一个青色块的行数。
  • 第 22 行通过累计来计算 p(n, k) 的值。ps[c] 对应 u0 为素数的情况(黄色的行),ps[c + 1] 对应 u0 为合数的情况(青色的块)。

这个程序的运行时间是 22.1 秒。

优化

实际上,IsPrime 和 PrimeCount 方法都调用了 Array.BinarySearch 方法,这两次调用可以合并为一次。因此,我们有以下优化后的 C# 程序:

 1: using System;
 2: using Skyiv.Utils;
 3: 
 4: static class E609
 5: {
 6:   static void Main()
 7:   {
 8:     int n = (int)1e8, m = 1000000007;
 9:     int[] ps = new int[64], a = Primes.GetPrimes(0, n);
10:     for (var i = 0; i < a.Length; i++)
11:       for (int c = 0, d = ((i < a.Length-1) ? (a[i+1]-1) : n) - a[i],
12:         j = i; (j = (j < 0) ? ~j : (j+1)) > 0; ps[c]++, ps[c+1] += d)
13:         if ((j = Array.BinarySearch(a, j)) < 0) c++;
14:     long z = 1;
15:     foreach (var p in ps) if (p != 0) z = (z * p) % m;
16:     Console.WriteLine(z);
17:   }
18: }

简要说明:

  • 第 12 和 13 行的 j 就是 π 序列中的元素,也就是 PrimeCount 的结果。
  • 第 13 行的 Array.BinarySearch 方法同时起到 PrimeCount 和 IsPrime 的作用。

这个程序的运行时间是 10.6 秒。

进一步优化

由于不同的 π 序列中有许多相同的元素,我们可以缓存这些相同的元素,进一步提高速度:

在第 9 行末尾的分号前增加数组 b 用作缓存:

, b = new int[a.Length + 1]

将第 13 行替换为以下语句以利用缓存:

{
  if (b[j] == 0) b[j] = Array.BinarySearch(a, j) + n;
  if ((j = b[j] - n) < 0) c++;
}

进一步优化后的程序运行时间是 2.6 秒。

参考资料

  1. 图灵社区:使用筛法生成素数
  2. Wikipedia: Prime-counting function
  3. MSDN: Array.BinarySearch Method