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

プロジェクトオイラーの問題を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)

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

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

問９７「大きな非メルセンヌ素数」
100万桁を超える初めての素数は1999年に発見された. これはメルセンヌ素数であり, 26972593-1 である. 実際, 2,098,960桁ある. 
それ以降も, より多くの桁になるメルセンヌ素数 (2p-1の形の数) が他にも発見されている.

しかし, 2004年に, 非常に大きな非メルセンヌ素数が発見された. 
これは2,357,207桁の数であり, 28433×27830457+1である.

この素数の末尾10桁を答えよ.






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






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

def f(a, n, p, b, m):
	k, t, c = n, 10**m, 1
	for s in str(p)[::-1]:
		c *= k**(int(s))%t
		k = k**10%t
	return (a*c+b)%t

def f2(a, p, b, m): return ((a<<p)+b)%(10**m)

#（ａ×（ｎのｐ乗）＋ｂ）の末尾m桁を返す
a, n, p ,b ,m =28433, 2, 7830457, 1,10
s = f(a, n, p ,b ,m)
print s

#（ａ×（２のｐ乗）＋ｂ）の末尾m桁を返す
a, p ,b ,m =28433, 7830457, 1,10
s = f2(a, p ,b ,m)
print s

１．関数f(a, n, p, b, m)
・a×np＋b　の末尾m桁を返します。

・for s in str(p)[::-1]:
　数値pをstr関数で文字型にして、文字のスライス機能の[開始:終了:刻み]での刻みを-1にすることで右から１文字ずつ、つまり１の位の数字から順に１つずつ取り出し、ループ変数sに設定します。

・c *= k**(int(s))%t
　k = k**10%t
　変数kは、ループ変数sの桁位置に相当するnの累乗値を設定します。
　初期値はnで、ループが回るにつれ、nの10回累乗、100回累乗、...と直前の累乗値を10回ずつ累乗していきます。
　変数cに、kのs乗をかけて、ためていきます。
　例えば、
　2の345乗ならば、(2の5乗) x ((2の10乗))の4乗) x ((2の100乗)の3乗)です。
　このときの「2の10乗」「2の100乗」が変数kで、各桁の数字が変数sです。

２．（別解）関数f2(a, p ,b ,m)
・a×2p＋b　の末尾m桁を返します。
　関数の中の累乗計算を、2の累乗に固定します。

・a<<p
　数値aをpビット左にシフトすることで、2進数でp桁繰り上がります。
　これは、数値aに対して2をp回かけることと同じです。

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

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

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

問９６「数独」
"数独"(日本語で数字を配置するという意味)とは人気があるパズルの名前である.
起源は不明だが, その評判はラテン語で"Squares"と呼ばれる同様な, 
そしてはるかに難しいパズルを考案したレオンハルト・オイラーの貢献によるものに違いない. 
しかしながら, "数独"パズルの目的はそれぞれの行, 列が3×3の枠を含む9×9の格子の空白（もしくは0）をそれぞれ1から9の数字で置き換えることである. 
下に, 一般的なパズルの開始状態とその解答の例がある.

うまく作られている"数独"パズルは, 選択肢を消去するために"仮定とテスト"方式を用いる必要があるかもしれないが, ただ一つの解を持ち, 論理によって解くことができる(これについては様々な意見がある).

探索の複雑さがパズルの難易度を決定する; 上に挙げた例は, 単純で直接的な推論によって解く事ができるため, 簡単であると考えられる.

6kバイトのテキストファイルsudoku.txt(右クリックで,"名前をつけてリンク先を保存") にはただ一つの解を持つ, 様々な難易度の50の"数独"パズルが含まれている(上の例題はこのファイルにおける最初のパズルである).

50すべてのパズルを解き, それぞれの解答の左上隅にある3桁の数の合計が求めよ; 
例えば483は上の解答例の左上隅の3桁の数である.






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







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

def d0(d):
	for x,y in Lxy():
		if not d[x,y]: d[x,y] = list(xrange(1,10))
	return d

def d1(d):
	for x,y in Lt(d):
		d[x, y] = sorted(set(d[x, y])-set(Ka(d,x,y)))
	return d

def d2(d):
	for x,y in Lt(d):
		for L in [Ltx(d,y), Lty(d,x), Ltm(d,x,y)]:
			M = [s for i,j in L for s in d[i,j]]
			N = [s for s in d[x,y] if M.count(s)==1]
			if N: d[x,y] = N
	return d

