Project Euler 200 / 連続する部分文字列に200を含む耐素数性スキューブ

問題概略

（, は異なる素数）で表せる数をスキューブ（sqube）と定義する。
, はスキューブである。
最初の 5 つのスキューブは 72, 108, 200, 392, 500 である。

200 はどの位の数字を変更しても素数とならない最小の数である。
この特徴をもつ数字を「耐素数性のある（prime-proof）数」と呼ぶ。
連続する部分文字列に「200」を含む次の耐素数性のあるスキューブは 1992008 である。

連続する部分文字列に「200」を含む 200 番目の耐素数性のあるスキューブを求めよ。

[https://projecteuler.net/problem=200

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

mathematicaで解く

番目の素数を であらわして とおきます。

, の上限を適当に定めて全探索すれば答えは出せそうですが，この方法だと と がアンバランスな場合を考えていないことになるので厳密には正解と言えないと思います。

そこで，「 の上限 を決める」→「 の上限を決める」→「 を動かしてチェック」という手順で解きました。

で のときを考えると の上限がわかります。 より

を固定したときの の上限も同様です。 より

連続する部分文字列に「200」を含むかどうかのチェックは簡単です。
文字列に直して，その部分文字列に「200」を含むかどうか調べるだけです。

    contain200Q[n_] := contain200Q[n] = ! StringFreeQ[ToString@n, "200"];

耐素数性のチェックは引数として受け取った数の 桁目を に変える関数を作って，すべての について調べました。
首位の数字を 0 に変えることはないので，下のコードでは Rest でその を除いています。

    primeproofQ[n_] := primeproofQ[n] = Module[{lst = IntegerDigits@n, t},
         t[{k_, x_}] := FromDigits@ReplacePart[lst, k -> x];

あとは全探索するだけ。問題文中に与えられている 1992008 が 7 桁なので， からはじめて条件をみたす数が 200 個以上になるまで「 を 10 倍して再計算」を繰り返しました。*1
計算時間は約 1.5 秒です。

    In[]:= AbsoluteTiming[
        nmax = 10^7;
        jmax = PrimePi[(nmax/4)^(1/3)];
        f[{i_, j_}] := f[{i, j}] = Prime[i]^2*Prime[j]^3;
        g[j_] := Module[{imax, lst},
          imax = PrimePi[Sqrt[nmax/Prime[j]^3]];
          lst = f[{#, j}] & /@ Select[Range@imax, # != j &];
          Select[lst, And[contain200Q@#, primeproofQ@#] &]];
        contain200Q[n_] := contain200Q[n] = ! StringFreeQ[ToString@n, "200"];
        primeproofQ[n_] := primeproofQ[n] = Module[{lst = IntegerDigits@n, t},
           t[{k_, x_}] := FromDigits@ReplacePart[lst, k -> x];
           NoneTrue[t /@ Rest@Tuples[{Range@Length@lst, Range[0, 9]}], 
            PrimeQ]];
        While[Length@(Flatten[g /@ Range@jmax]) < 200,
         nmax *= 10;
         jmax = PrimePi[(nmax/4)^(1/3)]];
        ans = (Sort@(Flatten[g /@ Range@jmax]))[[200]];]
       
    Out[]= {1.55457, Null}

pythonで解く

同じことを python でもやりました。計算時間は約 0.8 秒と mathematica の約半分で，思ったほど速くはなりませんでした。この解法は python らしくないんでしょうね。

    import time, math
    from sympy import isprime, prime, primepi, sieve
    from functools import lru_cache
    import itertools
    
    @lru_cache(maxsize=None)
    def f(pi, pj):
        return pi ** 2 * pj ** 3
    
    @lru_cache(maxsize=None)
    def contains200(n):
        return "200" in str(n)
    
    @lru_cache(maxsize=None)
    def is_prime_proof(n):
        nn = str(n)
        lst = list(itertools.product(
            list(range(1, len(nn) + 1)), list(map(str, range(0, 10)))))
        return all([not isprime(int(nn[:k - 1] + x + nn[k:])) for k, x in lst[1:]])
    
    def solve():
        nmax = 10 ** 7
        pjmax = prime(primepi(math.pow(nmax // 4, 1 / 3)))
    
        def g(pj):
            pimax = prime(primepi(math.sqrt(nmax // (pj ** 3))))
            return [f(pi, pj) for pi in sieve.primerange(pimax + 1)
                    if pi != pj and contains200(f(pi, pj)) and is_prime_proof(f(pi, pj))]
    
        while len(sum(list(map(g, sieve.primerange(pjmax + 1))), [])) < 200:
            nmax *= 10
            pjmax = prime(primepi(math.pow(nmax // 4, 1 / 3)))
    
        return sorted(sum(list(map(g, sieve.primerange(pjmax + 1))), []))[199]
    
    if __name__ == '__main__':
        start = time.time()
        print(solve())
        elapsed_time = time.time() - start
        print("elapsed_time: {0} [sec]".format(elapsed_time))

*1:101問以降なので答えの数値は省略。