【算法】验证哥德巴赫猜想

问题来源

Timus Online Judge 网站上有这么一道题目：1356. Something Easier。这道题目的输入是一组  2 到 10⁹ 之间整数，对于每个输入的整数，要求用最少个数的素数的和来表示。这道题目的时间限制是 1 秒。

问题解答

我们知道著名的哥德巴赫猜想是：

任何一个充分大的偶数都可以表示为两个素数之和

于是我们有以下的 C 语言程序(1356.c)：

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110

// http://acm.timus.ru/problem.aspx?space=1&num=1356 #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h>   // http://en.wikipedia.org/wiki/Prime_number_theorem #define PRIME_MAX 10000 #define PRIME_COUNT 1229   typedef unsigned long long U8; typedef char bool;   const bool true = 1; const bool false = 0;   // http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes bool* getSieve(int max) {   static bool sieve[(PRIME_MAX >> 1) + 1];   int i, j, imax = sqrt(max);   for (i = 3; i <= imax; i += 2)     if (!sieve[i >> 1])       for (j = i * i; j <= max; j += i << 1) sieve[j >> 1] = true;   return sieve; }   int* getPrimes(int max) {   static int primes[PRIME_COUNT + 1];   bool *sieve = getSieve(max);   int i, j = 0;   for (primes[j++] = 2, i = 3; i <= max; i += 2)     if (!sieve[i >> 1]) primes[j++] = i;   return primes; }   U8 modMultiply(U8 a, U8 b, U8 m) {   return a * b % m; }   U8 modPow(U8 a, U8 b, U8 m) {   U8 v = 1, p;   for (p = a % m; b > 0; b >>= 1, p = modMultiply(p, p, m))     if (b & 1) v = modMultiply(v, p, m);   return v; }   bool witness(U8 a, U8 n) {   U8 n1 = n - 1, s2 = n1 & -n1, x = modPow(a, n1 / s2, n);   if (x == 1 || x == n1) return false;   for (; s2 > 1; s2 >>= 1)   {     x = modMultiply(x, x, n);     if (x == 1) return true;     if (x == n1) return false;   }   return true; }   U8 random(U8 high) {   // http://www.cppreference.com/wiki/c/other/rand   return (U8)(high * (rand() / (double)RAND_MAX)); }   // http://en.wikipedia.org/wiki/Miller-Rabin_primality_test // n, an integer to be tested for primality // k, a parameter that determines the accuracy of the test bool probablyPrime(U8 n, int k) {   if (n == 2 || n == 3) return 1;   if (n < 2 || n % 2 == 0) return 0;   while (k-- > 0) if (witness(random(n - 3) + 2, n)) return false;   return true; }   bool isPrime(int n) {   return probablyPrime(n, 2); }   int outEven(int primes[], int n) {   int i, p, q;   for (i = 0; (p = primes[i]) != 0; i++)     if (isPrime(q = n - p))       return printf("%d %d", p, q);   return printf("error:%d", n); }   int main(void) {   int t, n, *primes = getPrimes(PRIME_MAX);   srand(time(NULL));   scanf("%d", &t);   while (t-- > 0)   {     scanf("%d", &n);     if (isPrime(n)) printf("%d", n);     else if ((n & 1) == 0) outEven(primes, n);     else if (isPrime(n - 2)) printf("2 %d", n - 2);     else printf("3 "), outEven(primes, n - 3);     puts("");   }   return 0; }

上述程序分析如下：

• 根据哥德巴赫猜想，充分大的偶数 n = p + q，这里 p <= q 是素数。我们猜测当 n <= 10⁹ 时，p < 10⁴。第 8 行就是定义 p 的最大值。
• 根据素数定理，我们知道 10⁴ 以内的素数有 1229 个。第 9 行就是定义程序中要用到的素数的个数。
• 第 17 到 26 行的 getSieve 函数用埃拉托斯特尼筛法筛选出素数。
• 第 28 到 36 行的 getPrimes 函数从筛中取出这些素数。
• 第 38 到 79 行的一系列函数最终是为了 probablyPrime 函数，用于检测素数。请参见我在2010年7月写的随笔：【算法】米勒-拉宾素性检验。
• 第 81 到 84 行的 isPrime 函数调用 probablyPrime 函数来检测素数。
• 第 86 到 93 行的 outEven 函数对大于 2 的偶数验证哥德巴赫猜想，即输出一对素数 p 和 q。
• 第 95 到 110 行是 main 函数。其中：
• 第 103 行处理 n 是素数的情况，直接输出该素数(包括素数 2，所以 outEven 函数处理的偶数肯定大于 2)。
• 第 104 行对大于 2 的偶数输出一对素数(通过调用 outEven 函数，强哥德巴赫猜想)。
• 第 105 行处理大于 5 的奇数能够分解为 2 和另外一个素数的和的情况(注意不要遗漏这个情形!)。
• 第 106 行处理大于 5 的奇数的其他情况，首先输出一个 3，然后调用 outEven 函数处理偶数 n - 3 (弱哥德巴赫猜想)。

上述程序在 Timus Online Judge 网站的运行时间是 0.015 秒。

参考资料

1. Wikipedia: Goldbach's conjecture
2. Wikipedia: Prime number theorem
3. Wikipedia: Sieve of Eratosthenes
4. Wikipedia: Miller–Rabin primality test
5. 【算法】米勒-拉宾素性检验

---

更多的 ACM 题的解法请参见：Timus 目录。