def d3(d):
	for n in xrange(1, 10):
		setn = set([n])
		for x,y in Lm3():
			N = [(i,j) for i,j in Ltm(d,x,y) if (n in d[i,j])]
			if len(N)!=2: continue
			if N[0][0]==N[1][0]:
				for i,j in Lty(d, N[0][0]):
					if (i,j) not in N: d[i,j] = sorted(set(d[i,j])-setn)
			if N[0][1]==N[1][1]:
				for i,j in Ltx(d, N[0][1]):
					if (i,j) not in N: d[i,j] = sorted(set(d[i,j])-setn)
	return d

def d4(d):
	for n in xrange(1, 10):
		setn = set([n])
		for s in xrange(9):
			L = [j for i,j in Lty(d,s) if (n in d[i,j])]
			if len(L)!=2: continue
			for i in sorted(set(xrange(9))-set([s])):
				N = [j for i,j in Lty(d,i) if (n in d[i,j])]
				if len(N)!=2: continue
				if sorted(L)!=sorted(N): continue
				for j in L:
					for k,j in Ltx(d,j):
						if k not in [s,i]: d[k,j] = sorted(set(d[k,j])-setn)
					
		for s in xrange(9):
			L = [i for i,j in Ltx(d,s) if (n in d[i,j])]
			if len(L)!=2: continue
			for j in sorted(set(xrange(9))-set([s])):
				N = [i for i,j in Ltx(d,j) if (n in d[i,j])]
				if len(N)!=2: continue
				if sorted(L)!=sorted(N): continue
				for i in L:
					for i,k in Lty(d,i):
						if k not in [s, j]: d[i,k] = sorted(set(d[i,k])-setn)
						
	return d

def d5(d):
	for x in xrange(9): d = f5(d, Lty(d, x), Vty(d, x))
	for y in xrange(9): d = f5(d, Ltx(d, y), Vtx(d, y))
	for x,y in Lm3(): d = f5(d, Ltm(d, x, y), Vtm(d, x, y))
	return d

def f5(d, Li, Lv):
	L = [s for s in Lv if len(s)==2]
	M = sorted(set([tuple(s) for s in L if Lv.count(s)==2]))
	M = [list(s) for s in M]
	for L in M:
		for n in L:
			setn = set([n])
			for i, j in Li:
				if d[i, j]!=L: d[i, j] = sorted(set(d[i, j])-setn)
	return d

