题目

Prolblem 578: Integers with decreasing prime powers

Any positive integer can be written as a product of prime powers: p1a1 × p2a2 × ... × pkak,
where pi are distinct prime integers, ai > 0 and pi < pj if i < j.

A decreasing prime power positive integer is one for which ai ≥ aj if i < j.
For example, 1, 2, 15=3×5, 360=23×32×5 and 1000=23×53 are decreasing prime power integers.

Let C(n) be the count of decreasing prime power positive integers not exceeding n.
C(100) = 94 since all positive integers not exceeding 100 have decreasing prime powers except 18, 50, 54, 75, 90 and 98.
You are given C(106) = 922052.

Find C(1013).

题解

第578题:素数的降幂构成的整数

任何正整数都可以写成素数的幂的乘积: p1a1 × p2a2 × ... × pkak,
这里 pi 是不同的素数,对于 i < j,ai > 0 且 pi < pj

如果对于 i < j 有 ai ≥ aj,则称为素数的降幂构成的正整数。
例如, 1, 2, 15=3×5, 360=23×32×5 和 1000=23×53 是素数的降幂构成的整数。

令 C(n) 是素数的降幂构成的不超过 n 的正整数的个数。
因为在不超过 100 的所有正整数中,除了 18, 50, 54, 75, 90 和 98 以外,都是由素数的降幂构成,所以 C(100) = 94。
我们有 C(106) = 922052.

请计算 C(1013).

暴力破解

根据题意,我们有以下 C# 程序:

 1: using System;
 2: using Skyiv.Utils;
 3: 
 4: static class E578A
 5: {
 6:   static int[] primes;
 7:   
 8:   static int C(int n, int i, int e)
 9:   {
10:     int p = primes[i], z = 0;
11:     if (p > n) return 1;
12:     for (int j = 0; j <= e && n >= 1; ++j, n /= p)
13:       z += C(n, i + 1, (j != 0) ? j : e);
14:     return z;
15:   }
16: 
17:   static void Main()
18:   {
19:     int n = (int)1e6;
20:     primes = Primes.GetPrimes(0, n, true);
21:     Console.WriteLine(C(n, 0, 64));
22:   }
23: }

简要说明:

  • 第 20 行的 Primes.GetPrimes 方法返回不超过n的所有素数,且尾部附加有一个哨兵,该哨兵的值大于n,但不一定是素数。该方法请参见参考资料[1]。
  • 第 8 至 15 行 C(n, i, e) 方法,就是通过递归调用自身来计算 C(n) 的。其中,i是各个素数的下标,e是各个素数的所允许的最大指数。
  • 第 21 行给出的e的初值是 64,这是因为第一个素数是 2,而n不可能大于 264
  • 第 11 行,如果素数p大于n,说明已经递归到底了,没有新的解,直接返回 1 (表示已经找到一个解)。
  • 第 12 行,对指数j进行循环,从 0 开始,直到最大指数e。指数每增 1,n就要除以p。如果达到最大指数之前,n已经降到 1,则停止循环。
  • 第 13 行,递归调用 C 方法自身,以完成计算。此时,往右移一个素数(i + 1)。如果当前指数j是 0,说明当前乘数是 1,即缺失当前位置的素数p,则指数限制为传入的指数e(即其左边的指数),否则指数限制为j,表示其右边的指数不能超过当前位置的指数。

这个程序计算出 C(106) = 922,052,运行时间是 0.1 秒。但是,它无法计算 C(1013)。

  • 首先,n = 1013 已经超过 int.MaxValue = 231 - 1,这个程序需要使用不超过 n 的所有素数,我们的计算机不可能有这么多的内存。
  • 其次,这个程序是一个个累加计算的,对于很大的数,递归调用会导致堆栈溢出。
  • 最后,对于很大的数,即使有无限的堆栈,运行时间也太长了。

改进

