题目111：寻找含有最多重复位的10位的素数

欧拉计划

Primes with runs

Considering 4-digit primes containing repeated digits it is clear that they cannot all be the same: 1111 is divisible by 11, 2222 is divisible by 22, and so on. But there are nine 4-digit primes containing three ones:

1117, 1151, 1171, 1181, 1511, 1811, 2111, 4111, 8111

We shall say that M(n, d) represents the maximum number of repeated digits for an n-digit prime where d is the repeated digit, N(n, d) represents the number of such primes, and S(n, d) represents the sum of these primes.

So M(4, 1) = 3 is the maximum number of repeated digits for a 4-digit prime where one is the repeated digit, there are N(4, 1) = 9 such primes, and the sum of these primes is S(4, 1) = 22275. It turns out that for d = 0, it is only possible to have M(4, 0) = 2 repeated digits, but there are N(4, 0) = 13 such cases.

In the same way we obtain the following results for 4-digit primes.

For d = 0 to 9, the sum of all S(4, d) is 273700.

Find the sum of all S(10, d).

题目：

考虑含有重复位的四位数素数。显然这四位数不能全部相同：1111 可以被 11 整除，2222 可以被 22 整除，等等。但是有 9 个四位的素数包含三个 1：

1117, 1151, 1171, 1181, 1511, 1811, 2111, 4111, 8111

我们用 M(n, d) 来表示一个 n 位素数中重复位的最大值，其中 d 代表重复的数字。用 N(n, d) 代表这样的素数的数量，S(n, d) 来代表这些素数之和。

所以 M(4,1) = 3，是四位的素数中重复位的最大值，其中 1 是重复的那一位数。共有 N(4,1)=9 个这样的素数，他们的和 S(4,1)=22275。事实证明对于  d=0，只能有 M(4,0)=2 个重复位，但是却有 N(4,0)=13 个这样的数。

用同样的方法我们可以得到对于四位数素数的如下结果：

对于 d=0 到 9，S(4,d) 之和为 273700。

求所有 S(10,d) 之和。

jerryxjr1220

当n=4时，程序已经验证成功了，而且速度也挺快 0.5秒不到
但是放到n=10的时候，感觉要等好久，看来还要优化算法。
273700
Time used: 0.34

1. # -*- coding: utf-8 -*-
   
2. import itertools as it
   
3. import math
   
4. import time
   
5. start = time.time()
   
