问题

素数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好用一点！

bruce天地不仁，以万物为刍狗

微信支付

支付宝

PayPal链接

https://paypal.me/xieweiking