プロジェクトオイラー 問６９

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



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

n	互いに素な数	φ(n)	n/φ(n)
2	1	1	2
3	1,2	2	1.5
4	1,3	2	2
5	1,2,3,4	4	1.25
6	1,5	2	3
7	1,2,3,4,5,6	6	1.1666...
8	1,3,5,7	4	2
9	1,2,4,5,7,8	6	1.5
10	1,3,7,9	4	2.5

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


－－－－－

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


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

#最大公約数を求める（エラトステネスのふるい）
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未満の素数すべてを素因数として１つずつ持つ値となります。

また、n/φ(n) の分母の値の増え方は１ずつですが、
nに素数因子を１つ持つごとのφ(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

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

１．h(n)
　引数nが素数ならばその数を、素数以外ならば1を返します。
　問５８で作成した素数判定を修正し、素数自身または1を返却することで、
　呼び出し元で簡単に素数だけ掛け算できるようにしました。
　詳細は問５８を参照してください。
　
・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はループ変数の増分です。

２．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未満の掛け算累積値を返します。

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


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