2013年11月17日日曜日

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

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

問87「3つの素数のべき乗」
素数の2乗と素数の3乗と素数の4乗の和で表される最小の数は28である. 
50未満のこのような数は丁度4つある.

28 = 22 + 23 + 24
33 = 32 + 23 + 24
49 = 52 + 23 + 24
47 = 22 + 33 + 24

では, 50,000,000未満の数で, 素数の2乗と素数の3乗と素数の4乗の和で表される数は何個あるか?





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






私の回答例は以下の通りです。
def p(n):
	L = [0,0]+[1]*(n-1)
	for i in xrange(n):
		if i*i>n: break
		if not L[i]: continue
		for j in xrange(i+i, n+1, i): L[j] = 0
	return [i for i in xrange(n+1) if L[i]]
def f(n):
	P, L = p(int((n-2**3-2**4)**.5)), []
	for i in P:
		for j in P:
			for k in P:
				s = k**2+j**3+i**4
				if n<=s: break
				L.append(s)
	return len(set(L))

n = 50000000
m = f(n)
print m


1.関数p(n)
・n以下の素数リストを返します。問37を流用しました。

2.関数f(n)
・素数の2乗と素数の3乗と素数の4乗の和で表されるn未満の数の個数を返します。

・Pに素数、Lには素数の2乗と素数の3乗と素数の4乗の和を格納します。
 ただし、使用する最大の素数は、自身の2乗と2の3乗と2の4乗の和がn未満になる最大素数なので、nから2の3乗と2の4乗を差し引いた値の平方根の整数部分です。
 この使用可能な最大素数までの素数をリストPに格納します。

・素数リストLから、3重のfor文でi,j,kとして順番に1つずつ、重複を許して取り出し、それぞれを2乗、3乗、4乗して、その和を変数sとします。

・変数sがnを超えたら、それ以上のkのループは不要なのでkのfor文を終了します。
 j,kのループにも同様の判定を入れることは可能ですが、1分ルールをクリアしたのでこのままとします。
 

・return len(set(L))
 set文でリストLの要素の値をユニークにして、len関数で個数を数えて戻します。

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

2013年11月16日土曜日

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

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

問86「直方体の道筋」
下に示す直方体は寸法が6×5×3である.
この直方体の1つの頂点Sにクモがいる. また反対の頂点Fにはハエがいる. 
SからFまでの壁に沿って直線移動する最短ルートは図に示す通りで, この長さは10である.

この最短ルートの候補は3本あるが, 最短のものがいつも整数長さとは限らない.
さて, M×M×M以下の寸法の直方体について, 最短ルートが整数である直方体の数を考える. 
M=100のとき, 条件を満たす直方体は2060個ある. このM=100は個数が2000を超える最小のMである. なお, M=99のときは1975個である.

100万個を超える最小のMを求めよ.





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






私の回答例は以下の通りです。
def f(n):
	m = c = 0
	while c<=n:
		m += 1
		k, L = m*m, []
		for i in xrange(2, m*2+1):
			s = (i*i+k)**.5
			if s==int(s): L.append(i)
		if not L: continue
		for i in L:
			c += i/2
			if m<i: c -= (i-m-1)
	return m, c

n = 1000000
m, c = f(n)
print "M:", m
print "num of routes:", c

反対側の頂点までの最短ルートの候補は、展開図で考えると以下の3つです。
直方体の各辺の長さをa,b,cとして、
・直角を挟む2辺が「a」「bとcの和」の直角三角形の斜辺:a2+(b+c)2 = a2+b2+c2+2bc
・直角を挟む2辺が「b」「cとaの和」の直角三角形の斜辺:b2+(c+a)2 = a2+b2+c2+2ca
・直角を挟む2辺が「c」「aとbの和」の直角三角形の斜辺:c2+(a+b)2 = a2+b2+c2+2ab

a≦b≦cとすると、最短ルートをsとして以下のとおりです。
s2 = a2+b2+c2+2ab = (a+b)2+c2 ・・・(A)

