プロジェクトオイラー 問１００

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

問１００「組み合わせの確率」
箱の中に15個の青い円盤と6個の赤い円盤からなる21個の色のついた円盤が入っていて, 無作為に2つ取り出すとき, 青い円盤2つを取り出す確率P(BB)は
P(BB) = (15/21) × (14/20) = 1/2
であることがわかる.

無作為に2つ取り出すとき, 青い円盤2つを取り出す確率がちょうど1/2となるような次の組み合わせは, 箱が85個の青い円盤と35個の赤い円盤からなるときである.

箱の中の円盤の合計が1012 = 1,000,000,000,000を超えるような最初の組み合わせを考える. 
箱の中の青い円盤の数を求めよ.






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






私の回答例は以下の通りです。

def f(n):
	x, y, i = 1, 1, 0
	while i<=n:
		x, y = x+y*2, x+y
		if x*y%2: i, j = (x+1)/2, (y+1)/2
	return i,j

n = 10**12
s, b = f(n)
print "blue : ", b
print "red  : ", s-b
print "sum  :", s

問題のまま実装するととても１分では処理が終わらず、プログラミングの前に数学的に解法を整理しておかないと解けないという、プロジェクトオイラーらしい問題です。

円盤総数をi個、青円盤数をj個とすると、定義から以下が成り立ちます。
j/i * (j-1)/(i-1) = 1/2
2j(j-1) = i(i-1)
i2-i - 2(j2-j) = 0
(i - 1/2)2-(1/4) - 2(j - 1/2)2 + (1/2) = 0
(i - 1/2)2 - 2(j - 1/2)2 + 1/4 = 0
4(i- 1/2)2 - 8(j - 1/2)2 + 1 = 0
(2i-1)2 - 2(2j-1)2 = -1  ...(a)
ここで、媒介変数xとyを
x = 2i-1, y = 2j-1 ...(b)
とおくと、
x2 -  2y2 = -1 ...(c)
これは、プロジェクトオイラーではおなじみのペル方程式で、
円盤総数i=(x+1)/2が、nを超えた最初の解のときの青円盤数j=(y+1)/2が求める値です。

問９４ より、ペル方程式の公式は以下の通り。
x2 - Dy2 = ±1　（Dは平方数以外。平方数ならば整数解が存在しない）
の解は、最小整数解を(x0, y0)、それ以降の解を順に、(x1, y1), (x2, y2)、...として、以下の連立方程式が成り立つ。
┌xn+1 = x0*xn + D*y0*yn
└yn+1 = y0*xn +   x0*yn ...(d)

最小整数解(x0, y0)は、式(c)から明らかで、x0=1, y0=1。また式(c)からD=2。
これらから連立方程式(d)は以下のとおり。
┌xn+1 = xn + 2yn
└yn+1 = xn + yn ...(e)
なお、式(a)より、xnとynは両方とも奇数。
これを実装します。

１．関数f(n)
・x, y, i = 1, 1, 0
　媒介変数xとyの初期値はともに1、総円盤数iの初期値は0としておきます。

・while i<=n:
　総円盤数iが指定値nを超えるまで処理を行います。

・x, y = x+y*2, x+y
　連立方程式(e)を実装します。

・if x*y%2: i, j = (x+1)/2, (y+1)/2
　媒介変数xとyの両方が奇数という判定は以下のとおりです。
　xとyの両方が奇数のときだけ、その積が奇数になるので、
　xとyの積を２で割った余りが１（論理判定ではTrue）になるかで判定します。
　そして、xとyが両方とも奇数のときに、式(b)から割り戻した値を
　円盤総数iと青円盤数jに退避しておきます。


解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
blue :  756872327473
red  :  313506783024
sum  : 1070379110497

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

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

問９９「最大のべき乗」
指数の形で表される2つの数, 例えば 211 と 37, の大小を調べることは難しくはない. 電卓を使えば, 211 = 2048 < 37 = 2187 であることが確かめられる.

しかし, 632382518061 > 519432525806 を確認することは非常に難しい (両者ともに300万桁以上になる).

各行に1組が書かれている1000個の組を含んだ22Kのテキストファイル base_exp.txt から, 最大の数が書かれている行の番号を求めよ.

注: ファイル中の最初の二行は上の例である.






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






私の回答例は以下の通りです。

# -*- coding:sjis -*-
def f(sIn):
	import math
	pIn = open(sIn)
	c = 0
	for i, r in enumerate(pIn.readlines()):
		a, b = r[:-1].split(",")
		v = int(b)*math.log10(int(a))
		if c<v: c, line, rec, digit = v, i+1, r[:-1], int(v)+1
	pIn.close()
	pIn = None
	return line, rec, digit
	
sIn = os.path.join(os.path.dirname(__file__),"p099_base_exp.txt")
line, rec, digit = f(sIn)
print "line :", line
print "rec  :", rec
print "digit:", digit

１．関数f(sIn)
・aのb乗をpとすると、
　p = ab
　両辺の10を底とする対数をとって、
　log(p) = b x log(a)
　なお、上記式の両辺の対数の底を10にすれば、
　この両辺の値を切り捨て＋１した値が桁数になります。
　この式を実装して最大値を判定します。
　
・for i, r in enumerate(pIn.readlines()):
　readline関数で1行ずつ読み込みます。

