题目60：找出一个5个质数的集合，该集合中每两个质数相连接都能产生一个素数

欧拉计划

Prime pair sets

The primes 3, 7, 109, and 673, are quite remarkable. By taking any two primes and concatenating them in any order the result will always be prime. For example, taking 7 and 109, both 7109 and 1097 are prime. The sum of these four primes, 792, represents the lowest sum for a set of four primes with this property.

Find the lowest sum for a set of five primes for which any two primes concatenate to produce another prime.

题目：

质数 3, 7, 109, 和 673 是值得注意的。将其中任意两个质数以任何顺序相连接产生的结果都是质数。例如，取 7 和 109，连接而成的 7109 和 1097 都是质数。这四个质数的和是 792，这也是满足这个性质的四个质数集合的最小总和。

找出满足这个性质的五个质数的集合中，集合数之和最小的。算出这个最小的和。

迷雾少年

Mark

jerryxjr1220

Result = 26033, in 36.276000 s

1. #!/usr/bin/env python
   
2. """
   
3. q060
   
4. Algorithm (by iwek7, from forum):
   
5. 1. Generate prime list up to N. Remove 2 and 5 from list since they are trivial.
   
6.    Call the list prime_list.
   
7. 2. Create an empty list all_list, to contain lists of prime pairs.
   
8. 3. For each prime, ci, in prime_list, do
   
9.    a. create an empty new_list to contain newly found lists of prime pairs
   
10.    b. for each list l in all_list, do
    
11.       i. Test if ALL primes in l form prime pair with ci.
    
12.       ii. If it does,
    
13.           1. Create a copy of list l, call it cl
    
14.           2. Append ci to this cl
    
15.           3. If length of cl >= LEN, solution found.
    
16.           4. Append cl to new_list
    
17.    c. Append all lists in new_list to all_list
    
18.    d. Append [ci] to all_list
    
19. Optimisation:
    
20. 1. Create a is_prime_table to store the flag if a given prime pairs are prime or not. This
    
21.    is to avoid primality test for a given prime pair.
    
22. 2. Use a faster prime generation algorithm.
    
23. 3. Use a faster primality test algorithm, rather than generating list of primes and test
    
24.    if a prime pair is in the list or not.
    
25. """
    
26. import time
    
27. import numpy as np
    
28. def getPrime(n):
    
29.         primes = [True]*n
    
30.         primes[0], primes[1] = False, False
    
31.         for (i,prime) in enumerate(primes):
    
32.                 if prime:
    
33.                         for j in range(i*i,n,i):
    
34.                                 primes[j] = False
    
35.         return primes
    
36. 
    
37. def is_prime(n):
    
38.         factor = 3
    
39.         if n == 2:
    
40.                 return True
    
41.         if (n%2 == 0):
    
42.                 return False
    
43.         while factor * factor <= n:
    
44.                 if n % factor == 0:
    
45.                         return False
    
46.                 factor += 2
    
47.         return True
    
48. 
    
49. def test_cat(prime_set, ci, i, lut):
    
50.     for j, cj in prime_set:
    
51.         if lut[i][j] == 0:
    
52.             cat1 = int(str(ci) + str(cj))
    
53.             cat2 = int(str(cj) + str(ci))
    
54.             isprime = is_prime(cat1) and is_prime(cat2)
    
55.             lut[i][j] = 1 if isprime else -1
    
56.             lut[j][i] = 1 if isprime else -1
    
57.         if lut[i][j] == -1:
    
58.             return False
    
59.     return True
    
60. 
    
61. def solve_q060(prime, LEN):
    
62.     all_list = list()
    
63.     for i, ci in enumerate(prime):
    
64.         new_list = list()
    
65.         for l in all_list:
    
66.             if test_cat(l, ci, i, is_prime_table):
    
67.                 new_cand = list(l)
    
68.                 new_cand.append((i, ci))
    
69.                 if len(new_cand) >= LEN:
    
70.                     return new_cand
    
