题目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

# -*- coding: utf-8 -*-
import itertools as it
import math
import time
start = time.time()
def getPrime(n):
        primes = [True]*n
        primes[0], primes[1] = False, False
        for (i,prime) in enumerate(primes):
                if prime:
                        for j in range(i*i,n,i):
                                primes[j] = False
        plist = []
        for (k,trueprime) in enumerate(primes):
                if trueprime:
                        plist.append(k)
        return plist
primes = getPrime(100000)
def isPrime(n):
        if n > 100000:
                for i in primes:
                        if n%i==0:
                                return False
                        if i*i > n:
                                break
                return True
        else:
                if n in primes:
                        return True
                else:
                        return False

def pr(n,d):
        pl = list(range(d))+list(range(d+1,10))
        for i in range(n-1,1,-1):
                cl = []
                com = it.combinations(pl,(n-i))
                for each in com:
                        tl = list(each)
                        for j in range(i):
                                tl.append(d)
                        cl.append(tl)
                plst = set()
                for ecom in cl:
                        al = it.permutations(ecom,n)
                        for a in al:
                                s = 0
                                for j in range(n):
                                        s = s*10+a[j]
                                if isPrime(s) and int(math.log10(s)) == n-1:
                                        plst.add(s)
                if len(plst) > 0:
                        break
        return plst

print (sum(sum(pr(4,k)) for k in range(10)))
print ('Time used: %.2f' % (time.time()-start))

jerryxjr1220

对算法进行了优化

# -*- coding: utf-8 -*-
from itertools import combinations

import time
start = time.time()
def getPrime(n):
        primes = [True]*n
        primes[0], primes[1] = False, False
        for (i,prime) in enumerate(primes):
                if prime:
                        for j in range(i*i,n,i):
                                primes[j] = False
        plist = []
        for (k,trueprime) in enumerate(primes):
                if trueprime:
                        plist.append(k)
        return plist
primes = getPrime(100000)
def isPrime(n):
        if n > 100000:
                for i in primes:
                        if n%i==0:
                                return False
                        if i*i > n:
                                break
                return True
        else:
                if n in primes:
                        return True
                else:
                        return False
 
def generateDigitLists(n,without, withZero=False):
    if n == 0:
        yield []
    else:
        start = 0 if withZero else 1
        for lst in generateDigitLists(n-1,without, True):
            for i in range(start,10):
                if i != without:
                    yield [i] + lst

def F(n,d):
    for k in range(n-1,0,-1):
        resultSet = set()
        start = 0 if d != 0 else 1
        for dInds in combinations(range(start,n), k):
            for dList in generateDigitLists(n-k,d, d != 0 and 0 in dInds):
                j=0
                lst = [d]*n
                for i in range(n):
                    if i not in dInds:
                        lst[i] = dList[j]
                        j += 1
                number = int(''.join(str(dig) for dig in lst))
                if isPrime(number) :
                    resultSet.add(number)
        if len(resultSet) > 0:
            return k,resultSet


def MNS(n,d):
    k,rs = F(n,d)
    return k,len(rs),sum(rs)

def e111(n):
    r = 0
    for d in range(10):
        m,ln,s = MNS(n,d)
        r += s
    return r

print (e111(10))

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

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <omp.h>
using namespace std;

const int nLen = 10;//讨论的素数位数

char cPrimes[100104];//cPrimes[i]=1表示i是一个素数
int nPrimes[9600];//具体记录素数
int nps = 0;//记录在0到100100之间素数的个数
int nCount[10];//nCount[i]记录N(nLen,i)的值
long long lTSum[10];//lTSum[i]记录S(nLen,i)的值
int nMaxCount[10];//nMaxCount[i]记录M(Len,i)的值

void dfs(int nCt, //当前未变动的数字的个数，也是取指定数字的个数
  int nPos, //填充的位置
  const char* str, //逐渐被填充的字符串（最后形成欲求的素数）
  int nChkInt)//指定的数字
{
  if (nCt == nMaxCount[nChkInt])//所有可变动的指标已经用完
  {
    //检查最后生成的数字是否是一个素数
    long long x = 0;
    for (int i = 0; i < nLen; ++i)
    {
      x *= 10;
      x += str[i];
    }
    int nD = sqrt((double)x);
    bool bIsPrime = true;
    for (int i = 0; nPrimes[i] <= nD; ++i)
    {
      long long p = (long long)nPrimes[i];
      if (x%p == 0)
      {
        bIsPrime = false;
        break;
      }
    }
    if (!bIsPrime)
      return;//不是素数
    ++nCount[nChkInt];//保存结果
    lTSum[nChkInt] += x;
    return;
  }
  if (nPos >= nLen)//超出可填充的位置
    return;
  for (int i = 0; i < 10; ++i)//继续枚举在第nPos位置上的数字
  {
    char str1[16];
    memcpy(str1, str, sizeof(str1));
    str1[nPos] = i;
    if (i == nChkInt)
      dfs(nCt, nPos + 1, str1, nChkInt);//如果是指定值则未占用变动指标
    else
      dfs(nCt - 1, nPos + 1, str1, nChkInt);//如果不是指定值则占用变动指标
  }
}

int main(void)
{
  double t = omp_get_wtime();
  //筛法求批量素数
  memset(cPrimes + 2, 1, sizeof(cPrimes) - 2);
  for (int i = 2; i <= (int)sqrt(100100.0); ++i)
  {
    if (cPrimes[i] == 0)
      continue;
    for (int j = i*i; j < 100100; j += i)
      cPrimes[j] = 0;
  }
  //保存具体的素数
  nPrimes[nps++] = 2;
  for (int i = 3; i < 100100; i += 2)
  {
    if (cPrimes[i] == 1)
      nPrimes[nps++] = i;
  }
  long long lSum = 0;
#pragma omp parallel for reduction(+:lSum)
  for (int nChkInt = 0; nChkInt < 10; ++nChkInt)//求M(Len,nChkInt)，N(Len,nChkInt)，S(Len,nChkInt)
  {
    nMaxCount[nChkInt] = nLen;
    while (nCount[nChkInt] == 0)
    {
      --nMaxCount[nChkInt];
      char str[16];//记录最后满足条件的素数
      memset(str, 0, sizeof(str));
      for (int i = 1; i <= 9; ++i)//枚举要生成的素数的第0位上的数字
      {
        int k = nLen;
        memset(str, nChkInt, nLen * sizeof(char));//把该数所有位置上的值都填成nChkInt
        str[0] = i;
        if (i != nChkInt)
          --k;//未取指定的值，可变动数字的个数减一
        dfs(k, 1, str, nChkInt);//逐步填充后面的位置上的数字
      }
    }
    lSum += lTSum[nChkInt];
  }
  //输出结果
  for (int i = 0; i < 10; ++i)
    cout << "M(" << nLen << "," << i << ") = " << nMaxCount[i] << "\tN(" << nLen << "," << i << ") = " << nCount[i] << "\tS(" << nLen << "," << i << ") = " << lTSum[i] << endl;
  cout << "Sum = " << lSum << endl;
  t = omp_get_wtime() - t;
  cout << "耗时：" << t << "秒" << endl;
  return 0;
}