经过观察,如果 p2 > n,则只可能再乘以一个素数(的 1 次方),所以有以下 C# 程序:

 1: using System;
 2: using Skyiv.Utils;
 3: 
 4: static class E578
 5: {
 6:   static int[] primes;
 7:   static PrimeHelper helper;
 8:   
 9:   static long C(long n, int i, int e)
10:   {
11:     long p = primes[i], z = 0;
12:     if (p * p > n) z = Math.Max(0, helper.π(n) - i) + 1;
13:     else for (int j = 0; j <= e && n >= 1; ++j, n /= p)
14:       z += C(n, i + 1, (j != 0) ? j : e);
15:     return z;
16:   }
17: 
18:   static void Main()
19:   {
20:     long n = (long)1e13;
21:     Console.WriteLine("n = {0:N0}", n);
22:     helper = new PrimeHelper(n);
23:     primes = helper.GetPrimes(true);
24:     Console.WriteLine("π(√n) = {0:N0}", primes.Length-1);
25:     Console.WriteLine("C(n) = {0:N0}", C(n, 0, 64));
26:   }
27: }

简要说明:

  • 第 22 行的 PrimeHelper 类请见参考资料[2]。
  • 第 23 行得到不超过 的所有素数,尾部附加有一个哨兵(该哨兵的值大于 ,但不一定是素数)。注意,这里只要求不超过 的素数,而上一小节的程序要求不超过 n 的素数。
  • 第 9 至 16 行的 C 方法基本上与上一小节的程序相同,除了第 12 行。
  • 第 12 行,当 p2 > n 时,只能乘以一个素数(的 1 次方)。所以返回 π(n) - π(p-1) + 1。
    • 此时,π(p-1) == i(注意,i 的初值是 0,而不是 1)。
    • 还要注意边界情况,当 p > n 时,π(n) - π(p-1) 有可能是负数,需要调整为 0。
    • 加 1 是因为除了素数以外,乘以 1 也算一个。

这个程序计算出 C(1013) ≈ 9.2 ×1012,运行时间是 74 秒。

运行过程分析

对于 n = 71,这个程序的运行结果如下所示:

n = 71
π(√n) = 4
C(n) = 68

运行过程(第一行蓝色格子里是i的值。p = 9 是哨兵):

上图中黄色格子的值由以下公式计算(红线标出):

  • 19 = 17+1+1, 17 = π(71)-4+1, 1 = π(4)-4+1, 1 = Max(0,π(1)-4)+1
  • 24 = 19+4+1, 4 = π(14)-3+1, 1 = Max(0,π(2)-3)+1
  • 36 = 24+8+3+1, 8 = π(23)-2+1, 3 = π(7)-2+1, 1 = Max(0,π(2)-2)+1
  • 11 = 9+2, 9 = π(35)-3+1, 2 = π(7)-3+1
  • 15 = 11+4, 4 = π(11)-2+1
  • 9 = 6+2+1, 6 = π(17)-2+1, 2 = π(5)-2+1, 1 = Max(0,π(1)-2)+1
  • 68 = 36+15+9+4+2+1+1, 4 = π(8)-1+1, 2 = π(4)-1+1, 1 = π(2)-1+1, 1 = Max(0,π(1)-1)+1

更快的程序

在欧拉计划的论坛中,俄罗斯的 aleksey 给出了下面 OCaml 程序:

I wrote n as

where ei > 1. The left part (with powers) is easy to enumerate with brute force, and the right part can be counted as

where the sum is over all numbers with prime factors less than pj+1 in the first power, and others in the second. 26 seconds in OCaml.

这个程序在 aleksey 的机器上的运行时间是 26 秒,在我的机器上的运行时间是 38 秒。

附加题

  1. 请问以下极限是否存在:
  2. 如果存在,它是代数数还是超越数?

参考资料

  1. 图灵社区:使用筛法生成素数
  2. 图灵社区:素数计数函数
  3. Stack Overflow Problems