・a, b = r[:-1].split(",")
　readline関数で末尾の改行コードは取り除かないので、末尾の1バイトを削除し、半角カンマで基数と指数に分けます。

・v = int(b)*math.log10(int(a))
　vは、aのb乗の値の、10を底とする対数値です。

・if c<v: c, line, rec, digit = v, i+1, r[:-1], int(v)+1
　対数値vの最大値をcにためます。
　iは0から始まる行カウンタなので、行番号はi+1です。
　readline関数で末尾の改行コードは取り除かないので、行内容recは末尾の1バイトを削除します。
　digitは桁数です。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
line : 709
rec  : 895447,504922
digit: 3005316

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

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

問９８「アナグラム平方数」
CARE という単語の各文字をそれぞれ 1, 2, 9, 6 に置き換えることによって, 
平方数 1296 = 362 ができる. 注目すべきことに, 同じ数字の置換をつかうことにより, アナグラムの RACE も平方数 9216 = 962 をつくる. 
CARE (と RACE) を平方アナグラム単語対と呼ぼう. 先頭のゼロは許されず, 
異なる文字が同じ数字をもつこともないとする.

約 2,000 個の一般的な英単語を含む 16K のテキストファイルwords.txt 

(右クリックして "名前をつけてリンク先を保存") を用いて, 
平方アナグラム単語対をすべて求めよ (回文となる単語はそれ自身のアナグラムとはみなさない).

そのような対のメンバーから作られる最大の平方数は何か?

注: 作られるアナグラムは, すべて与えられたテキストファイルに含まれている.






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






私の回答例は以下の通りです。

def f(sIn):
	import itertools
	W = {}
	pIn = open(sIn)
	for s in pIn.read().replace('"', '').split(","):
		t = "".join(sorted(s))
		W[t] =W.get(t, [])+[s]
	pIn.close()
	pIn = None
	
	for k,v in W.items():
		if len(v)<2: W.pop(k)
		
	N, m = {}, max([len(s) for s in W])
	for i in xrange(int((10**m)**.5)+1):
		s = str(i*i)
		t = "".join(sorted(s))
		N[t] = N.get(t, [])+[s]
	
	for k,v in N.items():
		if len(v)<2: N.pop(k)
		
	nm = 0
	for s in W:
		for wp in itertools.permutations(W[s], 2):
			
			for t in N:
				if len(s)!=len(t): continue
				if len(set(s))!=len(set(t)): continue
				
				for np in itertools.permutations(N[t], 2):
					p = dict([(u, v) for (u, v) in zip(wp[0], np[0])])
					if "".join([p[u] for u in wp[1]])!=np[1]: continue
					i = int(max(np))
					if nm<i: nm, wpm, npm = i, wp, np
					
	return nm,wpm,npm,int(nm**.5)
	
sIn = os.path.join(os.path.dirname(__file__), "p098_words.txt")
s = f(sIn)
print s

１．関数f(sIn)
・単語ファイルsInを読み込んで、問題に合う最大平方数nm、その場合の単語ペアwpmと二乗値ペアnpm、最大平方数の＋の平方根を返します。

・単語辞書Wとして、キーにソート済み構成文字、値にそのアナグラム単語のリストを持つ辞書を作成します。
　アナグラムにならない単語、つまり単語辞書で値リストの要素数が２つ未満の場合は単語辞書からpop関数で削除します。
・同様に二乗値の方も、
　二乗値辞書Nとして、キーにソート済み構成数字、値にそのアナグラム二乗値のリストを持つ辞書を作成します。
　ただし、候補単語の最大文字列長mから最大桁数を計算し、二乗値辞書にはその桁数以下の値を対象にします。
　またアナグラムにならない二乗値、つまり二乗値辞書で値リストの要素数が２つ未満の場合は二乗値辞書からpop関数で削除します。

・単語辞書の値から、アナグラム単語ペアを順列（順番が変われば別扱いとした全組合せ）で作り出します。
　アナグラム単語が２つならば2P2=2通り、３つならば3P2=6通りのペアを作り出します。
・同様に二乗値の方も、アナグラム二乗値ペアを作り出します。
・単語と二乗値とで、桁数と構成文字種類数が同じ場合に比較していきます。
　桁数はlen関数で取得し、構成文字種類数はset関数で構成文字をユニークにしてからlen関数で個数を取得します。

　
・単語ペアwpと二乗値ペアnpのそれぞれ1つ目の値から、文字と数字の一字対応辞書pを作成します。
　この一字対応辞書で単語ペアの2つ目を翻訳して二乗値ペアの2つ目になるかをチェックします。
　チェックNGの場合は次の候補値へ進みます。
　単語ペアと辞書ペアを順列で作り出しているのは、一字対応辞書を作成元をペア1つ目に固定して、
　そのチェック対象をペア2つ目に固定しているためです。
　順番無関係の組合せcombinationではなく、順番が関係する順列permutationsを使用するのはこのためです。


解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
(18769, ('BOARD', 'BROAD'), ('17689', '18769'), 137)