71.                 new_list.append(new_cand)
    
72.         all_list.extend(new_list)
    
73.         all_list.append([(i, ci)])
    
74.     return None
    
75. 
    
76. primes = getPrime(100000)
    
77. UPPER = 10000
    
78. LEN = 5
    
79. tic = time.time()
    
80. prime_cand = []
    
81. pl = getPrime(UPPER)
    
82. for (m,trueprime) in enumerate(pl):
    
83.         if trueprime:
    
84.                 prime_cand.append(m)
    
85. prime_cand.remove(2)
    
86. prime_cand.remove(5)
    
87. prime_len = len(prime_cand)
    
88. res_list = []
    
89. is_prime_table = np.zeros((prime_len, prime_len))
    
90. res_list = solve_q060(prime_cand, LEN)
    
91. res_list = [c[1] for c in res_list]
    
92. toc = time.time()
    
93. if res_list:
    
94.     res = sum(res_list)
    
95.     print "Result = %d, in %f s" % (res, toc - tic)
    
96. else:
    
97.     print "Not found, try increase prime numbers."

复制代码

jerryxjr1220

递归算法，20秒不到可以出结果：
13
5197
5701
6733
8389
[Finished in 19.3s]

1. def is_prime(n):
   
2.     factor = 2
   
3.     while factor * factor <= n:
   
4.         if n % factor == 0:
   
5.             return False
   
6.         factor += 1
   
7.     return True
   
8. 
   
9. def is_valid_pair(n, m):
   
10.     return is_prime(int(str(n) + str(m))) and is_prime(int(str(m) + str(n)))
    
11. 
    
12. def check_comb(nb_chosen, indices):
    
13.     if nb_chosen == NB_TO_CHOOSE:
    
14.         for i in indices:
    
15.             print(prime_list[i])
    
16.         exit()
    
17.     start = -1
    
18.     if nb_chosen > 0:
    
19.         start = indices[nb_chosen - 1]
    
20.     for i in range(start + 1, nb_primes):
    
21.         valid = True
    
22.         for j in range(nb_chosen):
    
23.             if not is_valid_pair(prime_list[i], prime_list[indices[j]]):
    
24.                 valid = False
    
25.                 break
    
26.         if valid:
    
27.             indices[nb_chosen] = i
    
28.             check_comb(nb_chosen + 1, indices)
    
29. 
    
30. UPPER_BOUND = 10 ** 4
    
31. NB_TO_CHOOSE = 5
    
32. prime_list = []
    
33. prime = [True] * UPPER_BOUND
    
34. 
    
35. for i in range(2, UPPER_BOUND):
    
36.     if prime[i]:
    
37.         prime_list.append(i)
    
38.         for j in range(i + i, UPPER_BOUND, i):
    
39.             prime[j] = False
    
40. 
    
41. nb_primes = len(prime_list)
    
42. indices = [-1] * NB_TO_CHOOSE
    
43. check_comb(0, indices)

复制代码

jerryxjr1220

另外还有一种字典法，不过用时也不短，关键生成字典的时间比较长。
方法多多

芒果加黄桃

1. # encoding:utf-8
   
2. # 寻找特殊质数  3,7,109,673已任何形式组合比如37  73  7109  1097 都是质数
   
3. # 寻找满足条件的  和最小的五个质数，满足该条件
   
4. from time import time
   
5. from math import sqrt
   
6. def getPrimes(N):
   
7.     primes = [True] * N
   
8.     primes[0], primes[1] = False, False
   
9.     for i, prime in enumerate(primes):
   
10.         if prime:
    
11.             for j in range(i * i, N, i):
    
12.                 primes[j] = False
    
13.     return [k for k, p in enumerate(primes) if p ]
    
14. def isPrime(n, l_pr):
    
15.     # 不考虑1
    
16.     if n < 4:
    
17.         return True
    
18.     if n % 2 == 0 or n % 3 == 0 :
    
19.         return False
    
