【Project Euler】Problem 12: 約数の多い三角数

• 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題 12. 約数の多い三角数

原文 Problem 12: Highly divisible triangular number

問題の要約：約数の数が初めて500を超える三角数を求めよ

まず三角数(Wikipedia)にあるように$n$番目の三角数$T_n$は$1$から$n$までの和なので、

T_n = n(n+1)/2

まずsympyのdivisorsを使うと、とても簡単プログラムになります。

import sympy
import itertools

for n in itertools.count(1):
  dn = len(sympy.divisors(n*(n+1)//2))
  if dn>500:
    break
print(f" n={n}, Tn= {n*(n+1)//2}, divisors={dn}")

ただし、これだとすべての三角数の約数を求めるので速いとは言えません。

ここで約数の個数(Wikipedia)を求める公式を使います

Nの素因数分解を \\
\large N = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k} \\
とし、約数の個数を d(N)とすると\\
d(N)=(a_1+1)(a_2+1) \dots (a_k+1)

要は素因数分解さえできれば約数を全部列挙しなくても数だけを求めることができるということです。素因数分解はsympyのfactorintでDICT形式で求めることが出来ます。

from sympy import factorint
for i in [4,5,6]:
  print(sympy.factorint(i))
#{2: 2}
#{5: 1}
#{2: 1, 3: 1}

約数の個数の公式を元にfactorintで素因数分解されたDICT形式から約数の個数を求める関数fct2divnは以下のようになります。

import numpy as np
# Factors to number of disivors
def fct2divn(fct):
  return np.prod(np.array(list(map(lambda x: x+1,list(fct.values())))))

print(fct2divn({2: 1, 3: 1}))
# 4

そこで問題の三角数$T_n$の約数の数ですがポイントは以下になります。

• $n$と$n+1$は共通の素因数がない（互いに素）
• 従って$d(n \times (n+1)) = d(n) \times d(n+1)$
• $n$と$n+1$はどちらかが偶数なのであらかじめ２で割っておける
• 例えば$n$が偶数の時$d(Tn) = d(n/2) \times d(n+1)$となる

これをプログラムに実装すると以下のようになります。divisorsとfactorintのプログラムの実行時間の比較表を見るとかなり速くなっているのが分かります。

dn = 1
for n in itertools.count(1):
  n1 = (n+1)//2 if (n+1)%2 == 0 else n+1 
  dn1 = fct2divn(factorint(n1))
  dnt = dn * dn1
  if dnt>500: 
    break
  dn = dn1
print(f" n={n}, Tn= {n*(n+1)//2}, divisors={dnt}")

N	500	1000	2000	3000	5000
divisors	0s	5s	133s	x	x
factorint	0s	0s	8s	13s	90s