2013年1月20日日曜日

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

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



問69「n/φ(n) が最大となる n の値を求めよ」
オイラーのトーティエント関数, φ(n) [時々ファイ関数とも呼ばれる]は,
n と互いに素な n 未満の数の数を定める.
たとえば, 1, 2, 4, 5, 7, そして8はみな9未満で9と互いに素であり, φ(9)=6.


n互いに素な数φ(n)n/φ(n)
2112
31,221.5
41,322
51,2,3,441.25
61,523
71,2,3,4,5,661.1666...
81,3,5,742
91,2,4,5,7,861.5
101,3,7,942.5


n ≦ 10 では n/φ(n) の最大値は n=6 であることがわかる.
n ≦ 1,000,000で n/φ(n) が最大となる値を見つけよ.


-----

まず、問題文のとおりにプログラミングしてみました。
n未満の整数を1つずつnとともに「エラトステネスのふるい」にかけることで
その最大公約数が1ならば互いに素と判定してカウントし、n/φ(n)を計算してみました。
でもこれではn=1,000,000どころか、n=10,000でさえ1分を超え、1分ルールを守れません。


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



#最大公約数を求める(エラトステネスのふるい)
def gcd(a, b):
	if a<b: a, b = b, a
	if b: return gcd(b, a%b)
	else: return a
#n未満でnと互いに素な数の個数
def g(n):
	if n: L = [i for i in xrange(n) if gcd(n, i)<=1]
	else: L = [0]
	return len(L)
def f(n):
	L = [float(i)/g(i) for i in xrange(n+1)]
	m = max(L)
	return L.index(m), m

問題文のまま実装したのではロジック改良でクリアできるレベルではないので、
プログラミング以前に、まずこの問題を数学的に別の表現に言い換えることを考えました。

n/φ(n)を最大値にするためには、nをなるべく大きくしφ(n)はなるべく小さします。

問題文中の表をよくみると、まずnが素数以外であることがわかります。
nが素数の場合、1以上n未満の整数はすべて互いに素なので
φ(n)=n-1になり、他の値と比べてnに対するφ(n)の値が明らかに大きすぎます。

さらに、nが偶数のときにφ(n)の値が小さいことがわかります。
理由は、nが偶数ならば、n未満の偶数すべてnとは互いに素にならないためです。
同様にnが3の倍数ならば、n未満の3の倍数はすべてnと互いに素になりません。
つまり、pを素数とすると、
nがpの倍数ならば、n未満のpの倍数はすべてnと互いに素にならずφ(n)にカウントされません。
また、nの素因数にpを持てば、pの倍数は因数として持たなくてよいです。
つまり、φ(n)の最小値となるnは、
n未満の素数すべてを素因数として1つずつ持つ値となります。

また、n/φ(n) の分母の値の増え方は1ずつですが、
nに素数因子を1つ持つごとのφ(n)の減り方は明らかに1個よりみ大きいので、
分子φ(n)の減少の方が分母nの増加よりも大きく効いてきます。

上記から、φ(n)が最小値となるnは、
素数を小さい順に掛け、その積がn未満の最大値になるような値となります。

ちなみに、問題文のn≦10では、素数は、2, 3, 5, 7,・・・なので、
2×3=6で、n=6のときにn/φ(n) が最大値になります。



私の解答例は以下です。畳んでいます。
def h(n):
	if n%2 or n==2:
		for i in xrange(3, n, 2):
			if i*i>n: return n
			if not n%i: return 1
		return n
	return 1
 
def f(n):
	r = 1
	for i in xrange(n):
		r *= h(i)
		if r>n: return r//i



2個の関数を使っています。

1.h(n)
 引数nが素数ならばその数を、素数以外ならば1を返します。
 問58で作成した素数判定を修正し、素数自身または1を返却することで、
 呼び出し元で簡単に素数だけ掛け算できるようにしました。
 詳細は問58を参照してください。
 
・if n%2 or n==2:
 「n%2」はnを2で割った余りで、論理的にTrueになるのは1の場合、つまり奇数です。
 if文の中は、nが2か奇数の場合の処理です。

・for i in xrange(3, n, 2):
 ループ変数iは3以上n未満の奇数です。2はループ変数の増分です。

2.f(n)
・r = 1
 求める数をクリアします。掛け算していくので初期値は1です。

・for i in xrange(n):
 ループ変数iは0以上n未満の整数です。

・r *= h(i)
 関数h(i)の戻り値を累積して掛け算していきます。

・if r>n: return r//i
 掛け算累積値がnを超えたら、最後に掛けた数で割って
 n未満の掛け算累積値を返します。

3.関数の外
・1,000,000を関数f(n)に渡し、戻り値を表示します。


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

0 件のコメント:

コメントを投稿