20.     t = int(sqrt(n))
    
21.     for i in l_pr:
    
22.         if i > t:
    
23.             return True
    
24.         if n % i == 0:
    
25.             return False
    
26.     return True
    
27. def findNums(l_pr, l_primes, N, l_result):
    
28.     if N == 1 and len(l_result) < 1:
    
29.         l_result.append(min(l_primes))
    
30.     else:
    
31.         if len(l_result) < N:
    
32.             findNums(l_pr, l_primes, N - 1, l_result)
    
33.             for i in [x for x in l_primes if x > max(l_result)]:
    
34.                 flag = True
    
35.                 for j in l_result:
    
36.                     if not (isPrime(int(str(j) + str(i)), l_pr)  
    
37.                             and isPrime(int(str(i) + str(j)), l_pr)):
    
38.                         flag = False
    
39.                         break
    
40.                 if flag:
    
41.                     l_result.append(i)
    
42.                     return         
    
43. def euler060(N=10000):
    
44.     l_primes = getPrimes(N)
    
45.     l_result = [3]
    
46.     l_pr = [x for x in l_primes]
    
47.     findNums(l_pr, l_primes, 5, l_result)
    
48.     while len(l_result) < 5:
    
49.         l_primes = [i for i in l_pr if i > max(l_result)]
    
50.         l_result.remove(max(l_result))
    
51.         findNums(l_pr, l_primes, 5, l_result)
    
52.     print(l_result)
    
53. if __name__ == '__main__':
    
54.     start = time()
    
55.     euler060()
    
56.     print('cost %.6f sec' % (time() - start))

复制代码

[13, 5197, 5701, 6733, 8389]
cost 11.987199 sec

najin

用的matlab
结果是
          13        5197        5701        6733        8389

时间已过 3.885805 秒。
>>


第二组的话
           7        1237        2341       12409       18433
          13        5197        5701        6733        8389

时间已过 16.375198 秒。
>>

583164028

是谁说 mathematica代码好理解的   下面这段代码，能理解算我输  答案26033

1. pp[a_, b_] := FromDigits[IntegerDigits@a~Join~IntegerDigits@b]
   