1.関数f(n)
・直方体のうち、各辺の長さと反対の頂点までの最短距離が整数長さのm×m×m以下の寸法で、その個数がn個を超える最小のmを返します。

・カウントすべき直方体の個数を変数cにカウントアップします。mとcはゼロクリアしておきます。

・whileループはカウントすべき直方体の個数cが限界値nを超えたら終了します。

・数え上げる直方体の最長辺がmなので式(A)は、a≦b≦c≦mという条件が付きます。最長辺の長さcをmとして、c2を変数kとします。

・for i in xrange(2, m*2+1):
 a+bの値を変数iとして1ずつ増やしていきます。
 この変数iの最小値はa=b=1のときで2、最大値はa=b=mのときなのでm*2です。
 xrange関数は第1引数から(第2引数-1)の値を第3引数(初期値1)刻みで発生させます。

・s = (i*i+k)**.5
 if s==int(s): L.append(i)
 式(A)を実装して、最短ルート長さsを求め、これをint関数で小数点以下切り捨てても同じ値かということで整数判定します。このsが整数ならばリストLにためます。

・if not L: continue
 上記ロジックでリストLが空の場合は整数長さの最短ルートがないので、現在のmでの探索は終了です。

・for i in L:
 c += i/2
 if m<i: c -= (i-m-1)
 これまでのロジックでリストLには題意に合うa+bの値をためていますので、1つずつ取り出して、a,bが整数になる組み合わせの個数を変数cに累計していきます。
 リストLから取り出した値iが2つの整数の和である組み合わせは、
 (1, i-1), (2, i-2), ..., (i-2, 2), (i-1, 1)
 の(i-1)通りで、真ん中以降は同じ組み合わせなので、異なる組み合わせとしては半分の(i-1)/2個です。i-1が奇数の場合は切り捨てます。

 リストLから取り出した変数iがm以下ならば、aとbの組み合わせは全部採用ですが、リストLには最大m*2までの値があるので、a,bの片方がmを超えてしまう分を差し引きます。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
M: 1818
num of routes: 1000457

2013年11月10日日曜日

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

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

問85「長方形の数え上げ」
注意深く数えると, 横が3, 縦が2の長方形の格子には, 18個の長方形が含まれている.

ぴったり2,000,000個の長方形を含むような長方形の格子は存在しない. 一番近い解を持つような格子の面積を求めよ.






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






私の回答例は以下の通りです。
def f(n):
	xmax, sk, xk, yk = int((-1+(n*8+1)**0.5)/2), 1, 1, 1
	for x in xrange(1, xmax):
		t = x*(x+1)
		ys = int(((-t)+(t*t+16*t*n)**0.5)/2/t)
		for y in [ys, ys+1]:
			s = x*(x+1)*y*(y+1)/4
			if s==n: return x*y, x, y, s
			if abs(n-s)<abs(n-sk): sk, xk, yk = s, x, y
	return xk*yk, xk, yk, sk

n = 2000000
a, x, y, m = f(n)
print "area:", a
print "grid:", x,"x",y
print "num of rectangular:", m



格子が縦1個横x個の場合、含まれる長方形の数nは、1からnまでの和であり、
n = x(x+1)/2
これをxの二次方程式として解くと、
x2 + x - 2n = 0
x = (-1 ± √1-4x1x(-2n)) / 2
x>0より
x = (-1 + √1+8n) / 2 ・・・(A)

このxは、縦1個のときの横の個数で、横の個数の最大値xmaxとなります。
こうして、xのループを1からxmaxの値まで1つずつ変化させます。

縦x個横y個の格子に含まれる長方形の個数nは、
縦方向にx(x+1)/2、横方向にy(y+1)/2の積です。
n = x(x+1)/2 × y(y+1)/2
n = x(x+1)y(y+1)/4 ・・・(B)
ここで、t = x(x+1)と置くと、
n = ty(y+1)/4
n = (ty2 + ty)/4
ty2 + ty - 4n = 0

