2015年1月3日土曜日

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

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

問102「三角形の包含」
3つの異なる点が -1000 ≦ x, y ≦ 1000 かつ三角形となるように, デカルト平面上にランダムに与えられる.

以下の2つの三角形を考える.

A(-340,495), B(-153,-910), C(835,-947) 
X(-175,41), Y(-421,-714), Z(574,-645)
三角形ABCが原点を内部に含み, XYZは原点を内部に含まないことが確かめられる.

27Kのテキストファイルtriangles.txt(右クリックしリンク先を保存して欲しい) にランダムな1000個の三角形が適当なフォーマットのもと含まれている. 内部に原点を含む三角形の数を答えよ.

注: ファイル中の最初の二つは三角形ABC, XYZである.






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






私の回答例は以下の通りです。
def f(fn):
	c = 0
	p = open(fn)
	for s in p.readlines():
		t, L, M = s.split(","), [], []
		for i in xrange(0, len(t), 2): L.append((int(t[i]), int(t[i+1]))) 
		L.append(L[0])
		for i in xrange(len(t)//2):
			a, b = L[i], (L[i+1][0]-L[i][0], L[i+1][1]-L[i][1])
			M.append(a[0]*b[1] - a[1]*b[0])
		if abs(sum(M))==sum([abs(u) for u in M]): c += 1
	p.close()
	p = None
	return c

s = "p102_triangles.txt"
t = os.path.join(os.path.dirname(__file__), s)
print f(t)



原点が三角形の各辺の左右同じ側にあるかという位置関係で判定します。
ベクトルの外積を利用します。

対象となる点(当問題では原点)が三角形の内部にあれば、
その点は三角形を一周して常に同じ側にあることになります。
つまり、対象となる点(当問題では原点)から各頂点へのベクトルと、
その頂点から次の頂点へのベクトルの外積の符号がすべて同じになります。
なお、平面上の話なので、外積はベクトルにならずスカラー量です。

aベクトル:原点から三角形のi番目の頂点へのベクトル
bベクトル:i番目の点からi+1番目の点へのベクトル
として、
aベクトルとbベクトルの外積の符号が、
正ならばaベクトルが右側、つまり原点はbベクトルの左側で、
負ならばaベクトルが左側、つまり原点はbベクトルの右側です。
3つの外積がすべて正ならば原点を右回りに一周し、
3つの外積がすべて負ならば原点を左回りに一周することになります。

外積の公式は以下のとおりです。
aベクトル=(Ax, Ay)、bベクトル=(Bx, By)とすると、
a×b = Ax*By - Ay*Bx

以上を実装しました。

1.関数f(fn)
・入力ファイルのフルパスを受け取り、問題にある内部に原点を含む三角形の数を返します。
・ループ変数sは、入力ファイルの1行分で、6つの数字のcsv文字列です。
・変数tは、6つの数字(まだ文字型)を持つリストです。
・一周して計算する都合上、最初の座標を最後にもう一度追加し
座標を4つ分にしてから外積計算を3回行います。
リストLは、平面座標値(タプル型)を最初3つ持ち、最初の1つを末尾に追加して
計4つの座標値を持ちます。
・変数a,bはそれぞれ上記の説明にあるaベクトルとbベクトルです。
・リストMには外積値をためていきます。
・外積の符号一致判定は、外積値の和の絶対値=外積値の絶対値の和かどうかで行います。



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

2015年1月1日木曜日

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

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

問101「最適多項式」
数列のk個の項を与えられたときに, 次の項を確実に求めることは不可能である. 
その数列に合うような多項式が無限個存在するからである.

例として, 立方数の数列を考えよう. これは生成関数 un = n3 で定義され, 1, 8, 27, 64, 125, 216, ...となる.

この数列の最初の2項のみが与えられているとしよう. 
"Simple is best"の法則にのっとり, 線形の関係があると仮定し, 3つ目の項が15であると予想する (差分が7). もし最初の3項のみが与えられていたとしても, 同じ原則により, 二次の関係があると仮定して次の項を予測する.

数列の最初のk項を生成できる最適な多項式のn項を OP(k, n) で表すことにする. 
明らかに, n ≦ k について OP(k, n) は正しい. 
最初の異なる項 (First Incorrect Term, FIT) は OP(k, k+1) であろう. 
これを bad OP (BOP) と呼ぶことにする.

原則より, 最初の項しか与えられていない場合には, 定数項とするのが理に適っているだろう;即ち, n ≧ 2, OP(1, n) = u1.

従って, 立方数の数列について以下のOPを得る.

OP(1, n) = 1 1, 1, 1, 1, ...
OP(2, n) = 7n−6 1, 8, 15, ...
OP(3, n) = 6n2−11n+6      1, 8, 27, 58, ...
OP(4, n) = n3 1, 8, 27, 64, 125, ...

明らかに, k ≧ 4 のときにはBOPは存在しない.
BOPのFIT (上の例ではで示されている) の和は, 1 + 15 + 58 = 74 である.

以下の10次多項式からなる生成関数を考える:
un = 1 - n + n2 - n3 + n4 - n5 + n6 - n7 + n8 - n9 + n10
BOPのFITの総和を求めよ.






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






私の回答例は以下の通りです。
def g(Y):
	x, y = len(Y), 0
	Xi = [i for i in xrange(x)]
	for i in xrange(x):
		a, b = Y[i], 1
		for j in xrange(x):
			if i==j: continue
			a *= (x - Xi[j])
			b *= (Xi[i] - Xi[j])
		y += a//b
	return y

def f(s):
	L, FIT, BOP, t, n = [], 0, 0, -1, 0
	while FIT!=t:
		BOP += FIT
		n += 1
		FIT, t = g(L), eval(s)
		L.append(t)
	return BOP

u = "1 - n + n**2 - n**3 + n**4 - n**5 + n**6 - n**7 + n**8 - n**9 + n**10"
BOP = f(u)
print "BOP:", BOP



この問題は解法が思いつかず、ネット検索したところ、「ラグランジュの補間多項式」の問題と判明。
ラグランジュ補間(wiki) を参考にロジックを組みました。

1.関数g(Y)
・ラグランジュの補間多項式を組み込んだ関数です。
 数値リストYから要素数より1つ少ない次元の補間多項式による、次の要素を返します。
・数値xはリストYの要素数で、求める要素の要素番号です。
 リストXiは、リストYの要素番号のリストです。
・ループ変数iは補間式内の項の番号です。
 ループ変数jは補間式の各項内の掛け算要素の番号です。
・変数a, bは、その掛け算要素の分子と分母です。
 項1つ分の分子と分母の計算が終わったら、その商を変数yに累積していきます。

2.関数f(s)
・文字型で数式を受取り、問題文にあるBOPのFITの総和を返します。
 数式の変数は1種類、nだけとします。
・リストLは、受け取り数式sによる実値をためていきます。
・whileループ内は、現在までの実値リストに基づく推定値FITと
 数式による実値tとが不一致の間、処理を実行します。
 また、初回n=1のときは一致することがあるので、実値リストLがカラのときは処理を行うようにします。
・受け取り数式sの変数nに1から順に値を代入していき、関数g()を使って推定値FITを求めます。
 文字列である数式をeval関数で処理する数式として認識させ実行します。



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

2014年12月27日土曜日

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

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

問100「組み合わせの確率」
箱の中に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


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

円盤総数を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が求める値です。

問94 より、ペル方程式の公式は以下の通り。
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は両方とも奇数。
これを実装します。

1.関数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の積を2で割った余りが1(論理判定ではTrue)になるかで判定します。
 そして、xとyが両方とも奇数のときに、式(b)から割り戻した値を
 円盤総数iと青円盤数jに退避しておきます。


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

2014年12月20日土曜日

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

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

問99「最大のべき乗」
指数の形で表される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


1.関数f(sIn)
・aのb乗をpとすると、
 p = ab
 両辺の10を底とする対数をとって、
 log(p) = b x log(a)
 なお、上記式の両辺の対数の底を10にすれば、
 この両辺の値を切り捨て+1した値が桁数になります。
 この式を実装して最大値を判定します。
 
・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

2014年12月10日水曜日

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

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

問98「アナグラム平方数」
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


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

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

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

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


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

2014年11月24日月曜日

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

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

問97「大きな非メルセンヌ素数」
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)

#(a×(nのp乗)+b)の末尾m桁を返す
a, n, p ,b ,m =28433, 2, 7830457, 1,10
s = f(a, n, p ,b ,m)
print s

#(a×(2のp乗)+b)の末尾m桁を返す
a, p ,b ,m =28433, 7830457, 1,10
s = f2(a, p ,b ,m)
print s


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

・for s in str(p)[::-1]:
 数値pをstr関数で文字型にして、文字のスライス機能の[開始:終了:刻み]での刻みを-1にすることで右から1文字ずつ、つまり1の位の数字から順に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です。

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

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

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

2014年11月4日火曜日

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

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

問96「数独」
"数独"(日本語で数字を配置するという意味)とは人気があるパズルの名前である.
起源は不明だが, その評判はラテン語で"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マスの位置を示す座標タプルで、
値は候補値の配列(リスト型)、または確定値(整数型)です。

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

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

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












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











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






















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





















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















8.その他
部品となる関数を以下の通りです。
・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