2. Total[FindClique[
   
3.     Graph[#[[1]] <-> #[[2]] & /@
   
4.       Select[Prime@Range@PrimePi@10000~Subsets~{2},
   
5.        PrimeQ[#[[1]]~pp~#[[2]]] && PrimeQ[#[[2]]~pp~#[[1]]] &]]][[
   
6.    1]]]

复制代码

4444567

感觉很麻烦，但还是出来的不算太慢

1. def prime(n):
   
2.     if n==2 or n==3:
   
3.         return True
   
4.     else:
   
5.         if n%2==0:
   
6.             return False
   
7.         else:
   
8.             for i in range(3,int(n**0.5)+1,2):
   
9.                 if n%i == 0:
   
10.                     return False
    
11.             return True
    
12. 
    
13. def lj(n,m):#连接两个数字
    
14.     return prime(int(str(n)+str(m))) and prime(int(str(m)+str(n)))
    
15. 
    
16. l1 = list(i for i in range(2,10000) if prime(i))
    
17. count = 0
    
18. for i in l1:
    
19.     for j in l1:
    
20.         if lj(i,j):
    
21.             for k in l1:
    
22.                 if lj(i,k) and lj(j,k):
    
23.                     for l in l1:
    
24.                         if lj(i,l) and lj(j,l) and lj(k,l):
    
25.                             for m in l1:
    
26.                                 if lj(i,m) and lj(j,m) and lj(k,m) and lj(l,m):
    
27.                                     count += 1
    
28.                                     print(i,j,k,l,m)
    
29.                                 if count==1:
    
30.                                     exit()

复制代码

guosl

1. /*
   
2. 答案：26033
   
3. 耗时：30.2571秒（单核）
   
4.       5.08559秒（六核）
   
5. */
   
6. #include <iostream>
   
7. #include <vector>
   
8. #include <cstring>
   
9. #include <cstdlib>
   
10. #include <cmath>
    
11. #include <omp.h>
    
12. using namespace std;
    
13. 
    
14. const int N = 30000;//搜索范围一定要大于那个最小和才行，否则从理论上来讲不能得到全局最小和的
    
15. const int M = 1000000;
    
16. const int INF = 0x7fffffff;
    
17. 
    
18. char cp[M + 4];
    
19. vector<int> vp;
    
20. 
    
21. bool chkbp(long long n)//判别数n是否是一个素数
    
22. {
    
23.   if ((n & 1) == 0)
    
24.     return false;
    
25.   int d = (int)sqrt((double)n);
    
26.   for (int i = 1; i < (int)vp.size() && vp[i] <= d; ++i)
    
27.   {
    
28.     if (n%vp[i] == 0)
    
29.       return false;
    
30.   }
    
31.   return true;
    
32. }
    
33. 
    
34. bool chk(int p1, int p2)//检查p1p2与p2p1是否是素数
    
35. {
    
36.   char str1[8], str2[8];
    
37.   _itoa_s(p1, str1, 10);
    
38.   _itoa_s(p2, str2, 10);
    
39.   char str[12];
    
40.   strcpy_s(str, str1);
    
41.   strcat_s(str, str2);
    
42.   long long n = atoll(str);
    
43.   if (n < M)
    
44.   {
    
45.     if (cp[n] == 0)
    
46.       return false;
    
47.   }
    
48.   else
    
49.   {
    
50.     if (!chkbp(n))
    
51.       return false;
    
52.   }
    
53.   strcpy_s(str, str2);
    
54.   strcat_s(str, str1);
    
55.   n = atoll(str);
    
56.   if (n < M)
    
57.   {
    
58.     if (cp[n] == 0)
    
59.       return false;
    
60.   }
    
61.   return chkbp(n);
    
62. }
    
63. 
    
64. int main(void)
    
65. {
    
66.   double t = omp_get_wtime();
    
67.   memset(cp + 2, 1, (M + 2) * sizeof(char));
    
68.   for (int i = 2; i <= (int)sqrt((double)M); ++i)
    
69.   {
    
70.     if (cp[i] == 0)
    
71.       continue;
    
72.     for (int j = i * i; j <= M; j += i)
    
73.       cp[j] = 0;
    
74.   }
    
75.   vp.push_back(2);
    
76.   for (int i = 3; i <= N; i += 2)
    
77.   {
    
78.     if (cp[i] == 1)
    
79.       vp.push_back(i);
    
80.   }
    
81.   volatile int nMinSum = INF, k = 1;
    
82.   volatile bool bContinue = true;
    
83. #pragma omp parallel shared(k,bContinue,vp,nMinSum)
    
84.   while (bContinue)
    
85.   {
    
86.     if (k >= nMinSum || k >= (int)vp.size())
    
87.     {
    
88.       bContinue = false;
    
89. #pragma omp flush(bContinue)
    
90.       break;
    
91.     }
    
92.     int i0;
    
93. #pragma omp critical(c1)
    
94.     {
    
95.       i0 = k;
    
96.       ++k;
    
97.     }
    
98.     for (int i1 = i0 + 1; i1 < (int)vp.size(); ++i1)
    
99.     {
    
100.       if (vp[i0] + vp[i1] >= nMinSum)
     
101.         break;
     
102.       if (!chk(vp[i0], vp[i1]))
     
103.         continue;
     
104.       for (int i2 = i1 + 1; bContinue && i2 < (int)vp.size(); ++i2)
     
105.       {
     
106.         if (vp[i0] + vp[i1] + vp[i2] >= nMinSum)
     
107.           break;
     
108.         if (!(chk(vp[i0], vp[i2]) && chk(vp[i1], vp[i2])))
     
109.           continue;
     
110.         for (int i3 = i2 + 1; bContinue && i3 < (int)vp.size(); ++i3)
     
111.         {
     
112.           if (vp[i0] + vp[i1] + vp[i2] + vp[i3] >= nMinSum)
     
113.             break;
     
114.           if (!(chk(vp[i0], vp[i3]) && chk(vp[i1], vp[i3]) && chk(vp[i2], vp[i3])))
     
115.             continue;
     
116.           for (int i4 = i3 + 1; bContinue && i4 < (int)vp.size(); ++i4)
     
117.           {
     
118.             if (vp[i0] + vp[i1] + vp[i2] + vp[i3] + vp[i4] >= nMinSum)
     
119.               break;
     
120.             if (!(chk(vp[i0], vp[i4]) && chk(vp[i1], vp[i4]) && chk(vp[i2], vp[i4]) && chk(vp[i3], vp[i4])))
     
121.               continue;
     
122.             int nS = vp[i0] + vp[i1] + vp[i2] + vp[i3] + vp[i4];
     
123.             if (nMinSum > nS)
     
124. #pragma omp critical(c2)
     
125.               nMinSum = nS;
     
126.           }
     
127.         }
     
128.       }
     
129.     }
     
130.   }
     
131.   t = omp_get_wtime() - t;
     
132.   cout << nMinSum << endl << t << endl;
     
133.   return 0;
     
134. }

复制代码

B1tetheDust

1. import time as t
   
2. import math
   
3. 
   
4. start = t.perf_counter()
   
5. 
   
6. 
   
7. def check_prime(a_num):
   
8.     is_prime = True
   
9.     for prime_i in range(2, int(math.sqrt(a_num) + 1)):
   
10.         if a_num % prime_i == 0:
    
11.             is_prime = False
    
12.     return is_prime
    
13. 
    
14. 
    
15. def is_exchangeable(prime_list):
    
16.     length_list = len(prime_list)
    
17.     for each_prime1 in range(length_list - 1):
    
18.         for each_prime2 in range(each_prime1 + 1, length_list):
    
19.             new_num_1 = int(str(prime_list[each_prime1]) + str(prime_list[each_prime2]))
    
20.             new_num_2 = int(str(prime_list[each_prime2]) + str(prime_list[each_prime1]))
    
21.             if not check_prime(new_num_1) or not check_prime(new_num_2):
    
22.                 return False
    
23.     return True
    
24. 
    
25. 
    
26. def p60():
    
27.     primes = []
    
28.     prime_range = 10000
    
29.     for numbers in range(2, prime_range):
    
30.         if check_prime(numbers):
    
31.             primes.append(numbers)
    
32.     primes.remove(2)
    
33.     primes.remove(5)
    
34.     test_range = len(primes)
    
35.     for prime1 in range(test_range):
    
36.         for prime2 in range(prime1 + 1, test_range):
    
37.             cur_plist = [primes[prime1], primes[prime2]]
    
38.             if is_exchangeable(cur_plist):
    
39.                 for prime3 in range(prime2 + 1, test_range):
    
40.                     cur_plist = [primes[prime1], primes[prime2], primes[prime3]]
    
41.                     if is_exchangeable(cur_plist):
    
42.                         for prime4 in range(prime3 + 1, test_range):
    
43.                             cur_plist = [primes[prime1], primes[prime2], primes[prime3], primes[prime4]]
    
44.                             if is_exchangeable(cur_plist):
    
45.                                 for prime5 in range(prime4 + 1, test_range):
    
46.                                     cur_plist = [primes[prime1], primes[prime2], primes[prime3],
    
47.                                                  primes[prime4], primes[prime5]]
    
48.                                     if is_exchangeable(cur_plist):
    
49.                                         return cur_plist
    
50. 
    
51. 
    
52. answer = p60()
    
53. print(answer)
    
54. print(sum(answer))
    
55. print("It costs %f s" % (t.perf_counter() - start))

复制代码

[13, 5197, 5701, 6733, 8389]
26033
It costs 51.725813 s