yの二次方程式として解くと、
y = (-t ± √(t2-4t(-4n))) / 2t
y>0より
y = (-t + √(t2+16nt)) / 2t
yは整数値なので、上記yの前後の整数、つまり上記yをint関数で切り捨てた値と、これに+1した値のいずれかが横の個数です。
上記の縦横の個数の候補の組み合わせで含まれる長方形の個数を計算し、ぴったりn個ならば終了、それ以外はnとの差の絶対値が小さいものをとっておきます。


1.関数f(n)
・n個の長方形を含む格子に一番近い解を持つ格子の面積、及びそのときの縦の個数、横の個数、長方形の数を返します。
 ただし、縦の個数は最も小さい値とします。

・xmax, sk, xk, yk = int((-1+(n*8+1)**0.5)/2), 1, 1, 1
 初期値設定します。縦個数の最大値xmaxとして上記式(A)を実装します。
 平方根はmathモジュールにsqrt関数もありますが、
 累乗演算子**で0.5乗しても同じです。

・ys = int(((-t)+(t*t+16*t*n)**0.5)/2/t)
 横個数の最大値ysとして上記式(C)を実装します。

・s = x*(x+1)*y*(y+1)/4
 含まれる長方形の数sとして上記式(B)を実装します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
area: 2772
grid: 36 x 77

num of rectangular: 1999998

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

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

問84「モノポリーの確率」 
モノポリーの標準的な盤面は以下である:

GO A1 CC1 A2 T1 R1 B1 CH1 B2 B3 JAIL
H2   C1
T2   U1
H1   C2
CH3   C3
R4   R2
G3   D1
CC3   CC2
G2   D2
G1   D3
G2J F3 U2 F2 F1 R3 E3 E2 CH2 E1 FP

各プレイヤーはGOのマスから開始し, 2個の6面サイコロを用いて時計回りに進む. 
他のルールが無いとすれば, 各マスに止まる確率は全て等しく, 2.5%である. 
しかし, G2J (Go To Jail), CC (Community Chest, 共同基金), CH (Chance, チャンス) のマスによってこの確率は変えられてしまう.

G2Jに止まる, または, CCやCHのマスに止まると引くカードのうちそれぞれ1枚によって, プレイヤーはJAILのマスに飛ばされる.
またプレイヤーが連続して3回ゾロ目を出すと, 3投目の結果のマスに進むのではなく, 直接JAILのマスに飛ばされる.
(筆者注:ぞろ目が連続して出る場合、連続が3回以内続く分の目の合計を進むのではなく、さいころを振る都度止まるマスの指示に従う)

ゲーム開始前にCCカードとCHカードはシャッフルされる.
プレイヤーがCCまたはCHマスに止まった場合, プレイヤーはCCカードまたはCHカードの山の一番上からカードを1枚引く. 
カードの指示に従ったのち, そのカードは山の一番下に戻される.
それぞれのカードは16枚あるが, 今回は問題を進み方に限定するので, 移動の指示があるカードのみを考える.
移動の指示が無いカードに関しては何もせずカードをそのまま山の下に戻す.
プレイヤーはそのままCC/CHマスに残るものとする.

Community Chest (16枚中2枚が移動カード) 
GOへ進め 
JAILへ進め
Chance (16枚中10枚が移動カード) 
GOへ進め 
JAILへ進め 
C1へ進め 
E3へ進め 
H2へ進め 
R1へ進め 
次のRへ進め (Rはrailway company, 鉄道会社の意) 
次のRへ進め 
次のUへ進め (Uはutility company, 公共事業会社の意) 
3マス戻れ
今回考えるのは, どのマスに止まりやすいかである. 
即ち, サイコロを投げた後に止まる確率である. 
より正確には, サイコロを1回振ってカードやマスによる移動を終えた後に各マスに止まる確率を求めたい. 
従って, G2Jに止まる確率は0であり, CHマスに止まる確率はその次に少ない(16枚中10枚が移動カードなので).
またJAILマスにたまたま止まることとJAILマスに送られることを区別しない.
またJAILマスから抜けるルール (自分のターンにゾロ目を2回出す) を無視する.
つまり必ず保釈金を払ってJAILマスから進むものとする.

