2013年1月30日水曜日

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

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

問71「真既約分数を昇順に並べ上げよ」
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
3/7のすぐ左の分数は2/5である.

d ≦ 1,000,000について真既約分数を大きさ順に並べたとき, 
3/7のすぐ左の分数の分子を求めよ.

-----


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





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

def gcd(a, b):
	if a<b: a, b = b, a
	if b: return gcd(b, a%b)
	else: return a

def f(n, a, b):
	L = [0,0,1]
	s = 1.0*a/b
	for i in xrange(1, n+1):
		j = int(i*s)
		if not j: continue
		while gcd(i, j)>1:
			j -= 1
		u = s - 1.0*j/i
		if 0<u<L[2]: L = [j, i, u]
	return L

n, a, b = 1000000, 3, 7
print f(n, a, b)


1.gcd(a,b)
 「ユークリッドの互除法」を使用して最大公約数を返します。
 説明は問5を参照してください。

2.g(n,a,b)
 n以下の分母を持つ真規約分数を大きさ順に並べたときに、
 a/bのすぐ左の分数についての、分子、分母、その分数のa/bと差を返します。
 
 まず、基準値となるa/bを計算しsとします。
 1.0を掛けているのは整数型ではなく浮動小数点型で計算するためです。
 float関数で型変換するよりも速いです。
 ループ変数iは求める分数の分母で、jは分子です。
 分子は、分母と基準値を掛け、int関数で整数部をとることで算出します。
 jが0の場合、分数にならなので分母を1つ進めます。
 また、j/iが真既約分数にならない場合、真既約分数になるまでjを1つずつ小さい値にします。
 真既約分数の判定は、分母子の最大公約数が1であることとしました。
 こうして分数j/iが定まったところで、基準値との差分を求めuとします。
 uがここまで処理したうちの最小値ならば、分子j、分母i、差分uをリストLの値と差し替えてためます。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
[428570, 999997, 1.4285757138354782e-07]

2013年1月26日土曜日

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

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


問70「φ(n) が n の置換となる n を調べ上げよ」
オイラーのトーティエント関数 φ(n) (ファイ関数とも呼ばれる) とは, 
n 未満の正の整数で n と互いに素なものの個数を表す. 
例えば, 1, 2, 4, 5, 7, 8 は9未満で9と互いに素であるので, φ(9) = 6 となる. 
1 は全ての正の整数と互いに素であるとみなされる. よって φ(1) = 1 である.

面白いことに, φ(87109)=79180 であり, 87109は79180を置換したものとなっている.

1 < n < 10 7 で φ(n) が n を置換したものになっているもののうち, 
n/φ(n) が最小となる n を求めよ.

-----


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






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

def p(n): L = [0,0]+[1]*(n-1) i = 2 while i*i<=n: while not L[i]: i += 1 for j in xrange(i+i, n+1, i): L[j] = 0 i += 1 return [i for i in xrange(n+1) if L[i]] def f(n): N = [n, n, []] P = p(n//2) w = int(n**(0.5)) for i, u in enumerate(P): if w<u: break z = n//u for j in xrange(len(P)): if z<P[j]: break for v in P[j: i: -1]: s, t = u*v, (u-1)*(v-1) if sorted(str(s))==sorted(str(t)): a = 1.0*s/t if a<N[1]: N = [s, a, [u,v]] break return N n = 10000000 print f(n)


問69の発展問題ということで、同様にはなく、数学的に別の表現ができないかを考えます。
まず、問69とは逆に、n/φ(n)を最小にするためには、nをなるべく小さくしφ(n)はなるべく大きくします。
また、問69の検討を考慮すると、
・nよりもφ(n)が大きく効いてくるので、φ(n)を大きくする。
・φ(n)の最大値は、nが素数のときのφ(n)=n-1なので、なるべく大きな素数を持つ。

ここで、nとφ(n)が置換の関係との条件から、素数では問題に合わないことがわかります。
それは、nとn-1は必ず1文字は異なり、置換の関係にならないからです。

そこで素数2つの積で、φ(n)との置換の関係を満たすnということになります。
それは、nが複数の素数の積ならば、その因数である素数の倍数以外には割りきれないため、φ(n)が大きい値のままですみます。
また、この場合の因数として使用する素数の件数なるべく少ない方がよいということになります。

1.関数p(n)
・n以下の素数リストを返します。問37と同じです。

2.関数f(n)
 n未満の整数で、問題に合うn、n/φ(n)、nを構成する素数2つのリストを返します。
 候補の素数の範囲は、2つの素数の積がn未満ということから、n//2以下で十分です。
 素数候補のリストから、
 ループ変数uで小さい方から1つずつ取り出し、
 ループ変数vで大きい方から1つずつ取り出します。
 素数uと素数vの積とそのφ値をs,tとします。
 置換の関係のチェックはs,tを1文字ずつ取り出してそのソート状態が一致することとします。
 リストNに関係情報をためて、n/φ(n)の小さい値がきたら入替えます。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
[8319823, 1.0007090511248113, [2339, 3557]]

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