黄昏より暗きもの、血の流れより赤きもの

読者です 読者をやめる 読者になる 読者になる

黄昏より暗きもの、血の流れより赤きもの

自分の好きな事を好きなように書いて行きます。

約数の個数が初めて500を超える三角形数は?:Project Euler Problem 12 を解き、約数の個数をグラフ化してみた

序章

今日はProject Eulerの問題12を解決できたので、そのやり方について紹介したい。まず問題12は以下の通りである。

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...Let us list the factors of the first seven triangle numbers:
1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,1/
21: 1,3,7,21
28: 1,2,4,7,14,28
We can see that 28 is the first triangle number to have over five divisors.
What is the value of the first triangle number to have over five hundred divisors?(Problem 12 - Project Euler 12より)

第一章:問題を解く上での予備知識

約数と倍数*1

自然数a,bにおいて「aがbで割り切れる」ときaはbの倍数、bはaの約数と言う。例えば(a,b)=(6,2)のときa ÷ b = 6 ÷ 2 = 3で割り切れるから、6は2の倍数、2は6の約数となる。

2つの自然数a,bに対する最大公約数GCD(a,b)と最小公倍数LCM(a,b)

自然数a,bに対して、cがa,bの公倍数であるとは、cがaの倍数であり、同時にcがbの倍数である事を言う。公約数の場合も同様に、cがa,bの公倍数であるとは、cがaの約数であり、同時にcがbの約数である事を言う。

ここで公倍数のうち最小の物を最小公倍数と言い、以下LCM(a,b)と表す。同様に公約数のうち最大の物を最大公約数と言い、以下GCD(a,b)と表す。ここでGCD(a,b)=1となるa,bをaとbは互いに素と言う。

素因数分解と約数の個数

まず素因数分解とは、自然数nを素数の積で表す事であった。ここで、

2以上の自然数N素数p_1,p_2,…..,p_nとその指数r_1,r_2,….,r_nを用いて、N=p_1^{r_1}*p_2^{r_2}*……*p_n^{r_n}と一意に素因数分解できる。(整数論の基本定理)*2

このとき約数の個数S(N)p_nを含まない場合を考え、S(N)=(r_1+1)(r_2+1)….(r_n+1)(指数+1)と表せる。例えば24の約数の個数S(24)は、24=2^3*3^1素因数分解できるためS(24)=(3+1)(1+1)=8個となる。

第二章:プログラムを実行してみて、その結果を検証してみた

1:まずは自分のプログラムの実行時間を計測

実際にプログラムを実行してみた。まず三角形数とは1から順々に自然数を足し合わせた数の事を言い、自然数nにおける三角形数の総和はs(n)=\frac{n(n+1)}{2}となる。本プログラムはn=12,375のとき、s(12,375)=76,576,500 を記録している。この時の約数の個数S(76,576,500)=576(個)となる。

一方実行時間は0.6346[s]であった*3。尚解答作成に使ったプログラムは以下である。


2:次にnにつれて、三角形数の約数がどう変化するか見てみた

今度はJavaプログラムを改造し、n=1,2,3....のときの三角形数の約数の個数をcsvファイルに書き出してみた。*4その結果をR言語に読み込み、色々と検証してみた。

三角形数の約数の個数のグラフ

まず約数の個数をplot関数でグラフ化すると以下のようになる。

plot(data[,2],type="l",col=rgb(0,0.7,0.3))


f:id:program_study:20150304231113p:plain

三角形数の約数の個数のヒストグラム

次にヒストグラムで約数の個数の登場する頻度について調べてみた。すると指数曲線的に頻度は少なくなって行く事が分かった。約数の個数が100個超える事もまれなことも分かった。

plot(data[,2],type="l",col=rgb(0,0.7,0.3))


f:id:program_study:20150304231145p:plain

ちなみにその四分位数、平均値、母標準偏差をそれぞれ計算すると以下の通りとなる。n<=13,000の範囲での平均値は42.41311、母標準偏差は45.90164となる。四分位数を見ると16個間隔で約数の個数が分布している事もわかる。

> data <- read.csv("test.csv") #csvファイルの読み込み
#csvファイルの2列目を扱いたいときはdata[,2]のように,を付けてベクトルに変換
> summary(data[,2]) #四分位数
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 16.00 32.00 42.41 56.00 576.00
> mean(data[,2]) #平均値
[1] 42.41311
> varp = function(x){var(x)*(length(x)-1)/length(x)} #母分散を計算
> sqrt(varp(data[,2])) #母標準偏差を計算
[1] 45.90164

最終章

今回は既に解答例が出回っていたので、三角形数の約数の個数の分布まで調べてみた。Project Eulerの問題を解くだけでなく、グラフで図示してみるなどすると他の問題に応用できる糸口が見つかるかも知れない。

おまけ:他の皆さんはどう解いているの?

今回の問題の 解答作成にあたり、各種関連サイトを参考にさせてもらった。今回の問題解決の鍵は約数の個数をどうやって計算するかですが、下の[★1]のブログさんのようにRuby素因数分解の関数を使ってそこから約数の個数を求めている。 又[★2]のブログのように、整数の性質を上手く使い、三角形数の約数に関する漸化式を作って高速に問題を解いているようだ。

Ruby

Python

Haskell

*1:研文書院「大学への数学IA」より説明を借りてきた

*2:整数論の基本定理の証明はこちら

*3:計測方法としてUNIXのtimeコマンドのuserの時間に着目している。プログラムを処理する時間を5回にわたって計測し、その標本平均値を最終的な実行時間としている。timeコマンドの詳細は「 timeコマンドでプログラムの実行時間を知る(Qiita)」

*4:csvファイル生成に当たって「CSVファイルを書き出す:Javaちょこっとリファレンス」を参考にした。