GOマスを00とし番号を00-39と順番に振る. これにより各マスを2桁の数で表すことが出来る.
統計的には, 3つのマスに止まりやすいことを示せる.
JAIL (6.24%) = 10番目, E3 (3.18%) = 24番目, GO (3.09%) = 00番目である.
従ってもっとも止まりやすいマスを6桁で表せて102400となる.

2つの6面サイコロではなくて, 2つの4面サイコロを用いた場合の, もっとも止まりやすいマスを6桁で表せ.





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






私の回答例は以下の通りです。
def f(n, d):
	import random
	#Board List
	BL = "GO A1 CC1 A2 T1 R1 B1 CH1 B2 B3 JAIL C1 U1 C2 C3 R2 D1 CC2 D2 D3 \
		FP E1 CH2 E2 E3 R3 F1 F2 U2 F3 G2J G1 G2 CC3 G3 R4 CH3 H1 T2 H2".split()
	BD = dict([(s, i) for i,s in enumerate(BL)])

	#Next Railway company
	L = [i for i,s in enumerate(BL) if s[0]=="R" and s[1:].isdigit()]
	NxR = [i for i,j in zip(L, [0]+L[:-1]) for k in xrange(i-j)]
	NxR += [NxR[0]]*(len(BL)-len(NxR))

	#Next Utility company
	L = [i for i,s in enumerate(BL) if s[0]=="U" and s[1:].isdigit()]
	NxU = [i for i,j in zip(L, [0]+L[:-1]) for k in xrange(i-j)]
	NxU += [NxU[0]]*(len(BL)-len(NxU))

	#community chest
	cc = [i for i,s in enumerate(BL) if s[:2]=="CC" and s[-1].isdigit()]
	ccL = [BD["GO"], BD["JAIL"]]+[""]*14
	random.shuffle(ccL)

	#chance
	ch = [i for i,s in enumerate(BL) if s[:2]=="CH" and s[-1].isdigit()]
	chL = [BD["GO"], BD["JAIL"], BD["C1"], BD["E3"], BD["H2"], BD["R1"], \
			"NxR", "NxR", "NxU", "bk3"]+[""]*6
	random.shuffle(chL)

	iCD = cntcc = cntch = p = 0
	B = dict([i, 0] for i in xrange(len(BL)))
	for i in xrange(n):
		d1, d2 = random.randrange(1, d+1), random.randrange(1, d+1)
		p = (p+d1+d2)%len(BL)

		#consecutive dobles
		if d1==d2:
			iCD += 1
			if iCD>=3: iCD, p = 0, BD["JAIL"]
		else: iCD = 0

		while True:
			#G2J
			if p==BD["G2J"]: p = BD["JAIL"]; break
			#JAIL
			elif p==BD["JAIL"]: break
			#comunity chest
			elif p in cc:
				t = ccL[cntcc]
				p = {True:t, False:p}[isinstance(t, int)]
				cntcc = (cntcc+1)%len(ccL)
				if p not in ch: break
			#chance
			elif p in ch:
				t = chL[cntch]
				t = {"NxR":NxR[p], "NxU":NxU[p], "bk3":p-3}.get(t,t)
				if isinstance(t, int): p = t
				cntch = (cntch+1)%len(chL)
				if p not in cc: break
			#etc.
			else: break

		B[p] += 1

	L = [t for t in sorted(B.items(), key=lambda x:x[1], reverse=True)]
	r = "".join([str(L[i][0]).zfill(2) for i in xrange(3)])
	P = [str(L[i][1]*100.0/n)[:4] for i in xrange(3)]
	return L, r, P

L, r, (p0, p1, p2) = f(1000000, 4)
print r, str(p0)+"%", str(p1)+"%", str(p2)+"%"


