Project Euler 012

012

三角数の数列は自然数の和で表わされ, 7番目の三角数は 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28 である. 三角数の最初の10項は:
1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...
となる.
最初の7項について, その約数を列挙すると, 以下のとおり.
1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28
これから, 7番目の三角数である28は, 5個より多く約数をもつ最初の三角数であることが分かる.
では, 500個より多く約数をもつ最初の三角数はいくつか.

->次の問題

この問題は割と力押しでやっており、処理時間が1秒を切れていません…

考え方

自然数nの約数の数は、nを構成する素因数の乗数から求めることができます。
例えばn=60のとき、

60 = 2^2 \times 3^1 \times 5^1

となります。この乗数にそれぞれ+1したものの積が約数の数になります。

(2+1) \times (1+1) \times (1+1) = 12

実際に60の約数は[1,2,3,4,5,6,10,12,15,20,30,60]の12個で合ってますね。
関数の呼び出し回数が増えると速度が遅くなってしまいますが、可読性と再利用のため約数の数を返す関数を作ります。以前作った素因数を返すprime_factor_generatorを使っています。

euler012.py

from math_functions import prime_factor_generator

def divisors_count(n: int) -> int:
    """
    自然数nの約数の数を返す関数
    素因数の乗数+1の総積が約数の数になるという方法を使用
    :param n: int 自然数 
    :return: int 約数の数
    """
    # nの素因数のlistを作成 n=60なら[2, 2, 3, 5]
    prime_list = list(prime_factor_generator(n))  
    # countを使ってlistの中に何個を同じ素数があるか計算しそれに+1する
    counts = [prime_list.count(p) + 1 for p in set(prime_list)]  
    result = 1
    for c in counts:
        result *= c  # countsの総積を計算
    return result

math_functions

import numpy as np
from collections import deque


def prime_factor_generator(n: int) -> int:
    """
    素因数分解するgenerator。素因数を順番に返す
    :param n: int
    :return: int
    """
    prime_que = deque(get_prime_arr(int(n ** 0.5) + 2))
    prime = 0
    while prime ** 2 < n:
        prime = prime_que.popleft()
        while n % prime == 0:
            n //= prime
            yield prime
    if n != 1:
        yield n


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となっている数のみを返す

これを使って約数の個数を出していき、500を超えるまで検証します。力技です。
三角数は階差数列なので一般項から出すこともできますが、どちらにしろ小さい数から順に検証していくことになるのでシンプルに足していきます。
ちなみにsympyには約数の数を返す関数sympy.divisor_countがあります。

コード

euler012.py

from math_functions import prime_factor_generator


def main():
    n = 500
    triangle_number = 1  # 初項 1
    difference = 2  # 階差数列の差
    while divisors_count(triangle_number) <= n:
        triangle_number += difference
        difference += 1  # 差を1つずつ増やしていく
    print(triangle_number)


def divisors_count(n: int) -> int:
    """
    自然数nの約数の数を返す関数
    素因数の乗数+1の総積が約数の数になるという方法を使用
    :param n: int 自然数 
    :return: int 約数の数
    """
    # nの素因数のlistを作成 n=60なら[2, 2, 3, 5]
    prime_list = list(prime_factor_generator(n))  
    # countを使ってlistの中に何個を同じ素数があるか計算しそれに+1する
    counts = [prime_list.count(p) + 1 for p in set(prime_list)]  
    result = 1
    for c in counts:
        result *= c  # countsの総積を計算
    return result


if __name__ == '__main__':
    main()

結果 76576500 main処理時間: 1.590303897857666 sec.
このコードでは2秒前後、sympy.divisor_countを使っても0.5秒くらいかかってしまいます。
現在下から順番に総当りしている部分にどうにか高速化できる余地があるはずですが、思いつかないです…