Project Euler 010を解いてみる。「素数の和」

Project Euler 010

010

10以下の素数の和は 2 + 3 + 5 + 7 = 17 である.
200万以下の全ての素数の和を求めよ.

->[次の問題]

Sympyにありました。

n = 200 * 10000
# sympy.primerange(a, b)でa以上b未満の素数を生成するgeneratorを返す
prime_generator = sympy.primerange(2,n+1)
print(sum(prime_generator))

ただ、この書き方だと下記のエラトステネスの篩のほうが早いです。

考え方

シンプルな問題ですね。
エラトステネスの篩で200万以下の素数を計算して合計すれば良さそうです。
篩のプログラムは以前に作ってmath_functions.pyに入れてあるのを使います。
なにかに応用できるかなと思ってboolean配列で返す関数と素数列で返す関数に分けてあります。
配列の数字を付け加えたり消したりといった処理やifによる分岐は重くなるようなので、boolean配列のTrue, Falseを切り替えることで値を表し、最後に数列に変換しています。
また、2,3以外の素数は6k±1であるので、少し探索範囲を減らしています。

boolean_arr[4::2] = False
#これはboolean_arrインデックス4から2個ずつ(4,6,8,...)をFalseにするという処理。
#bookean_arr[i*2::i] = Falseとすることでiより大きいiの倍数のインデックスをFalseにできる=非素数を消せる
#if文を使わないので多分早い

コード

math_functions.py

import numpy as np


def get_prime_arr_bool(n: int) -> np.ndarray:
    """
    自然数n未満の素数をboolean配列で返す関数
    2,3,5 -> [False, False, True, True, False, True]
    :param n: int
    :return: numpy.ndarray
    """
    boolean_arr = np.ones((n,), dtype=np.bool)  # n個の要素をもつのboolean配列(0~n-1に対応)を生成、全てTrueとする
    boolean_arr[0:2] = False  # [0,1]をFalseとする
    boolean_arr[4::2] = False  # 2以外の偶数をFalseとする
    boolean_arr[6::3] = False  # 3以外の3の倍数をFalseとする
    for i in range(5, ceil(sqrt(n)), 6):  # iは初項5、公差6の等差数列とする
        if boolean_arr[i]:  # 6k-1に対応。iに一致する要素がTrue=素数のとき
            boolean_arr[i * 2::i] = False  # iの倍数の要素(非素数)をFalseとする
        if boolean_arr[i + 2]:  # 6k+1に対応。
            boolean_arr[(i + 2) * 2::(i + 2)] = False
    return boolean_arr  # boolean配列を返す


def get_prime_arr(n: int) -> np.ndarray:
    """
    自然数n未満の素数の配列を返す関数
    :param n: int
    :return: numpy.ndarray
    """
    return np.arange(n)[get_prime_arr_bool(n)]  # 0~n-1の配列を生成、boolean配列でTrueとなっている数のみを返す

euler010.py

import numpy as np
from math_functions import get_prime_arr


def main():
    n = 200 * 10000
    prime_arr = get_prime_arr(n+1)  # 素数配列を生成(n未満を生成なので一応+1しておく)
    print(np.sum(prime_arr)) # 合計


if __name__ == '__main__':
    main()

結果 142913828922 main処理時間: 0.020266056060791016 sec.
できればforループの中のifもなくしたいところですが、どうやったらifなしで実現できるのか思いつかないので現状はこれが精一杯です。