共同基金やチャンスカードがあるので数学的に確率を求めるのは難しいと思い、
PC上で十分に多い回数を試行することにしました。
6面さいころで試行したところ、さいころを100万回振ると、結果として、問題文中の値との差異が0.05%程度以内に収まりました。
そこで、4面さいころを100万回振ることで解を求めました。

1.関数f(n,d)
・モノポリーでd面さいころをn回振ることを試行し、止まった回数の多いマスのベスト3を6桁表示した値とそのマスに止まった割合を返します。

・最初に固定データをリストや辞書にします。
 盤面位置リストBL、マス名称から盤面位置を取得できる辞書BD、
 各マスごとに「次の鉄道会社」の位置のリストNxR、
 各マスごとに「次の公共事業会社」の位置のリストNxU、
 共同基金の位置リストcc、共同基金カードリストccL、
 チャンスの位置リストch、チャンスカードリストchL

・import random
 pythonにある乱数に関するいろいろな関数を使えるようにしておきます。

・ramdom.shuffle(ccL)
 共同基金カードリストccLをシャッフルします。最初にカードをよく切り混ぜるイメージです。チャンスカードも同様によく切り混ぜます。

・カウンタなどを初期設定しておきます。
 ぞろ目(consecusive doubles)の連続回数iCD、共同基金カードの現在使用位置cntcc、チャンスカードの現在使用位置cntch、盤面の現在位置pはゼロクリアしておきます。
 マス位置別に止まった回数カウンターを用意し辞書Bとし、各カウンタはゼロクリアしておきます。

・for i in xrange(n):
 さいころを振る回数iのループです。n回振ります。

・d1, d2 = random.randrange(1, d+1), random.randrange(1, d+1)
 d1とd2それぞれがさいころの目です。
 randrange関数で1からdまで整数を乱数発生させ、さいころの目とします。

・ぞろ目判定では、ぞろ目のときにぞろ目連続回数iCDをカウントアップし、ぞろ目以外の場合必ずゼロクリアすることで連続3回か判定しています

・共同基金とチャンスでは、出るカードによって進む先が再び共同基金やチャンスのことがあるので、連続してカードを引けるように無限ループに入れます。
どちらもこれ以外が出ればループから抜けます。

・共同基金カードは使用した位置をcntchで保持して1つずつ加算することで、カードの山の上から順番に使って使用後はカードの山の一番下に戻すイメージを実現しました。チャンスカードも同様です。

・L = [t for t in sorted(B.items(), key=lambda x:x[1], reverse=True)]
 sorted関数で並び替えた値を取り出して内包表記でリストLにします。
 sorted関数の引数は以下のとおりです。
・B.items()は、マス位置別に止まった回数カウンター辞書Bから、キー(マス位置)と値(止まった回数)の組のタプルを要素として取り出したリストです。
・key=lambda x:x[1]
 lambdaは無名関数で引数xを受け取り、x[1]、配列位置1つまり2番目の要素を返します。ここでは辞書Bをリスト化したデータを受け取り、止まった回数を返します。これがsorted関数での並び替えるキーとなります。
・reverse=Trueでsorted関数での並び順は降順、ここでは止まった回数の多い順となります。

・r = "".join([str(L[i][0]).zfill(2) for i in xrange(3)])
 止まった回数の多いベスト3の6桁表示値を求めます。
 止まった回数順に並び替えたリストLから、ループiで配列位置0から2番目の要素を取得します。ここでは止まった回数の多いベスト3の要素が取得されます。
 L[i][0]は対象のマス位置の数字なので、str関数で文字列にして、zfill関数で0埋めで2桁でします。
 さらに、このベスト3のマス位置の2桁値を、join関数で区切り文字なしで文字列連結します。

・P = [str(L[i][1]*100.0/n)[:4] for i in xrange(3)]
 止まった回数の多い方からベスト3の回数のパーセントを求めて、順に格納したリストです。


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

101524 7.04%% 3.65% 3.29% ※ベスト3の割合はたまたま実行したときの値です。