问题
素数3、7、109和673十分特别。在它们当中取出任意2个数，前后连接在一起，都会得到一个素数。举例，取出7和109，连接起来得到7109和1097都是素数。这个集合之和是792，这也是具有这种性质的4个素数集合的和的最小值。
请找出取出任意2个数、前后连接后都能产生一个素数的5个素数集合的和的最小值。
* 传送门
分析
我靠，这个题目如果是直接暴力解题，假设搜索上界是小于N的素数，N是一个k位，那么预先生成的素数表（用于快速判断是否素数）则要以\(N \cdot 10^k\)作为上界。假设小于N的素数有x个，要在里面任意找5个素数出来进行两两拼接，则需要做\(C_x^5\)这么多次迭代，每次迭代还得再做\(C_5^2\)次判断，总共需要做\(C_x^5 \cdot C_5^2 = \frac{x!}{5!(x-5)!} \cdot \frac{5!}{2!(5-2)!} = \frac{x!}{2!3!(x-5)!} = \frac{x(x-1)(x-2)(x-3)(x-4)}{12} < \frac{x^5}{12}\)次操作，即暴力解题是一个时间复杂度达到\(O(n^5)\)的算法……更麻烦的是，不能用短路代码——遇到第一组满足条件的就当作是结果，因为题目要求的是和最小，即便我们从最小的素数开始往后遍历，即便找到满足条件的素数，也有可能不是和最小的！
我核算过，在小于1000的时候找那个只有4个素数且和最小的集合，备选素数只有168个，总共大概需要2.7秒，即我的CPU和JS能在140纳秒做完一次两个前后连接是否为素数的判断；而上升到小于10000的时候，备选素数只有1229个，但耗时将会惊人的超过1年！这已经是天文数字……
这就是欧拉计划第一个难度达到20%的题目。敬畏它吧！
更可笑的是，原本以为写好暴力方法就去解题论坛找个更优雅的方法。但我实在不想花超过一年时间就为等这一个答案，所以我便没办法去逛那个解题论坛（没有答对的人不能逛那一题的论坛）……那就只好Google一下偷看人家的答案了……
结果还是挺好笑的，我发现自己又忽略了重要的细节——重复判断其实是可以优化的！即便继续使用暴力解题的形式，还是可以在较短（长）时间内得到结果的！
解题论坛的法国佬Begoner又提出了一个十分有用的观察：
我以素数们的每位之和模3差别对待它们。这样可能有一定的优化，但也可能会使得程序有点慢，但是：
每个素数（除了3），它们的每位之和模3要么余1、要么余2。如果你将余为1的素数和余为2的素数连接起来，你肯定会得到一个每位之和能整除3的
素数。所以我创建了两个数组——一个用来存放模3余为1的素数（数组1）、另一个则存放余为2的素数（数组2）。3同时存在于两个数组中，因为不论任何其它素数和3相连接，产生新的数的每位之和肯定都不能整除3。最后，所有能满足题目要求的数要么都在数组1中、要么都在数组2中，这样子就只需要搜索更少的素数（组合）。
我有好几次遇到难题都会碰到这个法国佬的答案，实在令人佩服！
基于这个观察，我也可以将含有2和5的组合撇除，因为以2或5结尾的肯定不是素数、是合数。
答案
千万不要执行这个脚本！！！即便它是对的！！除非你肯等一辈子！！！
// answer1: brute force
const start = Date.now()
function* primes_below(n) {
const m = Math.floor((n - 2) / 2)
const sieve = new BitSet
for (let i = 1; i < m; ++i) {
let j = i, p = i + j + 2*i*j
for (; p <= m; ++j, p = i + j + 2*i*j)
sieve.set(p, 1)
}
yield 2
for (let i = 1; i <= m; ++i)
if (sieve.get(i) == 0)
yield (2*i + 1)
}
const limit = 10000
const L = (n) => Math.floor(Math.log10(n)) + 1
const scale = (n) => Math.pow(10, L(n))
const is_prime = new BitSet
const primes = []
for (const p of primes_below(limit * scale(limit))) {
is_prime.set(p, 1)
primes.push(p)
}
const concat = (n_1, n_2) => n_1 * scale(n_2) + n_2
function can_concat_as_prime(...numbers) {
const last = numbers.length - 1
const last_n = numbers[last]
for (let i = 0; i < last; ++i) {
const n = numbers[i]
if (!is_prime.get(concat(n, last_n)) || !is_prime.get(concat(last_n, n)))
return false
}
return true
}
let result = Infinity
for (let i_1 = 0, p_1 = primes[i_1]; p_1 < limit; p_1 = primes[++i_1])
for (let i_2 = i_1 + 1, p_2 = primes[i_2]; p_2 < limit; p_2 = primes[++i_2])
if (can_concat_as_prime(p_1, p_2))
for (let i_3 = i_2 + 1, p_3 = primes[i_3]; p_3 < limit; p_3 = primes[++i_3])
if (can_concat_as_prime(p_1, p_2, p_3))
for (let i_4 = i_3 + 1, p_4 = primes[i_4]; p_4 < limit; p_4 = primes[++i_4])
if (can_concat_as_prime(p_1, p_2, p_3, p_4))
for (let i_5 = i_4 + 1, p_5 = primes[i_5]; p_5 < limit; p_5 = primes[++i_5])
if (can_concat_as_prime(p_1, p_2, p_3, p_4, p_5)) {
const sum = p_1 + p_2 + p_3 + p_4 + p_5
console.log([p_1, p_2, p_3, p_4, p_5].join('+') + '=' + sum)
if (result > sum)
result = sum
}
alert(result + ', ' + (Date.now() - start) + 'ms')
这个版本的暴力解题其实至少也要花几分钟。
// answer2: brute force with cache
const start = Date.now()
function* primes_below(n) {
const m = Math.floor((n - 2) / 2)
const sieve = new BitSet
for (let i = 1; i < m; ++i) {
let j = i, p = i + j + 2*i*j
for (; p <= m; ++j, p = i + j + 2*i*j)
sieve.set(p, 1)
}
yield 2
for (let i = 1; i <= m; ++i)
if (sieve.get(i) == 0)
yield (2*i + 1)
}
const limit = 10000 // 这个上界其实是碰运气找到的，上界越大，后面满足条件的素数越多，碰巧和最小的就在一万以内
const L = (n) => Math.floor(Math.log10(n)) + 1
const scale = (n) => Math.pow(10, L(n))
const is_prime = new BitSet
const primes = []
for (const p of primes_below(limit * scale(limit))) {
is_prime.set(p, 1)
if (p < limit)
primes.push(p)
}
const concat = (n_1, n_2) => n_1 * scale(n_2) + n_2
function create_pair_mask(i) { // 一次性将后面要配对的素数都求出结果
const mask = new BitSet
const p = primes[i]
for (let j = i + 1; j < primes.length; ++j) {
const p_1 = primes[j]
if (is_prime.get(concat(p, p_1)) && is_prime.get(concat(p_1, p)))
mask.set(p_1, 1)
}
return mask
}
const mask_cache = new Array(primes.length)
function get_pair_mask(i) { // 利用缓存避免重复比较
let pair_mask = mask_cache[i]
if (pair_mask === undefined)
pair_mask = mask_cache[i] = create_pair_mask(i)
return pair_mask
}
let result = Infinity
for (let i_1 = 1; i_1 < primes.length - 4; ++i_1) { // 2和任何素数连接肯定都不会满足要求，所以从1开始
const p_1 = primes[i_1]
if (p_1*5 >= result) // 提前判断是否会超出之前找到的和
break
const mask_1 = get_pair_mask(i_1)
for (let i_2 = i_1 + 1; i_2 < primes.length - 3; ++i_2) {
const p_2 = primes[i_2]
if (p_1 + p_2*4 >= result)
break
if (mask_1.get(p_2)) {
const mask_2 = get_pair_mask(i_2)
for (let i_3 = i_2 + 1; i_3 < primes.length - 2; ++i_3) {
const p_3 = primes[i_3]
if (p_1 + p_2 + p_3*3 >= result)
break
if (mask_1.get(p_3) && mask_2.get(p_3)) {
const mask_3 = get_pair_mask(i_3)
for (let i_4 = i_3 + 1; i_4 < primes.length - 1; ++i_4) {
const p_4 = primes[i_4]
if (p_1 + p_2 + p_3 + p_4*2 >= result)
break
if (mask_1.get(p_4) && mask_2.get(p_4) && mask_3.get(p_4)) {
const mask_4 = get_pair_mask(i_4)
for (let i_5 = i_4 + 1; i_5 < primes.length; ++i_5) {
const p_5 = primes[i_5]
if (p_1 + p_2 + p_3 + p_4 + p_5 >= result)
break
if (mask_1.get(p_5) && mask_2.get(p_5) && mask_3.get(p_5) && mask_4.get(p_5)) {
const sum = p_1 + p_2 + p_3 + p_4 + p_5
console.log([p_1, p_2, p_3, p_4, p_5].join('+') + '=' + sum)
result = sum
}
}
}
}
}
}
}
}
}
alert(result + ', ' + (Date.now() - start) + 'ms')
这个方法除了使用了模3进行测试以外，还使用了集合交集的方式减少代码，但运行起来速度并没有什么明显的提升！
// answer3: brute force with mod 3
const start = Date.now()
Set.prototype.intersection = function(setB) {
const intersection = new Set()
for (const elem of setB)
if (this.has(elem))
intersection.add(elem)
return intersection
}
Array.prototype.search = function(el, min_idx = 0, max_idx = 0) {
max_idx = max_idx || this.length - 1
while (min_idx <= max_idx) {
const curr_idx = (min_idx + max_idx) >> 1
const curr_el = this[curr_idx]
if (curr_el < el)
min_idx = curr_idx + 1
else if (curr_el > el)
max_idx = curr_idx -1
else
return curr_idx
}
return -1
}
function* primes_below(n) {
const m = Math.floor((n - 2) / 2)
const sieve = new BitSet
for (let i = 1; i < m; ++i) {
let j = i, p = i + j + 2*i*j
for (; p <= m; ++j, p = i + j + 2*i*j)
sieve.set(p, 1)
}
yield 2
for (let i = 1; i <= m; ++i)
if (sieve.get(i) == 0)
yield (2*i + 1)
}
const limit = 10000
const L = (n) => Math.floor(Math.log10(n)) + 1
const scale = (n) => Math.pow(10, L(n))
const is_prime = new BitSet
const primes_1 = [3], primes_2 = [3]
for (const p of primes_below(limit * scale(limit))) {
is_prime.set(p, 1)
if (5 < p && p < limit) {
const r = p % 3 // == digit_sum(p) % 3
if (r == 1)
primes_1.push(p)
else if (r == 2)
primes_2.push(p)
}
}
const concat = (n_1, n_2) => n_1 * scale(n_2) + n_2
function create_pairs(primes, i) {
const pairs = new Set()
const p = primes[i]
for (let j = i + 1; j < primes.length; ++j) {
const p_1 = primes[j]
if (is_prime.get(concat(p, p_1)) && is_prime.get(concat(p_1, p)))
pairs.add(p_1)
}
return pairs
}
const pairs_cache = {}
function get_pairs(primes, p) {
if (p in pairs_cache)
return pairs_cache[p]
else
return pairs_cache[p] = create_pairs(primes, primes.search(p))
}
let result = Infinity
function search_combinations(primes) {
for (const p_1 of primes) {
const sum_1 = p_1
if (sum_1 + p_1*4 >= result)
break
const pairs_1 = get_pairs(primes, p_1)
for (const p_2 of pairs_1) {
const sum_2 = sum_1 + p_2
if (sum_2 + p_2*3 >= result)
break
const pairs_2 = get_pairs(primes, p_2).intersection(pairs_1)
for (const p_3 of pairs_2) {
const sum_3 = sum_2 + p_3
if (sum_3 + p_3*2 >= result)
break
const pairs_3 = get_pairs(primes, p_3).intersection(pairs_2)
for (const p_4 of pairs_3) {
const sum_4 = sum_3 + p_4
if (sum_4 + p_4 >= result)
break
const pairs_4 = get_pairs(primes, p_4).intersection(pairs_3)
for (const p_5 of pairs_4) {
const sum_5 = sum_4 + p_5
if (sum_5 >= result)
break
console.log([p_1, p_2, p_3, p_4, p_5].join('+') + '=' + sum_5)
result = sum_5
}
}
}
}
}
}
search_combinations(primes_1)
if (result === Infinity) // 这里其实是想偷步知道结果
search_combinations(primes_2)
alert(result + ', ' + (Date.now() - start) + 'ms')
我感觉这类问题的确是JS不好处理的，于是写了个Java版，几乎同样逻辑的代码，几秒钟之内就可以执行完了：
/* brute force with mod 3 */
import java.util.BitSet;
import java.util.Iterator;
import java.util.List;
import java.util.ArrayList;
import java.util.Map;
import java.util.HashMap;
import java.util.Collections;
import java.math.BigInteger;
public class Problem60 {
final BitSet isPrime;
final Map<Integer, BitSet> cache;
int result = Integer.MAX_VALUE;
Problem60() {
final int limit = 20000;
final int count = Primes.estimatedCountBelow(limit);
isPrime = new BitSet(count);
cache = new HashMap<>();
final int size = count/2;
final List<Integer> primes1 = new ArrayList<>(size),
primes2 = new ArrayList<>(size);
primes1.add(3);
primes2.add(3);
for (final int prime : Primes.below(limit)) {
isPrime.set(prime);
if (5 < prime) {
final int r = prime % 3;
if (r == 1)
primes1.add(prime);
else
primes2.add(prime);
}
}
searchCombinations(primes1);
searchCombinations(primes2);
}
boolean isPrime(final int n) {
// 这里借用 JDK 自带的 Miller-Rabin 素数判断法
return isPrime.get(n) || BigInteger.valueOf(n).isProbablePrime(10);
}
int L(final int n) {
return ((int) Math.log10(n)) + 1;
}
int scale(final int n) {
return (int) Math.pow(10, L(n));
}
int concat(final int n1, final int n2) {
return n1 * scale(n2) + n2;
}
BitSet createPairs(final List<Integer> primes, final int i) {
final BitSet pairs = new BitSet();
final int prime = primes.get(i);
final int size = primes.size();
for (int j = i+1; j < size; ++j) {
final int prime1 = primes.get(j);
if (isPrime(concat(prime, prime1)) && isPrime(concat(prime1, prime)))
pairs.set(prime1);
}
return pairs;
}
BitSet getPairs(final List<Integer> primes, final Integer prime) {
return cache.computeIfAbsent(prime,
(p) -> createPairs(primes, Collections.binarySearch(primes, p)));
}
void searchCombinations(final List<Integer> primes) {
for (final int p1 : primes) {
if (p1*5 >= result)
break;
final BitSet pairs1 = getPairs(primes, p1);
for (final int p2 : new BitSetItr(pairs1, p1)) {
final int sum2 = p1 + p2;
if (sum2 + p2*3 >= result)
break;
// JDK 的 BitSet.and 方法会将结果存至左侧容器，所以是一个非线程安全、可修改容器的做法；
// 所以这里先拷贝再操作，下同。
final BitSet pairs2 = (BitSet) pairs1.clone();
pairs2.and(getPairs(primes, p2));
for (final int p3 : new BitSetItr(pairs2, p2)) {
final int sum3 = sum2 + p3;
if (sum3 + p3*2 >= result)
break;
final BitSet pairs3 = (BitSet) pairs2.clone();
pairs3.and(getPairs(primes, p3));
for (final int p4 : new BitSetItr(pairs3, p3)) {
final int sum4 = sum3 + p4;
if (sum4 + p4 >= result)
break;
final BitSet pairs4 = (BitSet) pairs3.clone();
pairs4.and(getPairs(primes, p4));
for (final int p5 : new BitSetItr(pairs4, p4)) {
final int sum5 = sum4 + p5;
if (sum5 >= result)
break;
System.out.println(p1 +"+"+ p2 +"+"+ p3 +"+"+ p4 +"+"+ p5 +"="+ sum5);
result = sum5;
}
}
}
}
}
}
public static void main(final String ... args) {
final long start = System.currentTimeMillis();
System.out.println(new Problem60().result + ", " + (System.currentTimeMillis() - start) + "ms");
}
}
class Primes implements Iterable<Integer> {
final int bits;
final BitSet sieve;
Primes(final int n) {
final int m = (n-2)/2;
bits = m+1;
sieve = new BitSet(bits);
for (int i = 1; i < m; ++i)
for (int j = i, c = oddComposite(i, j); c < bits; ++j, c = oddComposite(i, j))
sieve.set(c);
}
int oddComposite(final int i, final int j) {
return i + j + 2*i*j;
}
public Iterator<Integer> iterator() {
return new Iterator<>() {
int i = 0, next = 2;
public boolean hasNext() {
return 0 <= i && i < bits;
}
public Integer next() {
final int prime = next;
i = sieve.nextClearBit(i+1);
next = 2*i + 1;
return prime;
}
};
}
static Primes below(final int n) {
return new Primes(n);
}
static int estimatedCountBelow(final int n) {
// 利用素数定理辅助预先开辟内存空间
return (int) Math.ceil(n < 55 ? 1.25506*n/Math.log(n) : n/(Math.log(n)-4));
}
}
class BitSetItr implements Iterable<Integer> {
final BitSet bs;
final int start;
BitSetItr(final BitSet bs) {
this(bs, 0);
}
BitSetItr(final BitSet bs, final int start) {
this.bs = bs;
this.start = start;
}
public Iterator<Integer> iterator() {
return new Iterator<>() {
int i = bs.nextSetBit(start);
public boolean hasNext() {
return i >= 0;
}
public Integer next() {
final int n = i;
i = bs.nextSetBit(i+1);
return n;
}
};
}
}
对于复杂问题，还是Java好用一点！
打赏作者
