#本題:Project Euler 003
003

13195 の素因数は 5, 7, 13, 29 である.
600851475143 の素因数のうち最大のものを求めよ.

->次の問題

#　できました

cheat.py

from sympy import factorint

n = 600851475143
prime_factor = factorint(n)
print(max(prime_factor))

はい、sympy使うと一瞬ですね。実行速度もこれが一番早かったです。
さすがにずるいので素因数分解も自分で作ってみます。無駄な処理を回すと計算にものすごく時間がかかってしまうのでできるだけ無駄を省くようにアルゴリズムを考えてみます。

考え方

例えばn = 60を素因数分解するときどうするか考えてみます。
60は2で割れる　→　素因数のリストに2を加え[2,]、60を2で割る(=30) 30は2で割れる　→　　　　　　〃　　　　　[2,2,]、30を2で割る(=15) 15は2で割れない。次に大きな素数は3。 15は３で割れる　→　素因数のリストに3を加え[2,2,3,]、15を3で割る(=5) 5は3で割れない。次に大きな素数は5。 5は5で割れる。　→　素因数のリストに5を加え[2,3,3,5,]、5を5で割る(n=1) nが1になったら終了。素因数=[2,3,3,5]
nを小さい素数から順番に割っていき、n=1になるまで繰り返すことで素因数分解できそうです。
そのためには素数のリストが必要になるので素数のリストを生成する関数を作成します。

###追記
コメントで教えていただきましたが、素数を作って素因数を求めるやり方は素数の生成に時間がかかってしまい、2以上の整数で順番に割っていくやり方のほうが早かったという結果になりました＼(^o^)／

def simple_prime_factor_generator(n: int) -> int:
    """
    シンプルな素因数分解のgenerator
    2から順番に割っていくだけ、でも断然早い
    :param n:
    :return:
    """
    for i in range(2, ceil(sqrt(n))):
        if i ** 2 > n:
            break
        while n % i == 0:
            yield i
            n //= i
    if n != 1:
        yield n

結果 n = 600851475143 [71, 839, 1471, 6857] prime_factor_generator 処理時間: 0.08188509941101074 sec. [71, 839, 1471, 6857] simple_prime_factor_generator 処理時間: 0.0006339550018310547 sec.
素数を生成してる間に整数でさっさと割って割れるかどうか見るほうが早いということでしょうか。
さらにシンプルな方はrange(2, int(sqrt(n))で偶数も含めたまま計算をしているのでここを奇数だけにしたり、さらに3の倍数を省いたりするだけでさらに早くできそう。
素数生成はせっかくなので他のところで使ってみよう。

追記②

素数生成関数を改良し、素数生成を高速化しました。ようやくsimpleな方法より早くなりました。

math_function.py

import numpy as np
from math import ceil, sqrt

def prime_factor_generator(n: int) -> int:
    """
    素因数分解するgenerator。
    :param n: int
    :return: int
    """
    prime_list = list(get_prime_arr(ceil(n ** 0.5) + 1))
    prime = 0
    while prime ** 2 < n:
        prime = prime_list.pop(0)
        while n % prime == 0:
            n //= prime
            yield prime
    if n != 1:
        yield n

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

素数を作ってみる

euler003.py

import numpy as np

def prime_generator(n: int) -> int:
    """
    n以下の素数を返すgenerator
    :param n: int
    :return: int
    """
    # 2からnまでの配列を作成
    num_arr = np.arange(2, n + 1)
    while num_arr.size > 0:
        # 配列の先頭の数字(最も小さい)を素数として返す。
        prime = num_arr[0]
        yield prime
        # 配列の中から上で出した素数で割れるもの(=素数でないもの)を削除
        # このとき上で出した素数自身も割れるので削除される
        num_arr = num_arr[num_arr % prime != 0]

有名なエラトステネスの篩です。2から任意の自然数までの配列を作り、
先頭の数字を取り出す→取り出した数字で割れる数を配列から除外
を繰り返すことで素数の配列を得ます。

generatorはforのループや内包表記で使ったり、next(generator)とすることで中身を取り出せます。

generator = prime_generator(10)
for p in generator:
    print(p)
#->2 3 5 7

generator = prime_generator(10)#もう一度使うときはgeneratorを作り直す
prime_list = [p for p in generator]
print(prime_list) 
#->[2, 3, 5, 7]

generator = prime_generator(10)
print(next(generator))
#->2 
print(next(generator))
#->3
for p in generator:#<-generatorを初期化せずに使い回すと
    print(p)
#->5 7・・・上で2,3を取り出しているので5から始まる

##素因数分解にはどこまで大きさの素数が必要？
例えば169を素因数分解するときに169までの素数を作るのは無駄と考えられます。13までの素数があれば十分です。
ある数字nが素数ではないとき、n = a * bで表せます(aは素数,bは2以上の自然数, a<=b)。
割れるかどうかは小さい素数から順に検討するのでaが最も大きくなる場合を考えれば良さそうです。
aが最も大きくなるのは a=bの場合、つまりnがなにかの2乗の場合です。
よって、素因数分解するには、≦√nまでの素数があれば良いと考えます。≦√nまでの素数で一度も割り切れなければnは素数ということになります。
√nの計算にはmathライブラリーのmath.sqrtが使用できます。また、自作のprime_generatorでは引数は自然数である必要があるので、√nを切り上げるためのmath.ceilを使用します。

√n ← math.ceil(math.sqrt(n))
##素因数分解してみる
はじめに書いた、「nを素数で繰り返し割っていく」をコードにします。

euler003.py

def prime_factor_generator(n: int) -> int:
    """
    自然数nの素因数を返すgenerator
    :param n: int
    :return: int
    """
    generator = prime_generator(ceil(sqrt(n)))#√nまでの素数を生成するgenerator
    for prime in generator:#素数を生成しprimeに代入
        if prime ** 2 > n:# nは割られてどんどん小さくなるのでprimeの2乗がnを上回ったらループを抜ける
            break
        while n % prime == 0:#primeで割り切れる限りはprimeを素数として出力し、nを割る。
            n //= prime
            yield prime
    if n != 1:#ループを抜けたときにn=1の場合があるのでそのときは値を返さない。
        yield n#最後に残ったnは素数なので返す。

コード

素因数を返す関数ができたので、あとはmax()で最大値を返すのみ。

euler003.py

from math import ceil, sqrt
import numpy as np
from stop_watch import stop_watch


@stop_watch
def main():
    n = 600851475143
    prime_factor = set(p for p in prime_factor_generator(n))
    print(max(prime_factor))


def prime_factor_generator(n: int) -> int:
    """
    自然数nの素因数を返すgenerator
    :param n: int
    :return: int
    """
    generator = prime_generator(ceil(sqrt(n)))
    for prime in generator:
        if prime ** 2 > n:
            break
        while n % prime == 0:
            n /= prime
            yield prime
    if n != 1:
        yield int(n)


def prime_generator(n: int) -> int:
    """
    n以下の素数を返すgenerator
    :param n: int
    :return: int
    """
    num_arr = np.arange(2, n + 1)
    while num_arr.size > 0:
        prime = num_arr[0]
        yield prime
        num_arr = num_arr[num_arr % prime != 0]


if __name__ == '__main__':
    main()

結果 6857 main処理時間: 0.28844213485717773 sec.
素因数を返す関数にfor, while, ifがごちゃっとあって美しくない。もっと綺麗にできるはず。