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

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

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

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

3つの自然数の和が1,000のピタゴラス数は?:Project Eulerの問題9をC言語で解き、その実行時間を検証しました!

序章

今日のProject Eulerの問題は「自然数a,b,c (a<b<c)が、a+b+c=1000 かつ a^2+b^2 = c^2を満たすとき、abcの値を求めよ」と言うものです。簡単に言えば3つの自然数の和が1,000となるピタゴラス数を探す事が主題となります。

ピタゴラス数とは

ピタゴラスとは「3つの自然数a,b,cが、a^2+b^2 = c^2を満たすもの」を言います。ここで直角三角形の3辺の長さを思い浮かべてもらえば分かる通り、(a,b,c) =(3,4,5)(a,b,c)=(5,12,13)などは上の式を満たしますね。以下英語版の問題を見てみましょう。

今回の問題

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,a^2 + b^2 = c^2 For example, 3^2 + 4^2 = 9 + 16 = 25 = 5^2.There exists exactly one Pythagorean triplet for which a + b + c = 1000.

Find the product abc.(Special Pythagorean triplet:Project Euler.netより)

第一章:まずは数学的に問題を捉える

まず本問題はピタゴラス数のパラメーター表示を用いる事で手計算(人力)で解く事も可能です。人力で解ける事を読者に確認してもらう為に、以下の解答を掲載しようとおもいます。プログラムのみに興味のある方は読み飛ばしても差し支えありません。

不定方程式「a^2+b^2=c^2」-(a) と 「a+b+c=1000」 -(b)をともに満足する自然数の解(a,b,c)*但し0<a<b<cを考える。尚、不定方程式とは「方程式の個数<未知数の個数を表す方程式」の事で、解が定まらないものを言う。

式(a)はピタゴラス数を表す方程式である。このためa,b,cm,nを用いて以下の(I)(II)のようにパラメータ表示できる。*1

(I) (a,b,c) = (2mn,m^2-n^2,m^2 + n^2)のとき、a = 2mnb = m^2-n^2c=m^2 + n^2 とパラメータ表示する。さらに(b) =>2m^2+2mn=1000 m^2+mn=500 <=> m(m+n) = 2^2*5^3。これを満たす(m,m+n) = (1,500),(2,250),(4,125),(10,50),(20,25)より、(m,n) = (499,1),(248,2),(121,4),(40,10),(20,5)

ここで(m,n) = (499,1)のときのcの値を計算すると499^2+1^2>1000となる。
これは(m,n) = (248,2),(121,4),(40,10)の場合も同様で(b)に反する。このことから(m,n)=(20,5)の場合に限られ、(a,b,c)=(200,375,425)を得る。

(II) (a,b,c) = (m^2-n^2,2mn,m^2 + n^2)のとき、(I)と同様にすると(a,b,c) = (375,200,425)を得るが、a<bに反する。

以上より(a,b,c)=(200,375,425)が解となりabcを求めると、abc= 31875000(答)*2

第二章:全探索法のプログラムを書き、総当たりで第一章の答えを検証

第二章では第一章の検算をプログラムに行なわせ、第一章の答えが合っているかを調べます。調べ方は、条件を満たすa,b,cを総当たりで探していきます。これはパスワードを総当たりで探す要領(パスワードクラック)と同じです。

今回はa+b+c=1,000と3つの数の和が1,000である事から単純に、1≦a≦1,000,1≦b≦1,000,1≦c≦1,000の範囲で調べます。このとき考えうる(a,b,c)の組の総数は(10^3)^3 = 10^9通りとなり、プログラムの計算回数も10^9回となります。以下のソースコードより三重ループなので、全探索における計算量のオーダーは O(n^3)となります。

あとC言語採用の理由を言うと、JavaよりC言語慣れしてるからです。もっとこじつけるとコンパイル言語なので実行時間が速い点があります。問題に行列演算などが必要な場合は実装の面でnumpy(Python)を使いたい所ですが、今回のように不定方程式の解の候補を探す程度ならC言語かなと思いました。

調査用のプログラム

プログラムの実行結果と第一章の検算

プログラムの実行結果は以下のようになります。最初の行の(a,b,c)が第一章の解答の(I)の場合、次の行の(a,b,c)が第一章の解答の(II)の場合に当たる事が分かります。答えの数値も31,875,000であっている事が分かります。


f:id:program_study:20150221194314p:plain

第三章:プログラムの実行時間について調べる

第二章では第一章の解答と照らし合わせ、実際に答えがあっていることがわかりました。第三章ではこのプログラムの実行時間について調べていきます。今回は筆者の練習も兼ね、少々統計学的なやり方で検証します。まずは検証に必要な統計的な基本についておさらいしていきます。

母集団と標本

例えば今回の場合プログラムを実行する実験(試行)において、実行時間全体の集合を母集団と言います。この内n回実行しその実行時間をt_1,_t_2,t_3,….,t_nとしたとき、t_1〜t_nのような母集団から一部分を抜き出した部分集合を標本と言います。また今回のように母集団の一部(部分集合)を取り出して調査することを標本調査と言います。

標本平均、標本分散、標本標準偏差

今回の場合もこの標本調査に当たるが、各試行における標本の値をX_1,X_2,X_3,…….,X_nとしたときその標本平均\bar{x}標本分散s^2は以下のように定義されます。又標本分散の平方根を取った値を標本標準偏差sと言います。