def Lxy(): return [(x, y) for x in xrange(9) for y in xrange(9)]
def Lm(x, y): return [(x//3*3+i, y//3*3+j) for j in xrange(3) for i in xrange(3)]
def Lm3(): return [(x, y) for x in xrange(0,9,3) for y in xrange(0,9,3)]

def Lt(d): return [(x, y) for x,y in Lxy() if isinstance(d[x,y], list)]
def Ltx(d,y): return [(x, y) for x in xrange(9) if isinstance(d[x,y], list)]
def Lty(d,x): return [(x, y) for y in xrange(9) if isinstance(d[x,y], list)]
def Ltm(d,x,y): return [(x, y) for x, y in Lm(x,y) if isinstance(d[x,y], list)]

def Vx(d, y): return [d[x, y] for x in xrange(9)]
def Vy(d, x): return [d[x, y] for y in xrange(9)]
def Vm(d, x, y): return [d[x, y] for x,y in Lm(x, y)]
def Va(d, x, y): return Vx(d, y)+Vy(d, x)+Vm(d, x, y)

def Vtx(d, y): return [s for s in Vx(d, y) if isinstance(s, list)]
def Vty(d, x): return [s for s in Vy(d, x) if isinstance(s, list)]
def Vtm(d, x, y): return [s for s in Vm(d, x, y) if isinstance(s, list)]

def Kx(d, y): return [s for s in Vx(d, y) if isinstance(s, int)]
def Ky(d, x): return [s for s in Vy(d, x) if isinstance(s, int)]
def Km(d, x, y): return [s for s in Vm(d, x, y) if isinstance(s, int)]
def Ka(d, x, y): return sorted(set(Kx(d, y)+Ky(d, x)+Km(d, x, y)))

def chk(d):
	N = list(xrange(1, 10))
	for i in xrange(9):
		if sorted(Vx(d,i))!=N: return False
		if sorted(Vy(d,i))!=N: return False
	for x,y in Lm3():
		if sorted(Vm(d,x,y))!=N: return False
	return True

def f(sIn):
	pIn = open(sIn)
	L = pIn.read().split("\n")
	c = 0
	for i in xrange(0, len(L), 10):
		print "-------\n",L[i]
		d = {}
		for x,y in Lxy(): d[x,y] = int(L[i+y+1][x])
		d = d0(d)
		while not chk(d):
			dd = d.copy()
			for i in xrange(1, 6):
				for x, y in Lt(d):
					if len(d[x,y])==1: d[x,y] = d[x,y][0]
				d = eval("d"+str(i)+"(d)")
				if dd!=d: break
				
		for y in xrange(9): print y,[d[(x,y)] for x in xrange(9)]
		c += d[0,0]*100+d[1,0]*10+d[2,0]
		
	pIn.close()
	pIn = None
	return c
	
if __name__=="__main__":
	import os, sys
	sIn = os.path.join(os.path.dirname(__file__), "problem096_sudoku.txt")
	s = f(sIn)
	print s

数独を辞書dに格納します。
辞書dのキーは(0,0)-(8,8)の81マスの位置を示す座標タプルで、
値は候補値の配列（リスト型）、または確定値（整数型）です。

１．関数f(sIn)
・問題ファイルのパスを受け取り、リストLに各行を要素として全部読み込みます。
・リストLをヘッダを含めた10行ごとに取り出し、
　左上(0,0)-右下(8,8)の81マスの位置を示す座標タプルをキーとする辞書d
　に、値を１つずつ格納しておきます。
・最初に関数d0()で未設定値0を候補リストまたは確定値に変換しておきます。
・次に、関数chk()で数独辞書dが完成したかチェックOKになるまで、関数d1()から関数d5()を順に行います。
　ただし、関数d1()から関数d5()はこの順に基本的な操作から順に並べているので、マスの値が１つでも変化した時点で、改めて関数d1()から再処理します。
・数独が完成したら、左上３マスの値を３桁の数値にして、変数cに累積します。

２．関数d0():
・数独辞書dの未設定値0のマスに[1～9]を候補リストとして設定します。

３．関数d1()
・数独全体の候補リストについて、候補リストから関係項目に存在する確定値を削除します。

４．関数d2()
・行内、列内、3x3メッシュ内もそれぞれで、その全候補値の中に１つだけ存在する値があれば確定値とします。

５．関数d3()
・3x3メッシュ内に２つだけある候補値の位置が一列に並んでいるとき、
どちらかが確定値（例の○か◎）なので、
その並んでいる行内、または列内の他の候補からその値（例の×）を削除します。

６．関数d4()
・ある候補値が２つだけ存在する行が２つあり、しかもその列位置が同じ場合、
　これら４つの値のうち、いずれかの対角位置の候補値が確定値（例の○か◎）になるため、この２つの列内の候補リストから、この候補値（例の×）を削除します。
　列についても上記と同様に処理します。

７．関数d5()
・行内に、候補値が２つだけの同一の候補リストが２つのマスにある場合、
どちらの候補値もどちらかのマスのそれぞれの確定値（下図の○か◎）になるため、このマスを含む、列内、また3x3メッシュがあればそれも含めた候補リストからこの２つの値（例の×）を削除します。
列内も同様です。

８．その他
部品となる関数を以下の通りです。
・Lxy():全マスの項目番号リスト
・Lm(x, y)：数独辞書dのx列y行が含まれる3x3メッシュの項目番号リスト
・Lm3()   ：3x3メッシュの左上マスの項目番号リスト

・Lt(d)   ：数独辞書dの処理中の項目番号リスト
・Ltx(d,y)：数独辞書dのy行目の処理中の項目番号リスト（横方向）
・Lty(d,x)　：独辞書dのx列目の処理中の項目番号リスト（縦方向）
・Ltm(d,x,y)：数独辞書dのx列y行が含まれる3x3メッシュの処理中の項目番号リスト

・Vx(d, y)　　：数独辞書dのy行目の値リスト（横方向）
・Vy(d, x)　　：数独辞書dのx列目の値リスト（縦方向）
・Vm(d, x, y)：数独辞書dのx列y行目が含まれる3x3メッシュの値リスト
・Va(d, x, y)：数独辞書dのx列y行に関係する値リスト

・Vtx(d, y)　　：数独辞書dのy行目の値リスト（横方向）
・Vty(d, x)　　：数独辞書dのx列目の値リスト（縦方向）
・Vtm(d, x, y)：数独辞書dのx列y行目が含まれる3x3メッシュの値リスト

・Kx(d, y)　：数独辞書dのy行目の確定値リスト（横方向）
・Ky(d, x)　：数独辞書dのx列目の確定値リスト（縦方向）
・Km(d, x, y)：数独辞書dのx列y行目が含まれる3x3メッシュの確定値リスト
・Ka(d, x, y)：数独辞書dのx列y行に関係する確定値リスト

・chk(d):#数独辞書dの縦、横、3x3メッシュの全てが[1～9]であるかの最終チェック


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

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

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

問９５「友愛鎖」
ある数の真の約数とは, それ自身を除く約数すべてである. 例えば, 28 の真の約数は 1, 2, 4, 7, 14 である. これらの約数の和は 28 に等しいため, これを完全数と呼ぶ.
面白いことに, 220 の真の約数の和は 284 で, 284 の真の約数の和は 220 となっており, 二つの数が鎖をなしている. このため, 220 と 284 は友愛数と呼ばれる.

さらに長い鎖はあまり知られていないだろう. 例えば, 12496 から始めると, 5 つの数の鎖をなす.
12496 → 14288 → 15472 → 14536 → 14264 (→ 12496 → ...)
この鎖は出発点に戻っているため, 友愛鎖と呼ばれる.

いずれの要素も 1,000,000 を超えない最長の友愛鎖の最小のメンバーを求めよ.






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






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

def d(n):
	c, m = 1, int(n**.5)
	for i in xrange(2, m+1):
		if not n%i:
<			c += i
<			j = n//i
<			if i<j: c += j
	return c

def f(n):
	L, N = [d(i) for i in xrange(n+1)], []
	for i in xrange(n+1):
		j, M = L[i], []
		while (j not in M) and (j<=n):
<			M.append(j)
<			j = L[j]
		if M and i==M[-1] and len(N)<len(M): N = [i]+M[:-1]
	return len(N), min(N), N

n = 1000000
n,m,L = f(n)
print "min:",m
print "num:",n
print "chain:",L

１．関数d(n)
・関数d(n)はnの友愛数を返します。問２１を改良しました。
・数値nの「真の約数」の上限は、その定義からnの正の平方根です。
　n**(0.5)　が、0.5乗ということで、nの正の平方根です。
・変数cに真の約数を累積していきます。
　//演算は割り算の商（整数）で、%演算は割り算の余りです。
　変数iのforループで、nがiで割り切れるとき、変数cにiを累積します。
　割り切れるときは余り0で、%演算の結果がFalseになるので、not演算子でこのような場合を捉えます。
　またこのときの商j=n//iも約数なので、i<jならこれも変数cに累計します。
　nが平方数でiがその平方根の場合、i=jなのでiだけ累積しjは累積しません。

２．関数f(n)
・関数f(n)は、n以下の友愛数リストを取得し、その友愛数をたどり、最長の友愛鎖を返します。
・リストLは、その位置のオフセット値の友愛数を持つリストです。例）L[220]=284
・リストNには、最長の友愛鎖を格納します。
・ループ変数iは、友愛鎖をたどる最初の数値です。
・変数jは次の友愛数が格納されます。
・リストMには、数値iから始まる友愛鎖を格納します。
・while文のループ条件は、
　友愛数jが友愛鎖リストMになく、上限値以下であることとしています。
・友愛鎖リストMを最長友愛鎖リストNに保存する条件は、
　リストMに値があり、友愛鎖たどりの最初の数値iが友愛鎖の最後の値と一致し、
　リストNよりも長い（最長の友愛鎖）こととしています。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
min: 14316
num: 28
chain: [19116, 31704, 47616, 83328, 177792, 295488, 629072, 589786, 294896, 358336, 418904, 366556, 274924, 275444, 243760, 376736, 381028, 285778, 152990, 122410, 97946, 48976, 45946, 22976, 22744, 19916, 17716, 14316]

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

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

問９４「擬正三角形」
一辺の長さが整数の正三角形は面積が整数にならないことを示すのは簡単である. 
しかし, 5-5-6の辺を持つ殆ど正三角形に近い擬正三角形 (almost equilateral triangle) は面積が12で整数である.
以降, 二等辺三角形で, 3つめの辺の長さが他と1つしか違わないもの (5-5-6, 5-5-4等) を, 擬正三角形と呼ぶ.

さて, 周囲の長さが1,000,000,000以下の面積が整数になる擬正三角形を考え, その周囲の長さの総和を求めよ.






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






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

def f0(n):
	s, c = n//3, 0
	d = {}; dd = {}
	for i in xrange(1, s):
		ii = i*i
		d[i] = ii; dd[ii] = 0
	for i in xrange(3, s, 2):
		t = i//2; u = d[i]
		if (u-d[t]  ) in dd: c += i*3-1
		if (u-d[t+1]) in dd: c += i*3+1
	return c

n = 10**9
print f0(n)

１．関数f(n)
・関数f(n)は周囲の長さがn以下で面積が整数になる擬正三角形の周囲の長さの和を返します。
・変数sは１辺の最大値、変数cは周囲の長さの累積です。
・擬正三角形の二等辺に挟まれた角からもう１辺に垂線hを下ろすと、
　垂線の足が直角である、同じ面積の直角三角形２つに分割できます。
　（以下、半分直角三角形と呼びます）
　ここで半分直角三角形の３辺は以下のように置きます。
　h：擬正三角形の二等辺で挟まれた角からの垂線。
　a：擬正三角形の二等辺であり、半分直角三角形の斜辺。
　b：hと直角を構成する。擬正三角形の底辺の半分の長さ。

・最初のforループで二乗値の準備をします。
　辞書dには二乗値を格納し、辞書ddには整数二乗値をキーにその正の平方根を格納します。

・２つ目のforループは、ループ変数iが擬正三角形の二等辺の長さです。
　これは半分直角三角形の斜辺でもあり、問７５より、３辺が整数の直角三角形は、
> 直角を挟む辺a,bは奇数と偶数の組で、斜辺cは奇数です。
　なので、ループ変数iは3から１辺最大値sまで2飛びで値を増加させます。
・変数tは二等辺の長さを2で割って切り捨てることで、
擬正三角形の３辺がi,i,i-1の場合の半分直角三角形の直角を挟む１辺の長さです。
そして、t+1が擬正三角形の３辺がi,i,i+1の場合の半分直角三角形の直角を挟む１辺の長さになります。

・二乗値辞書dから取得した、上記iとtの二乗値の差は、三平方の定理より
　擬正三角形の垂線hの長さの二乗値と等しく、この値が辞書ddのキーにあれば整数直角三角形なので、擬正三角形の周囲の長さを変数cに累積していきます。

・・・と、ここまできて、１分ルールが守れず長考に入ってしまいました・・・。

def f2(n):
	r, s, t, x, y = 0, 0, 0, 1, 0
	while s<=n or t<=n:
		x, y = 2*x+3*y, x+2*y
		
		s = x*x*4
		if s<=n: r += s
		
		t = x+3*y
		t = 2*t*t
		if t<=n: r += t
		
	return r

n = 10**9
print f2(n)

改めて、問７５より、３辺が整数の直角三角形では、整数i,j（0<j<i）を使って、
> 斜辺がi2+j2、直角を挟む２辺はi2-j2, 2ij

よって、上記の半分直角三角形に当てはめると、
a = i2+j2　：(A)
(b, h) = (i2-j2, 2ij) or (2ij, i2-j2)　：(B)
擬正三角形の周長L = 2(a+b)　：(C)
擬正三角形の定義から、a = 2b±1　：(D)

・ケース１：b = i2-j2の場合
(D)に(A)と上記を代入して
i2 + j2 = 2(i2 - j2)±1
i2 + j2 = 2i2 - 2j2 ±1
i2 - 3j2 = ±1　：(E)
また、(C)に(A)と上記を代入して
L = 2((i2+j2)+(i2-j2))
  = 4i2　：(F)

・ケース２：次に上記(B)で、b = 2ijの場合
(D)に(A)と上記を代入して
i2 + j2 = 2・2ij ±1
i2 + j2 = 4ij ±1
(i - 2j)2 - 3j2 = ±1　：(G)
また、(C)に(A)と上記を代入して
L = 2(i2+j2 +2ij)
  = 2(i+j)2　：(H)
ここで、t=i-2jと置くと、(G)、(H)は、
t2 - 3j2 = ±1　：(G')
L = 2((t+2j) +j)2
  = 2(t+3j)2　：(H')
　
・上記の２つの場合をまとめると、
x2 - 3y2 = ±1　：(J)
L = 4x2, 2(x+3y)2　：(K)

・さてここで、（J)の「±1」の「-1」の場合を考えてみます。
x2 - 3y2 = -1
x2 = 3y2 -1
xの二乗は、３の倍数に１足りないので、
xの二乗を３で割ると２余ることになります。
x = 3m, 3m+1, 3m+2　として確認していきます。
x = 3mの場合、x2 = (3m)2 = 6m、3で割ると余り0
x = 3m+1の場合、x2 = (3m+1)2 = 9m2+6m+1、3で割ると余り1
x = 3m+2の場合、x2 = (3m+2)2 = 9m2+12m+4、3で割ると余り1
上記より、二乗して３で割ったときに２余る整数は存在せず、式(J)は以下に改めます。
x2 - 3y2 = 1　：(J')

式(j')はペル方程式と呼ばれる式で、以下の公式があります。
x2 - Dy2 = 1　（Dは平方数以外。平方数の場合、整数解が存在しない）
の解は、最小整数解を(x0, y0)、それ以降、(x1, y1),(x2, y2), ...として、
┌xn+1 = x0xn + y0Dyn
└yn+1 = y0xn + x0yn

当問題に当てはめると、D=3、また、、
問題にある5-5-6の擬正三角形が最小整数解で、このときx0=2, y0=1 なので、
┌xn+1 = 2*xn + 3*yn
└yn+1 = 1*xn + 2*yn

上記式をループ内に実装し、上記(K)の擬正三角形の周の長さL = 4x2, 2(x+3y)2　を計算し累積します。




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

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

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

問題を再掲します。

問９３「算術式」
集合 {1, 2, 3, 4} の各数字をちょうど一度用い, また四則演算 (+, －, *, /) と括弧を使うことにより, 異なる正の整数を作ることができる.

例えば,
8 = (4 * (1 + 3)) / 2
14 = 4 * (3 + 1 / 2)
19 = 4 * (2 + 3) － 1
36 = 3 * 4 * (2 + 1)

12 + 34 のように数字をつなげることは許されないことに注意しよう.

集合 {1, 2, 3, 4} を使うと, 36 を最大とする 31 個の異なる数が得られる.
最初の表現できない数に会うまで, 1 から 28 の各数を得ることができる.

最長の連続した正の整数 1 から n の集合を得ることができる, 4 つの異なる数字 a < b < c < d を見つけよ. 答えを文字列 abcd として与えよ.






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






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

def f(n):
	import itertools
	d = {}
	for wi in xrange(1, n+1):
		for xj in itertools.combinations(xrange(1, 10), wi):
			zk = int("".join([str(i) for i in xj]))
			d[zk] = []
			for vm in range(wi//2):
				for pk in itertools.combinations(xj, vm+1):
					a = int("".join([str(i) for i in pk]))
					b = int("".join([str(i) for i in xj if i not in pk]))
					if a>b: continue
					for i in d[a]:
						for j in d[b]:
							d[zk] += [i+j, i-j, j-i, i*j, 1.0*i/j, 1.0*j/i]
							
			if d[zk]: d[zk] = set([s for s in d[zk] if s])
			else: d[zk] = [zk]
			
	vol = 0
	for i in d.keys():
		if i<10**(n-1): continue
		j = 1
		while j in d[i]: j += 1
		if vol<j-1: ans, vol = i, j-1
		
	return ans, vol
	
n = 4
ans, vol = f(n)
print "ans:",ans
print "vol:",vol

1桁少ない桁数までの計算結果を再利用する方法です。
数式に使用する数字の個数を順に増やしていきます。
総当りよりも速く、4桁以上の桁数でも動作します。
初期設定として、1桁の場合にその数自身だけの配列（リスト）を設定します。

桁数ループwで1桁ずつ増やしつつ、
1から9までの数字から重複無しで桁数分w個だけ数字を選びます。
この数字の組xのすべての並び替えパターンを1つずつ試していきます。
数字の組xの要素で構成された数値zが計算対象です。

この数値を表記上の前半と後半の2つに分けます。
前半の数aが後半の数bよりも大きいとすることで重複計算を防ぎます。
このa,bはzよりも必ず桁数が小さいので、すでに計算済みです。
そしてzの計算として、a,bの計算結果同士を総当りで四則演算します。
計算前に結果をためる配列を、数値zをキーにした連想配列（辞書）dの要素として準備しておきます。

足し算と掛け算は演算子の前後を入れ替えても同じですが、
引き算と割り算は演算子の前後を入れ替えた値も計算し配列にためていきます。

すべてのzについて計算が終わったら最長連続の判定をします。
n桁の数字iをキーにして、計算結果辞書dの要素、つまり数値iの計算結果配列d[i]について、
1から順にいくつまで連続して要素があるかを確認します。

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

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

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

問９３「算術式」
集合 {1, 2, 3, 4} の各数字をちょうど一度用い, また四則演算 (+, －, *, /) と括弧を使うことにより, 異なる正の整数を作ることができる.

例えば,
8 = (4 * (1 + 3)) / 2
14 = 4 * (3 + 1 / 2)
19 = 4 * (2 + 3) － 1
36 = 3 * 4 * (2 + 1)

12 + 34 のように数字をつなげることは許されないことに注意しよう.

集合 {1, 2, 3, 4} を使うと, 36 を最大とする 31 個の異なる数が得られる.
最初の表現できない数に会うまで, 1 から 28 の各数を得ることができる.

最長の連続した正の整数 1 から n の集合を得ることができる, 4 つの異なる数字 a < b < c < d を見つけよ. 答えを文字列 abcd として与えよ.






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






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

def f():
	import itertools
	r = 0
	#括弧の位置  [0]a-[1]b[2]-[3]c[4]-d[5]
	kakko = [("(((","",")","",")",")"), ("(","",")","(","",")"), ("","(","","",")","")]
	for t in itertools.combinations(xrange(1, 10), 4):
		D = {}
		for b in itertools.permutations(t):
			f = [str(i)+".0" for i in b]
			for m in itertools.product("+-*/", repeat=3):
				for k in kakko:
					s0= k[0]+f[0]+m[0]+k[1]+f[1]+k[2]+m[1]+k[3]+f[2]+k[4]+m[2]+f[3]+k[5]
					try: s1 = eval(s0)
					except ZeroDivisionError: continue
					if int(s1)==s1: D[int(s1)] = s0

		for i in xrange(1, max(D.keys())+1):
			if i not in D.keys(): break

		if r<i-1: p, r, w = t, i-1, [D[j] for j in xrange(1,i)]

	p = "".join([str(i) for i in p])
	w = [s.replace(".0","")+"="+str(int(eval(s))) for s in w]
	return p,r,w

ans, vol, pat = f()
print "ans:",ans
print "vol:",vol
print "pat:",pat

総当りです。１分ルールには間に合いました。

式の文字列として作り出し、eval関数で計算式として計算します。

・9個の数字から異なる４つを選ぶパターンは、重複なし組合せで、9C4 = 126通り。

・4つの異なる数字の並び順パターンは、重複なし順列で、4P4 = 4! = 24通り。

・四則演算を3回行うパターン数は、重複あり順列で、4の3乗で64通り。

・数字が4個のとき、カッコのパターン数は並1べ順替えを考慮すると以下の3通り。

例えば以下のとおり。

(((9 - 5)- 2)- 1) = 9-5-2-1 = 1

  (9 - 5)-(2 - 1) = 9-5-2+1 = 3

   9 -(5 - 2)- 1  = 9-5+2-1 = 5

126x24x64x3 = 580,608通りを総当りで試しました。

文字型で式を作成した後、eval関数で計算式そのものに変換して実行します。

割り算記号がどこに来ても切り捨てられないように、全ての数字の後ろに「.0」を追加しておきます。

割り算でゼロ除算になる場合は、try-exceptでつかまえて次のケースに進みます。

計算結果の小数有無判定はint関数結果と同じ値になるかで判定します。

計算結果が正の整数値の場合、辞書（連想配列）Dに式そのものをためていきます。

連続判定では、リストDのキーに1からの連続値があるところi番目の1つ前、i-1番目まで連続です。

最長判定では、変数rに連続パターン件数をとっておき、最長かどうかを判定します。

求められている数字の組のパターンの他、連続個数とその式も表示します。

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

ans: 1258

vol: 51

pat: ['(8-5)/(2+1)=1', '(8-5)-(2-1)=2', '(8-5)/(2-1)=3', '8-(5-2)-1=4','(((8-5)*2)-1)=5', '(8-5)*(2/1)=6', '(((8-5)*2)+1)=7', '(((8-5)+1)*2)=8', '(8-5)*(2+1)=9', '8+(5-2)-1=10', '8+(5-2)/1=11', '(8+5)-(2-1)=12', '(8+5)/(2-1)=13', '8+(5+2)-1=14', '8+(5+2)/1=15', '8+(5+2)+1=16', '8+(5*2)-1=17', '8+(5*2)/1=18', '8*(5/2)-1=19', '8*(5/2)/1=20', '8*(5/2)+1=21', '(8*2)+(5+1)=22', '8*(5-2)-1=23', '8*(5-2)/1=24', '8*(5-2)+1=25', '(8+5)*(2/1)=26', '(((8+5)*2)+1)=27', '(((8+5)+1)*2)=28', '(((8-2)*5)-1)=29', '8*(5-1)-2=30', '(((8-2)*5)+1)=31', '(((5-2)+1)*8)=32', '(((8-1)*5)-2)=33', '8*(5-1)+2=34', '(((8-2)+1)*5)=35', '(8-2)*(5+1)=36', '(((8*5)-2)-1)=37', '(8*5)-(2/1)=38', '(8*5)-(2-1)=39', '(8*5)/(2-1)=40', '(8*5)+(2-1)=41', '(8*5)+(2/1)=42', '(8*5)+(2+1)=43', '(((1/2)+5)*8)=44', '(((8+2)-1)*5)=45', '8*(5+1)-2=46', '(((8+1)*5)+2)=47', '(((5+2)-1)*8)=48', '(((8+2)*5)-1)=49', '8*(5+1)+2=50',  '(((8+2)*5)+1)=51']

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

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

問９２「桁の2乗による連鎖」
各桁の2乗を足し合わせて新たな数を作ることを, 同じ数が現れるまで繰り返す.

例えば

44 → 32 → 13 → 10 → 1 → 1
85 → 89 → 145 → 42 → 20 → 4 → 16 → 37 → 58 → 89

のような列である. どちらも1か89で無限ループに陥っている. 
驚くことに, どの数から始めても最終的に1か89に到達する.

では, 10,000,000より小さい数で89に到達する数はいくつあるか.





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






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

def f(n):
	P = [i*i for i in xrange(10)]
	s = len(str(n))
	m = 9*9*s
	
	L = [0]+[-1]*(m-1)
	L[1], L[89] = 0, 1
	for i in xrange(1, m):
		j = i
		while L[j]<0:
			j, r = 0, j
			while r:
				j += P[r%10]
				r //= 10
			L[i] = L[j]
			
	M = [0]*m
	for i in xrange(10): M[P[i]] = 1
	for u in xrange(2, s):
		W = [0]*m
		for i in xrange(m):
			for j in xrange(10):
				if P[j]<=i: W[i] += M[i-P[j]]
		M = W[:]
		
	return sum([i*j for i, j in zip(L, M)])
	
n = 10000000
print f(n)

1から10,000,000まで、桁の2乗和の連鎖を総当りで計算すると、とても1分ルールに合いません。

試行錯誤したところ、桁2乗和は１回計算するだけで値がずいぶんちいさくなることがわかりました。
桁数sのときの桁2乗和最大値mは9がs個並んだ数で、m = 9*9*s です。
7桁ならば、m = 9*9*7 = 567 です。

そこで、
・m以下の数は桁2乗和の連鎖を最後まで計算して89ループになるか判定しておく。
・桁2乗和は１回だけ計算しm以下にしたところで数える。
ことにしました。

・89ループ判定結果リストL
1からmまでの桁2乗和の連鎖を最後まで行い、89ループになる場合は1、その他は0を設定します。
-1で初期設定しているので、while L[j]<0として未設定が続く限り連鎖をたどります。
rのwhileループで、rの10で割った余りとして、下から1桁ずつ取り出して2乗和を累積します。
2乗の計算はいちいち行わず、2乗リストP[0,1,4,9,16...]を参照します。

・桁2乗和ごとの件数リストM
初期値として1桁数値の桁2乗和の位置(1,4,9...)に1を設定します。
2桁以上の数値については、１つ少ない桁数までの累積件数の、最上位の数値の2乗分だけ小さい位置にある値を、採用し累積していきます。
このようにしてだんだん桁数を増やしていきます。

こうして、
・89ループ判定結果リストL（89ループに達する数の位置に1が立ち、他は0が設定）
・桁2乗和を１回だけ計算した値ごとの件数リストM
の同じ位置の要素同士の積を求めることで89ループの部分を取り出します。
この総和が求める件数です。


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