2013年2月16日土曜日

プロジェクトオイラー 問72

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問72「真既約分数の集合は何個の要素を持つか?」

nとdを正の整数として, 分数 n/d を考えよう. 
n<d かつ HCF(n,d)=1 のとき, 真既約分数と呼ぶ.

d ≦ 8について既約分数を大きさ順に並べると, 以下を得る:
1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8
この集合は21個の要素をもつことが分かる.

d ≦ 1,000,000について, 真既約分数の集合は何個の要素を持つか.

-----


注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。





私の解答例は以下です。畳んでいます。

def f(n):
	L = range(n+1)
	for i in xrange(2, n):
		if L[i]<i: continue
		for j in xrange(i, n+1, i):
			L[j] /= i
			L[j] *= i-1
	return sum(L[2:])

n = 1000000
print f(n)


分母nの真既約分数の個数は「nと互いに素なn未満の数」の個数と同じです。
これは問69のトーティエント関数φ(n)の定義そのものです。
つまり求める真既約分数の個数は、φ(2)からφ(1000000)までの和です。

φ(n)の計算ですが、
「nと互いに素なn未満の数」は、言い換えれば、1からnまでのn個の数字の中から素因数で割り切れる数を外していき最後まで残った数字の個数です。

例えば
nの素因数に2があれば、1からnのうち、2,4,6,...のように偶数を外すので
個数はn*(1/2)個減り、残りはn*(1 - 1/2)=n*(1/2)になります。

nの素因数に3があれば、1からnのうち、3,6,9,...のよう3の倍数を外すので
個数はn*(1/3)個減り、残りはn*(1 - 1/3)=n*(2/3)になります。

このとき外す数からカウントすると、例えば6なら2の倍数と3の倍数で重複してカウントされてしまいます。そこで、残る数から考えると上記の累積値として求まります。

つまり、nの素因数が2,3,・・・,mの場合、素因数で割っていって残る数の個数φ(n)は、
 n*(1 - 1/2)*(1 - 1/3)*・・・*(1 - 1/m)
= n*(1/2)    *(2/3)    *・・・*(m-1)/m
となります。

そこで1からnまでの一つひとつの数字について、
素因数を求めて上記のようにφ(n)を計算し総和を求めてみましたが、
このロジックのまま実装すると1分ルールを守れませんでした。
素因数を何度も求める部分が重いようでした。

そこでさらにひと工夫です。
φ(n)をnの順に求めるのではなく素因数候補の方を順に反映していくようにしてみます。
つまり、mを素因数とするφ(i)すべてに対して、(m-1)/m倍を累積していきます。
素因数候補が2ならば、φ(2)から2刻みにφ(4), φ(6)・・・を途中計算として1/2倍します。
素因数候補が3ならば、φ(3)から3刻みにφ(6), φ(9)・・・を途中計算として2/3倍します。
素因数候補がmならば、φ(m)からm刻みにφ(m*2), ・・・を途中計算として(m-1)/m倍します。

1.関数f(n)
・φ(0)からφ(n)のリストLを返します。リストLのi番目の要素L[i]=φ(i)です。
・φ(n)の値のリストLとして、最初に、i番目の要素に値iをセットした値nまでのリスト[0, 1,..., n]を準備します。
・ループ変数iが素因数候補です。
 iを小さい値から順に処理したとき、
 L[i]==iならば、上記の(m-1)/mが反映されていないので素因数です。
 L[i]<iならば、すでにiより小さい値で(m-1)/mされてるので素因数ではありません。
・なお、(m-1)/m倍するとき、この分数を先に計算するとfloat型になります。
 int型のまま計算する方が処理時間が短くすむため、(m-1)倍と1/m倍は別々に計算しています。
・ループ変数jがφ(n)の引数nです。
・ループ処理がすべて終わったら、リストLの2番目以降の要素の和を求めます。


解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
303963552391