6. def getPrime(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.         plist = []
    
14.         for (k,trueprime) in enumerate(primes):
    
15.                 if trueprime:
    
16.                         plist.append(k)
    
17.         return plist
    
18. primes = getPrime(100000)
    
19. def isPrime(n):
    
20.         if n > 100000:
    
21.                 for i in primes:
    
22.                         if n%i==0:
    
23.                                 return False
    
24.                         if i*i > n:
    
25.                                 break
    
26.                 return True
    
27.         else:
    
28.                 if n in primes:
    
29.                         return True
    
30.                 else:
    
31.                         return False
    
32. 
    
33. def pr(n,d):
    
34.         pl = list(range(d))+list(range(d+1,10))
    
35.         for i in range(n-1,1,-1):
    
36.                 cl = []
    
37.                 com = it.combinations(pl,(n-i))
    
38.                 for each in com:
    
39.                         tl = list(each)
    
40.                         for j in range(i):
    
41.                                 tl.append(d)
    
42.                         cl.append(tl)
    
43.                 plst = set()
    
44.                 for ecom in cl:
    
45.                         al = it.permutations(ecom,n)
    
46.                         for a in al:
    
47.                                 s = 0
    
48.                                 for j in range(n):
    
49.                                         s = s*10+a[j]
    
50.                                 if isPrime(s) and int(math.log10(s)) == n-1:
    
51.                                         plst.add(s)
    
52.                 if len(plst) > 0:
    
53.                         break
    
54.         return plst
    
55. 
    
56. print (sum(sum(pr(4,k)) for k in range(10)))
    
57. print ('Time used: %.2f' % (time.time()-start))

复制代码

jerryxjr1220

对算法进行了优化

1. # -*- coding: utf-8 -*-
   
2. from itertools import combinations
   
3. 
   
4. import time
   
5. start = time.time()
   
6. def getPrime(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.         plist = []
    
14.         for (k,trueprime) in enumerate(primes):
    
15.                 if trueprime:
    
16.                         plist.append(k)
    
17.         return plist
    
18. primes = getPrime(100000)
    
19. def isPrime(n):
    
20.         if n > 100000:
    
21.                 for i in primes:
    
22.                         if n%i==0:
    
23.                                 return False
    
24.                         if i*i > n:
    
25.                                 break
    
26.                 return True
    
27.         else:
    
28.                 if n in primes:
    
29.                         return True
    
30.                 else:
    
31.                         return False
    
32. 
    
33. def generateDigitLists(n,without, withZero=False):
    
34.     if n == 0:
    
35.         yield []
    
36.     else:
    
37.         start = 0 if withZero else 1
    
38.         for lst in generateDigitLists(n-1,without, True):
    
39.             for i in range(start,10):
    
40.                 if i != without:
    
41.                     yield [i] + lst
    
42. 
    
43. def F(n,d):
    
44.     for k in range(n-1,0,-1):
    
45.         resultSet = set()
    
46.         start = 0 if d != 0 else 1
    
47.         for dInds in combinations(range(start,n), k):
    
48.             for dList in generateDigitLists(n-k,d, d != 0 and 0 in dInds):
    
49.                 j=0
    
50.                 lst = [d]*n
    
51.                 for i in range(n):
    
52.                     if i not in dInds:
    
53.                         lst[i] = dList[j]
    
54.                         j += 1
    
55.                 number = int(''.join(str(dig) for dig in lst))
    
56.                 if isPrime(number) :
    
57.                     resultSet.add(number)
    
58.         if len(resultSet) > 0:
    
59.             return k,resultSet
    
60. 
    
61. 
    
62. def MNS(n,d):
    
63.     k,rs = F(n,d)
    
64.     return k,len(rs),sum(rs)
    
65. 
    
66. def e111(n):
    
67.     r = 0
    
68.     for d in range(10):
    
69.         m,ln,s = MNS(n,d)
    
70.         r += s
    
71.     return r
    
72. 
    
73. print (e111(10))
    
74. 
    
75. print ("Time used: %.2f sec" % (time.time()-start))
    

复制代码

guosl

M(10,0) = 8     N(10,0) = 8     S(10,0) = 38000000042
M(10,1) = 9     N(10,1) = 11    S(10,1) = 12882626601
M(10,2) = 8     N(10,2) = 39    S(10,2) = 97447914665
M(10,3) = 9     N(10,3) = 7     S(10,3) = 23234122821
M(10,4) = 9     N(10,4) = 1     S(10,4) = 4444444447
M(10,5) = 9     N(10,5) = 1     S(10,5) = 5555555557
M(10,6) = 9     N(10,6) = 1     S(10,6) = 6666666661
M(10,7) = 9     N(10,7) = 9     S(10,7) = 59950904793
M(10,8) = 8     N(10,8) = 32    S(10,8) = 285769942206
M(10,9) = 9     N(10,9) = 8     S(10,9) = 78455389922
Sum = 612407567715
耗时：0.0164531秒
主要检验的数不能用循环来产生。而要根据规律用“填数”方法来产生。

guosl

1. #include <iostream>
   
2. #include <cstdio>
   
3. #include <cstring>
   
4. #include <cmath>
   
5. #include <omp.h>
   
6. using namespace std;
   
7. 
   
8. const int nLen = 10;//讨论的素数位数
   
9. 
   
10. char cPrimes[100104];//cPrimes[i]=1表示i是一个素数
    
11. int nPrimes[9600];//具体记录素数
    
12. int nps = 0;//记录在0到100100之间素数的个数
    
13. int nCount[10];//nCount[i]记录N(nLen,i)的值
    
14. long long lTSum[10];//lTSum[i]记录S(nLen,i)的值
    
15. int nMaxCount[10];//nMaxCount[i]记录M(Len,i)的值
    
16. 
    
17. void dfs(int nCt, //当前未变动的数字的个数，也是取指定数字的个数
    
18.   int nPos, //填充的位置
    
19.   const char* str, //逐渐被填充的字符串（最后形成欲求的素数）
    
20.   int nChkInt)//指定的数字
    
21. {
    
22.   if (nCt == nMaxCount[nChkInt])//所有可变动的指标已经用完
    
23.   {
    
24.     //检查最后生成的数字是否是一个素数
    
25.     long long x = 0;
    
26.     for (int i = 0; i < nLen; ++i)
    
27.     {
    
28.       x *= 10;
    
29.       x += str[i];
    
30.     }
    
31.     int nD = sqrt((double)x);
    
32.     bool bIsPrime = true;
    
33.     for (int i = 0; nPrimes[i] <= nD; ++i)
    
34.     {
    
35.       long long p = (long long)nPrimes[i];
    
36.       if (x%p == 0)
    
37.       {
    
38.         bIsPrime = false;
    
39.         break;
    
40.       }
    
41.     }
    
42.     if (!bIsPrime)
    
43.       return;//不是素数
    
44.     ++nCount[nChkInt];//保存结果
    
45.     lTSum[nChkInt] += x;
    
46.     return;
    
47.   }
    
48.   if (nPos >= nLen)//超出可填充的位置
    
49.     return;
    
50.   for (int i = 0; i < 10; ++i)//继续枚举在第nPos位置上的数字
    
51.   {
    
52.     char str1[16];
    
53.     memcpy(str1, str, sizeof(str1));
    
54.     str1[nPos] = i;
    
55.     if (i == nChkInt)
    
56.       dfs(nCt, nPos + 1, str1, nChkInt);//如果是指定值则未占用变动指标
    
57.     else
    
58.       dfs(nCt - 1, nPos + 1, str1, nChkInt);//如果不是指定值则占用变动指标
    
59.   }
    
60. }
    
61. 
    
62. int main(void)
    
63. {
    
64.   double t = omp_get_wtime();
    
65.   //筛法求批量素数
    
66.   memset(cPrimes + 2, 1, sizeof(cPrimes) - 2);
    
67.   for (int i = 2; i <= (int)sqrt(100100.0); ++i)
    
68.   {
    
69.     if (cPrimes[i] == 0)
    
70.       continue;
    
71.     for (int j = i*i; j < 100100; j += i)
    
72.       cPrimes[j] = 0;
    
73.   }
    
74.   //保存具体的素数
    
75.   nPrimes[nps++] = 2;
    
76.   for (int i = 3; i < 100100; i += 2)
    
77.   {
    
78.     if (cPrimes[i] == 1)
    
79.       nPrimes[nps++] = i;
    
80.   }
    
81.   long long lSum = 0;
    
82. #pragma omp parallel for reduction(+:lSum)
    
83.   for (int nChkInt = 0; nChkInt < 10; ++nChkInt)//求M(Len,nChkInt)，N(Len,nChkInt)，S(Len,nChkInt)
    
84.   {
    
85.     nMaxCount[nChkInt] = nLen;
    
86.     while (nCount[nChkInt] == 0)
    
87.     {
    
88.       --nMaxCount[nChkInt];
    
89.       char str[16];//记录最后满足条件的素数
    
90.       memset(str, 0, sizeof(str));
    
91.       for (int i = 1; i <= 9; ++i)//枚举要生成的素数的第0位上的数字
    
92.       {
    
93.         int k = nLen;
    
94.         memset(str, nChkInt, nLen * sizeof(char));//把该数所有位置上的值都填成nChkInt
    
95.         str[0] = i;
    
96.         if (i != nChkInt)
    
97.           --k;//未取指定的值，可变动数字的个数减一
    
98.         dfs(k, 1, str, nChkInt);//逐步填充后面的位置上的数字
    
99.       }
    
100.     }
     
101.     lSum += lTSum[nChkInt];
     
102.   }
     
103.   //输出结果
     
104.   for (int i = 0; i < 10; ++i)
     
105.     cout << "M(" << nLen << "," << i << ") = " << nMaxCount[i] << "\tN(" << nLen << "," << i << ") = " << nCount[i] << "\tS(" << nLen << "," << i << ") = " << lTSum[i] << endl;
     
106.   cout << "Sum = " << lSum << endl;
     
107.   t = omp_get_wtime() - t;
     
108.   cout << "耗时：" << t << "秒" << endl;
     
109.   return 0;
     
110. }

复制代码