円周率の中から連続する10桁の最初の素数を見つける

与えられた数列の中から連続するN桁の最初の素数を見つけるプログラム

def is_prime(n):
# nが素数かどうか判定するプログラム
    for p in range(2, int(math.sqrt(n)) + 1):
        if n % p == 0:
            return False # 素数でないならFalseを返す
    return True # 素数ならTrueを返す

def get_prime(str_num, N):
# 文字列として受け取ったstr_numからN桁の最初の素数を見つけるプログラム
    for i in range(0, len(str_num) - N):
        if is_prime(int(str_num[i : i + N])) is True:
            print(i,'-th,',  str_num[i:i + N] )
            break

実行結果

math.piを文字列にして上のプログラムを実行。

>> import math
>> str_pi2 = '3' + str(math.pi)[2:]
>> get_prime(str_pi, 10)
4 -th, 5926535897

追記

@shiracamusさんにもっとカッコイイ書き方を教えていただので追記。
元コードではforとifで判定していたが、allですっきりとかける。またsqrtもn**0.5にすることでmathをインポートせずにかける。

def is_prime(n):
    """nが素数かどうか判定する"""
    return all(n % p != 0 for p in range(2, int(n**0.5) + 1))

こちらも元コードではforとifで判定していたが、イテレータでワンライナーに。あと range(len(str_num) - N)のようにNを引くのを忘れていることに気づきました。(修正済み)

def get_prime(str_num, N):
    """文字列として受け取ったstr_numからN桁の最初の素数を見つける"""
    i = next(i for i in range(len(str_num) - N) if is_prime(int(str_num[i:i+N])))
    print(f'{i}-th, {str_num[i:i+N]}')

もっと長い桁数の素数を見つけたい場合

円周率の計算

math.piやnumpy.piの円周率では有効桁数が足りないので自分で計算する必要がある。そこで、ガウス＝ルジャンドルのアルゴリズムを用いてN桁の円周率を計算。また、decimalモジュールを使い正確に浮動小数点計算を行う。(じゃないとOverflowしちゃうよ！)

import math
from decimal import *
def get_pi(N):
# N桁の円周率をガウス＝ルジャンドルのアルゴリズムで計算するプログラム
    max_iter = 2*math.ceil(math.log2(N))
    getcontext().prec = N

    # 初期値の設定
    a_n = 1
    b_n = Decimal(1) / Decimal(2).sqrt()
    t_n = Decimal(1) / Decimal(4)
    p_n = Decimal(1) 

    for n in range(0, max_iter):
        a_next = Decimal((a_n + b_n) / 2)
        b_next = Decimal(a_n * b_n).sqrt()
        t_next = Decimal(t_n - p_n * (a_n - a_next) ** 2)
        p_next = Decimal(2 * p_n)

        a_n = a_next
        b_n = b_next
        t_n = t_next
        p_n = p_next

    return Decimal((a_n + b_n) ** 2) / Decimal(4 * t_n)

実行結果

10000桁の円周率を計算しその中から15桁の最初の素数を見つける。

>> pi = get_pi(10000)
>> str_pi = '3' + str(pi)[2:]
>> get_prime(str_pi, 15)
35 -th, 841971693993751