f:id:program_study:20150221193758p:plain

nを限りなく大きくすると、標本平均\bar{x} \rightarrow μと母平均に収束し、標本平均の分散はV(\bar{x}) \rightarrow 0に収束します。この事から試行回数nを大きくすればする程真の平均に近づきます*3

この事から計測した実行時間の標本平均から最終的な実行時間、標本標準偏差から各試行における実行時間の誤差を見て行きます。

プログラム実行時間の流れについて

第二章のプログラムの実行時間計測について話します。計測に使うPCは「macbook air(10.9.5) 1.6GHz Intel Core i5/2GB Memory」となります。実験の回数もとい試行回数(サンプルサイズ)は5回(n=5)とし、その標本平均、標本標準偏差をみていきます。各実行時間はC言語のclock関数を使い、ループにかかる時間を計測していきます。

プログラム実行時間の集計結果

以下Excelで集計した結果を示します。まずプログラムの実行時間の標本標準偏差は0.00759788より、小数点第4位の所で誤差が発生していると考えられます。この為標本平均を小数点第4位の所で四捨五入して考えると標本平均は3.086[秒]となり、プログラムの実行時間の目安は3.086[秒]となります。


f:id:program_study:20150221192649p:plain

第四章:第二章のプログラムを改良し、再度実行時間について調べる

今回のケースでは総当たりで不定方程式の解を探しました。しかし実際にはnが大きく計算が追いつかない場合も多く、プログラムの手順(アルゴリズム)を改善していく必要があります。この方法で本一冊位になったり、「TopCoder.com」のようなアルゴリズムを競うコンテストも世の中にはあったりします。

一応インターネットで計算量を下げる方法を探してみましたが見つかりませんでした。この為最も簡単なループを見直す事で、探索回数を減らす事を考えました。問題文のa<b<cを加味し、プログラムのループを以下のように改良しました。

改良後の実行時間について

そして改良後の実行時間は以下のようになります。


f:id:program_study:20150221192712p:plain

再度実行時間の標本平均を見てやると、約0.557[秒]である事が分かります*4。又改善前の実行時間の標本平均÷改善後の実行時間の標本平均 = 3.086[秒] ÷ 0.557[秒] = 約5.5[倍]の速度向上に成功しました。自画自賛ですが条件1つ変えるだけで、ここまで速度向上するとは思いませんでした!

最後は雑談で!

いきなりですがみなさんここまで読んでて疲れましたね。自分も楽しいんですけど書いてて疲れました。自分の疲れを癒す目的で、こんな物を作ってみました。いつもはアニメ「ラブライブ!」のキャラクターですが、今回はプログラムにちなんで「プロ生ちゃん(暮井慧)/CheeseDeHappy氏より」を起用してみました。*5


f:id:program_study:20150221192810j:plain

今回はプログラムの実行速度を考える上で理論より実践よりの記事になりました。実践や実験の好きな人向けの記事になったと自画自賛しています。自分の統計の復習のために実行時間の標本平均、標本標準偏差を取ってみました。合っていれば幸いですが…。

数学の話から離れますが最後のグラフを作ってて、自作PCベンチマーク結果のグラフを思い出します。自分の記憶だとPC自作雑誌のベンチマークなんかも3回計測して平均値を出してましたね。今こそプログラムの関係でmacですがDOS/V POWER REPORTやWinPCなんかを読んで、秋葉原に行ったあの頃が懐かしいです。それではみなさん次の記事まで。

参考にしたWebサイト

今回の記事作成にあたり、多くの大学のサイトの文献を参考にしました。大学卒業するとこういうサイトの存在が懐かしくなります。

参考文献/参考書籍

まずは日本に住む皆さんおなじみのこの書籍。大学で一度統計学を学び終えた人におすすめの一冊です。標本標準偏差などの数学的な意味を簡単に説明しています。ざらっと数学的意味合いを確認したい人向けです。

統計学入門 (基礎統計学)

統計学入門 (基礎統計学)

今回の記事作成にあたり用意したけど結局使わなかったこの書籍。数学的な説明を抑えめに、実際にどのような事なのかを挿絵等を使って分かりやすく説明しています。理論より実践に重きをおきたい人向けです。

Statistics

Statistics

*1:大学受験の問題で余り使った記憶はありませんが、数学の本だと度々登場するピタゴラス数のパラメータ表示。自然数a,b,cに対しa^2+b^2 = c^2を満たすとき、a = m^2 - n^2b = 2mnc = m^2+n^2とパラメータ表示できます。色々解説サイトもありますが、その中でおすすめとして「ある種のピタゴラス数と2次形式・課題研究から/富山県立砺波高等学校 片山 喜美:http://www1.tst.ne.jp/ja9nfo/math/kiyou.pdf」があります。

*2:この計算を手でやると結構大変です。まず200,375,425を素因数分解して10の倍数や2乗を外に出して計算してみて下さい。尚cを消去する事でパラメータ表示せずに解く事もできますが、計算が大変でした。

*3:詳細は「基礎統計学I 統計学入門(東京大学出版会)」のP183〜P184を参照して下さい

*4:小数点は改良前に合わせました

*5:利用規約を見てみましたが、本記事は非商用の記事ですので多分プロ生ちゃんの利用規約に反しないとは思うんですが...