Pythonではtimeitモジュールを使って簡単に実行時間を計測できる。以下に示すコードはPythonでは一般的なコードだと思う。たらい回し関数(竹内関数)と階乗を求める関数だ。これを実行した時の時間計測、さらにコードに含まれている関数の実行時間を簡単に測るにはどうすればよいか。
ここでtimeitモジュールが利用できる。まず、Pythonシェル環境では以下の通り。
main関数の測定は次のように行う。stmtで時間計測したい関数を指定して、setupで最初の1回だけ実行する文を指定する。numberは実行する回数となっている。正確を期するなら複数回を指定するべきだろう。
>>> import timeit
>>> timeit.timeit(stmt='test_timeit.main()', setup='import test_timeit', number=1)
12
3628800
2.663659573618423
たらい回し関数(tarai)のみの測定。
>>> timeit.timeit(stmt='test_timeit.tarai(12, 6, 0)', setup='import test_timeit', number=1)
2.6390463109480637
階乗関数(factorial)のみの測定。100万回実行。
>>> timeit.timeit(stmt='test_timeit.factorial(10)', setup='import test_timeit', number=1000000)
2.3606330638673825
さらに、コマンドラインでの計測も可能だ。
main関数の測定。コマンドオプションの-sは最初に1回だけ実行する文、-nは実行する文の回数、-rはタイマのリピート数になっている。
$ python -m timeit -n 1 -r 1 -s "import test_timeit" "test_timeit.main()"
12
3628800
1 loops, best of 1: 2.63 sec per loop
たらい回し関数(tarai)のみの測定。
$ python -m timeit -n 1 -r 1 -s "import test_timeit" "test_timeit.tarai(12, 6, 0)"
1 loops, best of 1: 2.63 sec per loop
階乗関数(factorial)のみの測定。100万回実行。
$ python -m timeit -n 1000000 -r 1 -s "import test_timeit" "test_timeit.factorial(10)"
1000000 loops, best of 1: 2.33 usec per loop
ところで、上に示した階乗関数については「車輪の再発明」なので、通常はmath.factorialを使った方が良いだろう。
2013年3月10日日曜日
Python: timeitモジュールを使ってお手軽に実行時間計測
2013年1月1日火曜日
2013年と素因数分解
あけましておめでとうございます。
新年を迎えて、ふと、今年の西暦である2013を素因数分解したらどんな数に分解されるのか気になったので考えてみた。桁をすべて足し合わせると 2+0+1+3=6 と3の倍数になり、3で割り切れることは明白だ。3で割ると671となる。671に対し桁を一つ飛ばしにして足し合わせた数の差が (6+1)-7=0 で11の倍数(0を含む)なので11で割り切れることも簡単にわかる。671を11で割ると61の素数となる。つまり、2013の素因数分解は 3×11×61 となる。
大きな素数同士による合成数の素因数分解は非常に困難であることが知られている。この性質を利用して多くの暗号アルゴリズムで大きな素数による合成数が利用されている。現在のところ、素因数分解アルゴリズムとしては、楕円曲線法(ECM; Elliptic Curve Method)や複数多項式二次ふるい法(MPQS; Multiple Polynomial Quadratic Sieve)がよく使われているらしい。で、手軽にこれらを利用できるライブラリがないかと軽く調べたところNZMATH(ニジマス)という数論のためのPythonモジュールがあることを知ったので、早速インストールして使ってみた。
ところで、大きな素数といえばメルセンヌ素数が有名だが、マラン・メルセンヌは、2n-1 に対してnが257以下のとき、n = 2, 3, 5, 7, 13, 17, 19, 31, 67, 127, 257のときにのみ素数となると主張した。しかし、一部に間違いがあり、n = 61, 89, 107のときも素数であり、また、n = 67, 257は素数ではなく合成数であった。
そこで、誤って素数としていたn = 67のときの 147573952589676412927 について楕円曲線法(ECM)で解かせると、すぐに 193707721×761838257287 と結果が出力された。
>>> import nzmath.factor.methods as methods
>>> methods.factor(147573952589676412927, method="ecm")
[(193707721L, 1), (761838257287L, 1)]
更に、メルセンヌが見落としていた9番目と10番目のメルセンヌ素数である 2305843009213693951 (261-1) と 618970019642690137449562111 (289-1) を掛け合わせた46桁の合成数 1427247692705959880439315947500961989719490561 を複数多項式二次ふるい法(MPQS)で解いてみたところ、30分ほどで正しく2つの素数に分解できた。
>>> methods.factor(1427247692705959880439315947500961989719490561, method="mpqs")
[(2305843009213693951L, 1), (618970019642690137449562111L, 1)]
ソースコードも簡単に確認できるし、NZMATHはなかなか楽しい。
それでは、本年もよろしくお願い申し上げます。
2012年9月15日土曜日
「フカシギの数え方」の問題を解いてみた
先日、「『フカシギの数え方』 おねえさんといっしょ！ みんなで数えてみよう！」という動画を見た。格子状のマスの左上から右下までの経路が何通りあるのかを調べて、格子が多くなればなるほど組み合わせの数が爆発的に増えることを教えてくれる動画だ。これは自己回避歩行(Self-avoiding walk)と呼ばれている問題らしい。
これだけ聞いてもそれほどインパクトはないのだが、動画に出てくるおねえさんの経路を調べあげる執念がもの凄く、ネット上でも結構な話題になっている。執念と言うよりも狂気に近い。しかし、話題になった割には動画内で言及されている高速なアルゴリズムを実装したという話を聞かなかったので、自分で確かめることにした。
動画のおねえさんは深さ優先探索によるプログラムを使っていると思われるが、それだとスパコンを使っても10×10マスの格子を解くのに25万年も掛かってしまう。そこで、高速化のためにゼロサプレス型二分決定グラフ(ZDD; Zero-Suppressed Binary Decision Diagram)と呼ばれるアルゴリズムを利用することにした。このアルゴリズムを開発したのは北大の湊先生で、ZDDによりすべての経路を見つけ出すアルゴリズムとしてクヌース先生のSIMPATHを使った。ZDDについてはクヌース先生も強い関心を持っていて、 The Art of Computer Programming Volume 4, Fascicle 1
(TAOCP 4-1)ではBDD/ZDDの詳細な解説を読むことができる。演習問題の解説だけで書籍の半分を使っていることからしても気合の入れようがわかるだろう。実際に自分のノートPCでZDDアルゴリズムを使ったコードを走らせたら、ほんの10秒程度で10×10マスの問題を解いてしまった。おねえさんがスパコンで25万年かかった問題をノートPCでたった10秒である。約8千億倍の高速化だ。これだけ劇的に変わるとやっぱり楽しい。そして、アルゴリズムの重要性を再認識させられた。
さて、以下におねえさんが利用したであろう自作の深さ優先探索(DFS)プログラムと高速な解法であるZDDアルゴリズムの両方を載せておく。ZDDについては4つのプログラムを1つのPythonスクリプトで統合している。これは、クヌース先生のSIMPATHおよびSIMPATH-REDUCEを利用しており、また、SGB (Stanford GraphBase)ライブラリでグラフを作成し、経路の数え上げにGMPを使ったC++で自作コードを組んだためである。最初はSIMPATHを参考に実装して一つのソースコードに統一しようかと思ったが、解説もちゃんと書かれているクヌース先生のコードをそのまま利用することにした。これらのプログラムの簡単な解説と使い方を載せておく。
まずはDFSによる実装を以下に示す。Node構造体から経路を辿って数え上げるだけの単純なプログラムである。
実行方法は以下の通り。引数の7は7×7のノードを表しており、マスで表現すれば6×6マスとなる。そして計算結果として、575,780,564通りの経路が探索されたことを示している。因みに7×7ノード(6×6マス)で約3分ほどの時間が掛かっており、それ以上の大きさだと現実的な時間で解くことが難しくなってくる。
% ./count_lattice_path_dfs 7 Count all paths from (0, 0) to (6, 6) in the 7x7 lattice graph: Count: 575780564
次に、ZDDを利用したプログラムのビルド方法を示す。
ビルドをする前にSGBおよびGMPを準備しておく。また、SIMPATHおよびSIMPATH-REDUCEについてはCWEBを利用しているので、それも入れておく。
最初に、SGBライブラリを用いて格子状のグラフを作成する。
% gcc gen_lattice.c -O3 -lgb -o gen_lattice
ソースコード内でグラフを下記のように作成している。board関数はチェスのような格子状のグラフを作成することができる。他にもグラフを作成する便利な関数が多数あるので、興味のある方はSGBライブラリのgb_basic.wあたりを参照して欲しい。
Graph* g = board(N, M, 0L, 0L, 1L, 0L, 0L);
次に、SIMPATHおよびSIMPATH-REDUCEを以下のようにビルドする。
% wget http://www-cs-faculty.stanford.edu/~uno/programs/simpath.w % ctangle simpath.w % gcc simpath.c -O3 -lgb -o simpath
% wget http://www-cs-faculty.stanford.edu/~uno/programs/simpath-reduce.w % ctangle simpath-reduce.w % gcc simpath-reduce.c -O3 -lgb -o simpath-reduce
因みに、ctangleでC言語への変換、cweaveでTeX形式への変換を行うことができるので、以下のようにTeX形式のドキュメントも一緒に作っておくことをお勧めする。また、これらのアルゴリズムを理解するために前述のTAOCP Vol. 4を読んでおくことが望ましい。
% cweave simpath.w % tex simpath.tex % dvipdf simpath.dvi # PDFファイルが必要なら. % cweave simpath-reduce.w % tex simpath-reduce.tex % dvipdf simpath-reduce.dvi # PDFファイルが必要なら.
最後に、ZDDから経路を数え上げるプログラムをコンパイルする。
% g++ count_path_mp.cpp -O3 -lgmp -o count_path_mp
経路を一つ一つ数えていては高速に数え上げるという点で意味がなくなってしまうので、ちゃんとメモ化再帰で数える。ZDDは高度に圧縮されているのでメモ化再帰でほとんど瞬時に数え上げることができる。
そして、上記のプログラムを統括するPythonスクリプトが以下となる。グルー言語としてもPythonは優秀だ。
実行方法は以下の通り。
% ./count_lattice_path.py 8 Count all paths from (0, 0) to (7, 7) in the 8x8 lattice graph: Count: 789360053252
DFSのプログラム同様、引数の8は8×8ノード(7×7マス)を表している。そして計算結果として、789,360,053,252通りの経路が探索されたことを示している。また、中間ファイルとして以下が出力される。
lattice_08.gb # SGBによるグラフデータ. lattice_08.out # SIMPATHの出力: not-necessarily-reduced BDD lattice_08.out.r # SIMPATH-REDUCEの出力: ZDD lattice_08.out.rc # SIMPATH-REDUCEの出力を再番付. lattice_08.log # 一連の実行ファイルのログ.
因みに、DFSで7×7ノード(6×6マス)を解くと3分ほどかかったが、ZDDでは0.1秒未満で完了した。また、ZDDで14×14ノード(13×13マス)は数分で終わることを確認したが、それ以上は途中で搭載メモリ以上に必要なメモリが大きくなるので実行していない。
2012年5月14日月曜日
テトレーション・フラクタル
以前の記事では、∞i が 0.43828+0.36059i に収束することを見つけたのだが、今回はそれに関連したフラクタルについて紹介したい。
テトレーション na の a には複素数を指定することができる。このとき、a = x + yi として、それを複素平面に置く。ここで正の整数nを大きくしていき、発散と判定されたnに対応した色で平面を色分けする。発散しない場合、予め決めておいたnの上限値を使う。これで得られる図をテトレーション・フラクタル(tetration fractal)と呼ぶ。
n(x + yi)=a + biとしたとき、n+1(x + yi)=a' + b'i は以下のように計算できる。
Pythonでの実装はActiveSate CodeにTetration Fractalとして載っている。今回はこれを少し修正したコード、C++で書いたコード、更にTBBで並列化したコードを用意して、それらの実行速度のベンチマークを取ってみた。
| Python | C++ | C++ with TBB |
| 258.732秒 | 16.999秒 | 0.976秒 |
ここで、複素平面の領域は(-1.5, 0)-(-0.75, 0.75)であり、最大繰り返し数は256としている。画像の大きさは1,024×1,024とした。実行環境はIntel Xeon X5650の2CPUで、12コア・24スレッドとなっている。
以下に、ソースコードを示す。詳細についてはコードを読んでいただきたい。
Pythonによる実装:
% ./tetration.py tetration.png -1.5 0.0 0.75 0.75 1024 1024
C++による実装:
% g++ tetration.cpp -o tetration -O3
% ./tetration tetration.out -1.5 0.0 0.75 0.75 1024 1024
% ./image_tetration.py tetration.out tetration.png
TBBを利用したC++による並列化実装:
% g++ tetration_tbb.cpp -o tetration_tbb -ltbb -O3
% ./tetration_tbb tetration.out -1.5 0.0 0.75 0.75 1024 1024
% ./image_tetration.py tetration.out tetration.png
C++の実行ファイルで生成された出力ファイルを画像ファイルに変換するスクリプトを以下に示す。
2012年3月31日土曜日
IPythonの埋め込みプロットが素晴らしい
先日、Tokyo.SciPy #3に参加して、SciPyの生みの親であるTravis Olipant氏のセッションを聞いた。その時に最近のIPythonでは埋め込みプロットができることを知ったので早速入れてみたところ、これがとてもクールでカッコよかったのでここで紹介したい。
埋め込みプロットはIPythonのバージョン0.11からできるようになっている。ただ、自分が使っているLinux環境では普通に持ってくるとバージョン0.10が入ってしまうので使えない。そこで、IPythonの公式サイトから最新版のバージョン0.12を持ってきて入れてみた。
ただ、埋め込みプロットができるのはターミナル版ではなくQt版のシェルなので、それに関連するライブラリなども一緒に入れなくてはならない。IPythonの起動時に何々のライブラリが足りないとかいろいろと文句を言われるが、足りないものをyumなりapt-getなりソースコードをコンパイルなりして順次入れていけば使えるようになると思う。
あと、IPythonの設定ファイルの構成が0.10から大きく変わった。以前は .ipython/ 内の pythonrc を利用していたのだが、それがなくなって、代わりに .ipython/profile_default/ 内の ipython_config.py と ipython_qtconsole_config.py が使われるようになった。デフォルトの設定ファイルは以下のコマンドで作ることができる。
$ ipython profile create
そして、埋め込みプロットを使うためのQt版IPythonの起動コマンドは以下の通り。
$ ipython qtconsole --pylab=inline
これでQt版シェルが立ち上がって埋め込みプロットを表示できるようになる。
せっかくなので実際に埋め込みプロットを使ってみる。ここでは自分のお気に入りのマンデルブロ集合をプロットしてみた。因みに、シェル上での複数行のソースコードの編集も改良されたので入力がかなり楽になった。冒頭にスクリーンショットを貼っておく。
rangex = arange(-2.0, 1.0, 0.01)
rangey = arange(-1.0, 1.0, 0.01)
(X, Y), Z = meshgrid(rangex, rangey), []
for y in rangey:
Z_ = []
for x in rangex:
c, z = complex(x, y), 0.0
for i in range(100):
if abs(z) > 2.0: break
z = z * z + c
Z_.append(i)
Z.append(Z_)
contourf(X, Y, Z, arange(100.0))
2012年2月19日日曜日
Android端末をサーバにしてHTML5を使ったお絵かきBBSを作成する
最近、HTML5に触れる機会があり、その良さが何となく伝わってきたので、何かしら簡単なコードを書いてみたくなった。そこで、ちっちゃいけどリッチ、というギャップを楽しむためにAndroid端末を使うことにした。具体的には、Android端末をSL4Aでウェブサーバにして、HTML5をインターフェイスにしたお絵かきBBSをPythonで書いてみた。
BBSでは絵とテキストを扱うことができて、それらはAndroid端末上でSQLiteのデータベースで管理される。利用者の利便を考えて、名前などはクッキーで保存する。3G回線や無線LANなどで接続することを考慮して、IPアドレスも取得できるようにした。また、書き込み時にサーバのAndroid端末が振動して書き込みがあったことを知らせてくれる。因みに、NTTドコモの3G回線で使うためには、グローバルIPが振られるmopera Uなどのサービスが必要で、spモードでは利用できない。
実際に作ってみたプログラムのスクリーンショットを冒頭に入れてみた。必要最低限の機能しかないが、それでも文章と画像を保存できるちゃんとした掲示板だ。草の根BBSのころを考えると隔世の感がある。
以下に今回作成したソースコードを載せておく。こんな短いコードでもちゃんと機能するのが不思議な感じだ。
image_bbs.py
# -*- coding: utf-8 -*-
import sys,os,cgi,sqlite3,datetime
import socket,fcntl
from wsgiref.simple_server import make_server
import android
droid=android.Android()
LIMIT=10 # 最大表示記事数.
DB_FILE='/sdcard/image_bbs.sqlite'
P=8080
con=sqlite3.connect(DB_FILE)
cur=con.cursor()
cur.execute('CREATE TABLE IF NOT EXISTS bbs (id INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT, user TEXT, datetime TEXT, image TEXT, text TEXT)')
INSERT_DB='INSERT INTO bbs VALUES(NULL,?,?,?,?)'
def ipconfig():
s=socket.socket(socket.AF_INET,socket.SOCK_DGRAM)
s.connect(('gmail.com',80))
return s.getsockname()[0]
def post(user,image,text):
if image=='' or text=='': return
if user=='': user='匿名'
dt=datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')
cur.execute(INSERT_DB,(cgi.escape(user.decode('utf-8')),dt,cgi.escape(image.decode('utf-8')).replace('\n','<br />'),cgi.escape(text.decode('utf-8')).replace('\n','<br />')))
con.commit()
droid.vibrate()
#droid.notify(user,text)
def bbs(environ,start_response):
global IP
if environ['PATH_INFO']=='/':
user=''
if environ.has_key('HTTP_COOKIE'):
cookie={}
for d in environ['HTTP_COOKIE'].strip(';').split(';'):
d=d.split('=')
cookie[d[0]]=d[1]
if cookie.has_key('BBSUSER'):
user=cookie['BBSUSER']
if environ['REQUEST_METHOD']=='POST':
fs=cgi.FieldStorage(fp=environ['wsgi.input'],environ=environ,keep_blank_values=1)
user=fs.getfirst('user','').strip()
image=fs.getfirst('bbsImage','')
text=fs.getfirst('text','').strip()
post(user,image,text)
data=u"""<!DOCTYPE html>
<html><head>
<meta http-equiv="Set-Cookie" content="BBSUSER=%s">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>Image BBS by Android</title></head><body>
<p><a href="http://%s:%d">http://%s:%d</a></p>
<form action="/" method="post" id="form">
<span>名前:<input type="text" name="user" size="20" maxlength="30" value="%s" /></span>
<div><canvas id="bbsCanvas" width="400" height="200" style="border:1px solid;"></canvas></div>
<textarea name="text" cols="40" rows="5"></textarea><br />
<input type="hidden" name="bbsImage" id="bbsImage" />
<span><input type="button" onclick="clearCanvas();" value="画像消去" /></span>
<span><input type="button" value="更新" onclick="location.reload(true);" /></span>
<span><input type="button" value="送信" onclick="setImage();this.form.submit();" /></span>
</form>
<script type="text/javascript">
var mouseDownFlag = false;
window.onload = function() {
draw();
setImage();
};
function draw() {
var canvas = document.getElementById('bbsCanvas');
if (!canvas || !canvas.getContext) return false;
var ctx = canvas.getContext('2d');
canvas.onmousemove = function(e) {
if (mouseDownFlag) {
var rect = e.target.getBoundingClientRect();
ctx.beginPath();
ctx.arc(e.clientX - rect.left, e.clientY - rect.top, 3, 0, Math.PI * 2, false);
ctx.fill();
}
}
canvas.onmousedown = function(e) { mouseDownFlag = true; }
canvas.onmouseup = function(e) { mouseDownFlag = false; }
canvas.onmouseout = function(e) { mouseDownFlag = false; }
}
function setImage() {
var canvas = document.getElementById('bbsCanvas');
if (!canvas || !canvas.getContext) return false;
document.getElementById("bbsImage").value = canvas.toDataURL();
}
function clearCanvas() {
var canvas = document.getElementById('bbsCanvas');
if (!canvas || !canvas.getContext) return false;
var ctx = canvas.getContext('2d');
ctx.clearRect(0, 0, 400, 200);
setImage();
}
</script>""" % (user,IP,P,IP,P,user)
cur.execute('SELECT * FROM bbs ORDER BY id DESC')
for i, row in enumerate(cur):
if i>=LIMIT: break
data+=('<div><p>%d <b>%s</b> %s</p><img src="%s" alt="Image" width="400" height="200" style="border:1px solid;" /><p>%s</p></div>' % row)
data+='</body></html>'
start_response('200 OK',[('Content-type','text/html;charset=utf-8')])
return [data.encode('utf-8')]
IP=ipconfig()
httpd=make_server('',P,bbs)
httpd.serve_forever()
2012年1月1日日曜日
2011年12月23日金曜日
円周率を1万桁まで求める
マチンの公式で円周率を1万桁まで求めるプログラムを、Python, Erlang, Haskell, C++で書いてみた。出力結果の最後の数桁はずれているかも。また、実行時間を測ったりもしているが、1万桁程度ならどのコードでも瞬時に求まる。
まずは素直にPythonで。実行時オプションで桁数指定、実行時間測定付き。オプションなしで1万桁まで求める。
pi.py
#!/usr/bin/env python
import sys, time
N = 10**10000
def arctan(m):
global N
c = N
a = b = c / m
m2 = m * m
s = k = 1
while c:
b /= m2
k += 2
c, s = b / k, -s
a += c * s
return a
def main(args):
global N
if len(args) > 1: N = 10**int(args[1])
t1 = time.time()
pi = str((arctan(5) * 4 - arctan(239)) * 4)
t2 = time.time()
print pi[0] + '.' + pi[1:]
print "Time: %f" % (t2 - t1)
if __name__ == "__main__": main(sys.argv)
実行:
$ ./pi.py
次にErlangで求めてみる。
pi.erl
-module(pi).
-export([pi/0]).
pi()->N=e(10,10000),(a(5,N)*4-a(239,N))*4.
e(B,N)->e(B,N,1).
e(_,0,R)->R;
e(B,N,R)->e(B,N-1,R*B).
a(X,N)->a(X,N div X,N div X,N,1,1).
a(_,A,_,0,_,_)->A;
a(M,A,B,_,S,K)->B_=B div(M*M),C=B_ div(K+2),S_=-S,A_=A+C*S_,a(M,A_,B_,C,S_,K+2).
実行:
$ erl
1> c(pi).
2> pi:pi().
そしてHaskellを使う。ソースコードは一行だ。
pi.hs
n=10^10^4;p=a 5 n*4-a 239 n;a m n=t m(n`div`m)(n`div`m)n 1 3;t _ a _ 0 _ _=a;t m a b c s k=t m(a-y*s)x y(-s)(k+2)where x=b`div`m^2;y=x`div`k
実行:
$ ghci
Prelude> :l pi
Prelude> p*4
最後にC++のコード。多倍長演算でちょっと長い。実行時間測定付き。
pi.cpp
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <cmath>
#ifdef _MSC_VER
#include <ctime>
inline double get_time()
{
return static_cast<double>(std::clock()) / CLOCKS_PER_SEC;
}
#else
#include <sys/time.h>
inline double get_time()
{
timeval tv;
gettimeofday(&tv, 0);
return tv.tv_sec + 1e-6 * tv.tv_usec;
}
#endif
typedef unsigned long long number;
const int FIG = 8;
const number BASE = 0xffffffffULL;
const number NN = 100000000ULL;
// C = 1.0 / (4 * log10(2)) / FIG;
#define C (0.830482023722/FIG)
#define M (X/FIG+2)
#define N (static_cast<int>(X*C+2))
template<int X> class PI {
private:
number ans[N];
public:
void add(number* a, number* b);
void sub(number* a, number* b);
void div(number* a, number d);
number intdiv(number x, number y);
void mul(number* a, number d);
void atan(number* a, number m);
number or_numbers(number* a);
bool has_bit(number* a);
void decimal(number* a, number* w);
void print();
void calc();
};
// a <- arctan(1 / m)
template<int X>
void PI<X>::atan(number* a, number m)
{
number b[N] = { 1 };
number c[N] = { 1 };
div(b, m);
std::copy(b, b + N, a);
const number m2 = m * m;
for (int i = 1; has_bit(c); i++) {
div(b, m2);
std::copy(b, b + N, c);
div(c, i * 2 + 1);
if (i & 1) sub(a, c);
else add(a, c);
}
}
// a <- a + b
template<int X>
void PI<X>::add(number* a, number* b)
{
for (int i = 0; i < N; i++) {
number x = a[i] + b[i];
if (x & ~BASE) {
a[i] = x & BASE;
int j;
for (j = i - 1; a[j] == BASE; j--) a[j] = 0;
a[j]++;
} else a[i] = x;
}
}
// a <- a - b (a > b)
template<int X>
void PI<X>::sub(number* a, number* b)
{
for (int i = 0; i < N; i++) {
if (a[i] < b[i]) {
a[i] = (BASE + 1) + a[i] - b[i];
int j;
for (j = i - 1; a[j] == 0; j--) a[j] = BASE;
a[j]--;
} else a[i] -= b[i];
}
}
// a <- a / d (d <= BASE)
template<int X>
void PI<X>::div(number* a, number d)
{
number res = 0;
for (int i = 0; i < N; i++) {
res <<= (FIG * 4);
number x = a[i] + res;
number q = x / d;
a[i] = q;
res = x - q * d;
}
}
template<int X>
number PI<X>::intdiv(number x, number y)
{
const number maxbit = ~((number)-1>>1);
number a = 0;
number b = 0;
int cnt = sizeof(number) * 8;
while (cnt--) {
a <<= 1ULL;
if (x & maxbit) a |= 1ULL;
x <<= 1ULL;
b <<= 1ULL;
if (a >= y) { ++b; a -= y; }
}
return b;
}
// a <- a * d (d <= BASE)
template<int X>
void PI<X>::mul(number* a, number d)
{
number q = 0;
for (int i = N - 1; i >= 0; i--) {
number x = a[i] * d + q;
a[i] = x & BASE;
q = x >> (FIG * 4);
}
}
// OR numbers
template<int X>
number PI<X>::or_numbers(number* a)
{
number ret = 0;
for (int i = 0; i < N; i++) ret |= a[i];
return ret;
}
template<int X>
bool PI<X>::has_bit(number* a)
{
for (int i = 0; i < N; i++) if (a[i]) return true;
return false;
}
// hex to decimal number
template<int X>
void PI<X>::decimal(number* a, number* w)
{
number b[N];
std::copy(a, a + N, b);
w[0] = b[0];
b[0] = 0;
for (int i = 1; i < M; i++) {
mul(b, NN);
w[i] = b[0];
b[0] = 0;
}
}
// print answer
template<int X>
void PI<X>::print()
{
number w[M];
decimal(ans, w);
std::cout << std::setw(4) << std::right << w[0] << ".";
for (int i = 1; i < M; i++) {
std::cout << std::setw(FIG) << std::setfill('0') << w[i] << " ";
if(i % 6 == 0) std::cout << std::endl << " ";
}
std::cout << std::endl;
}
// calculate PI
// Machin's formula: 4 * arctan(1/5) - arctan(1/239) = PI / 4
template<int X>
void PI<X>::calc()
{
number x[N];
atan(ans, 5);
mul(ans, 4);
atan(x, 239);
sub(ans, x);
mul(ans, 4);
}
int main()
{
std::ios_base::sync_with_stdio(false);
PI<10000> pi;
double start_time = get_time();
pi.calc();
double end_time = get_time();
pi.print();
std::cerr << "Time: " << end_time - start_time << std::endl;
return 0;
}
実行:
$ c++ pi.cpp -O3 -o pi
$ ./pi
2011年11月30日水曜日
好きな・嫌いなコンピュータ言語のアンケート
社内で好きなコンピュータ言語と嫌いなコンピュータ言語のアンケートを取ってみた。ことの発端は自分のツイートなのだが、思いがけず賛同を得ることができたので、実際にアンケートを社内で行うことにしたのだ。各自で社内Wikiを編集してもらってもよかったのだけど、プログラムを作る会社でもあるし、投票システムに興味もあったので、Google App Engine (GAE)で作ってみることにした。以前にGAEで掲示板などを作成していたので、それらを流用して時間をかけずにすぐに仕上げることができた。以下のリンクから実際に投票できる。
表示されているコンピュータ言語一覧から、最も好きな言語と最も嫌いな言語をそれぞれ一つずつ、それ以外にも好きな言語や嫌いな言語があれば、それらも選択(複数可)する。投票したい言語が一覧にない場合は「言語の追加」で自由に追加できるようになっている。言語を選択して投票ボタンを押すと、投票数とそのグラフが更新される。表示されるグラフには、「最も好きな・最も嫌いな言語」と「好きな・嫌いな言語」の2つがあり、青のグラフが「好きな言語」、赤のグラフが「嫌いな言語」を表している。「好きな・嫌いな言語」のグラフについては、「最も好きな・最も嫌いな言語」と「好きな・嫌いな言語」を足し合わせた票数となっている。
また、コメントも投稿できるようになっていて、自分の好きな言語を啓蒙するのも良し、嫌いな言語について文句を言うのもアリだ。自由に書き込める。
現時点での社内での票は以下のグラフに示す通りで、幅広く票が入っているように思う。傾向としては、C, C++, C#, ASM, Python, Rubyあたりが人気で、Visual Basic, Java, Perlが不人気みたいだ。
今回作成したコードの簡単な説明を書いておく。上述の「好きな・嫌いなコンピュータ言語」だけでなく、初期変数を設定することで、「好きな・嫌いな動物」「好きな・嫌いな食べ物」などのように、簡単に任意のアンケートを作成できるようになっている。グラフについてはGoogle Chart APIを利用しているが、1,000ピクセルを超えることができないので、表示数を1,000ピクセル以内に抑える処理を入れている。また、頻繁にデータベースにアクセスしないようにmemcacheも利用している。詳細についてはコードを読んで欲しい。
以下に今回作成したGAEのコードを示しておく。
main.py
#!/usr/bin/env python # -*- coding: utf-8 -*- """main.py Copyright (C) 2011 nox""" import os, cgi, re, math, datetime import urllib, Cookie import wsgiref.handlers from google.appengine.api import memcache from google.appengine.ext import webapp from google.appengine.ext import db from google.appengine.ext.webapp import template # 許可する投票タイトル. TITLES = { "test": u"投票テスト", "lang": u"好きな・嫌いなコンピュータ言語のアンケート" } # 表示メッセージ. MESSAGES = { "test": u"投票をお願いします。", "lang": u"最も好きな言語と最も嫌いな言語をそれぞれ1つ、他に好きな言語や嫌いな言語があればそれも選択(複数可)して投票をお願いします。" } # 投票する属性. (タイトル / love / like / hate / dislike) PROPS = { "test": (u"テスト", u"最高", u"良い", u"最低", u"悪い"), "lang": (u"言語", u"最も好き", u"好き", u"最も嫌い", u"嫌い") } # チャートのタイトル. (love and hate / like and dislike [including love and hate]) CHARTS = { "test": (u"最高・最低のテスト", u"良い・悪いテスト"), "lang": (u"最も好きな・最も嫌いな言語", u"好きな・嫌いな言語") } # 最大表示件数 TARGET_LIMIT = 300 COMMENT_LIMIT = 100 class MainPage(webapp.RequestHandler): def get(self): self.response.out.write(u""" <html> <head> <title>Opinion Poll</title> </head> <body> <div><a href="vote">投票システム</a></div> </body> </html>""") # 投票データ. class VoteData(db.Model): id = db.StringProperty(multiline=False) content = db.StringProperty(multiline=False) love = db.IntegerProperty() like = db.IntegerProperty() hate = db.IntegerProperty() dislike = db.IntegerProperty() # コメントデータ. class CommentData(db.Model): id = db.StringProperty(multiline=False) username = db.StringProperty(multiline=False) content = db.TextProperty() date = db.DateTimeProperty(auto_now_add=True) # 投票システム・表示. class Vote(webapp.RequestHandler): def get(self): global TITLES, TARGET_LIMIT id = self.request.get("id") if id not in TITLES.keys(): return # データベースまたはMemcacheからデータの読み込み. cache_votes_key = "votes_%s" % id try: votes = memcache.get(cache_votes_key) except: try: memcache.delete(cache_votes_key) except: pass votes = None if votes is None: votes = VoteData.all().order("content").filter("id =", id).fetch(TARGET_LIMIT) try: memcache.add(cache_votes_key, votes) except: pass cache_comments_key = "comments_%s" % id try: comments = memcache.get(cache_comments_key) except: try: memcache.delete(cache_comments_key) except: pass comments = None if comments is None: comments = CommentData.all().order("-date").filter("id =", id).fetch(COMMENT_LIMIT) try: memcache.add(cache_comments_key, comments) except: pass # チャート表示 (1,000ピクセル以下にデータ数を収める). target2, target4 = [], [] target2_max, target4_max = 0, 0 for v in votes: v.love = int(v.love) if v.love else 0 v.hate = int(v.hate) if v.hate else 0 v.like = int(v.like) if v.like else 0 v.dislike = int(v.dislike) if v.dislike else 0 if v.love or v.hate or v.like or v.dislike: target4.append((v.content, v.love + v.like, v.hate + v.dislike)) target4_max = max(target4_max, v.love + v.like, v.hate + v.dislike) if v.love or v.hate: target2.append((v.content, v.love, v.hate)) target2_max = max(target2_max, v.love, v.hate) target2_max = (target2_max - 1) / 10 if target2_max > 0 else 0 target2_max = (target2_max + 1) * 10 target4_max = (target4_max - 1) / 10 if target4_max > 0 else 0 target4_max = (target4_max + 1) * 10 n = 0 while True: target2 = filter(lambda t: max(t[1:]) > n, target2) if len(target2) <= 16: break n += 1 if target2: text = "" if n > 0: text = u" (%d票以上)" % (n + 1) fig = 10**(int(math.log10(target2_max))) if target2_max / fig < 5: fig /= 2 chart_scale = "|".join([str(i) for i in range(0, target2_max + 1, fig)]) chart_target2_url = u"http://chart.apis.google.com/chart?chs=%dx200&chd=t:%s|%s&chds=0,%d&cht=bvg&chco=4d89f9,f9894d&chxt=y,x&chxl=0:|%s|1:|%s&chtt=%s%s" % (len(target2) * 60 + 30, ",".join([str(t[1]) for t in target2]), ",".join([str(t[2]) for t in target2]), target2_max, chart_scale, "|".join([urllib.quote(t[0].encode('utf-8')) for t in target2]), CHARTS[id][0], text) else: chart_target2_url = "" n = 0 while True: target4 = filter(lambda t: max(t[1:]) > n, target4) if len(target4) <= 16: break n += 1 if target4: text = "" if n > 0: text = u" (%d票以上)" % (n + 1) fig = 10**(int(math.log10(target4_max))) if target4_max / fig < 5: fig /= 2 chart_scale = "|".join([str(i) for i in range(0, target4_max + 1, fig)]) chart_target4_url = u"http://chart.apis.google.com/chart?chs=%dx200&chd=t:%s|%s&chds=0,%d&cht=bvg&chco=4d89f9,f9894d&chxt=y,x&chxl=0:|%s|1:|%s&chtt=%s%s" % (len(target4) * 60 + 30, ",".join([str(t[1]) for t in target4]), ",".join([str(t[2]) for t in target4]), target4_max, chart_scale, "|".join([urllib.quote(t[0].encode('utf-8')) for t in target4]), CHARTS[id][1], text) else: chart_target4_url = "" # 挨拶. H = (datetime.datetime.utcnow() + datetime.timedelta(hours=9)).hour if H >= 5 and H < 11: greeting = u"おはようございます。" elif H >= 11 and H < 17: greeting = u"こんにちは。" elif H >= 17 or H < 5: greeting = u"こんばんは。" else: greeting = u"こんにちは。" greeting += MESSAGES[id] # ウェブ上で表示させるデータ. template_values = { "id": id, "title": TITLES[id], "target": PROPS[id][0], "love": PROPS[id][1], "like": PROPS[id][2], "hate": PROPS[id][3], "dislike": PROPS[id][4], "votes": votes, "comments": comments, "chart_target2_url": chart_target2_url, "chart_target4_url": chart_target4_url, "greeting": greeting, } path = os.path.join(os.path.dirname(__file__), "vote.html") self.response.out.write(template.render(path, template_values)) # 投票ターゲットの追加. class AddVote(webapp.RequestHandler): def post(self): vote = VoteData() vote.id = self.request.get("id") try: try: memcache.delete("votes_%s" % vote.id) except: pass except: pass if not self.request.get("content").strip(): self.redirect("/vote?id=%s" % vote.id) return content = cgi.escape(self.request.get("content").strip()[:100]) v = VoteData.all().filter("id =", vote.id).filter("content =", content).get() if not v: vote.content = content vote.put() self.redirect("/vote?id=%s" % vote.id) # 投票. class PostVote(webapp.RequestHandler): def post(self): vote = VoteData() vote.id = self.request.get("id") love = self.request.get_all("love") like = self.request.get_all("like") hate = self.request.get_all("hate") dislike = self.request.get_all("dislike") try: try: memcache.delete("votes_%s" % vote.id) except: pass except: pass if not (love or like or hate or dislike): self.redirect("/vote?id=%s" % vote.id) return for d in love: v = VoteData.all().filter("id =", vote.id).filter("content =", d).get() if not v.love: v.love = 1 else: v.love = int(v.love) + 1 v.put() for d in like: v = VoteData.all().filter("id =", vote.id).filter("content =", d).get() if not v.like: v.like = 1 else: v.like = int(v.like) + 1 v.put() for d in hate: v = VoteData.all().filter("id =", vote.id).filter("content =", d).get() if not v.hate: v.hate = 1 else: v.hate = int(v.hate) + 1 v.put() for d in dislike: v = VoteData.all().filter("id =", vote.id).filter("content =", d).get() if not v.dislike: v.dislike = 1 else: v.dislike = int(v.dislike) + 1 v.put() self.redirect("/vote?id=%s" % vote.id) # コメントの投稿. class PostComment(webapp.RequestHandler): def insert_whitespace(self, matchobj): return re.sub("(?<!<br) ", " ", matchobj.group(0).replace("\t", " ")) def post(self): global TARGET_LIMIT comment = CommentData() comment.id = self.request.get("id") try: try: memcache.delete("comments_%s" % comment.id) except: pass except: pass if not self.request.get("content").strip(): self.redirect("/vote?id=%s" % comment.id) return if self.request.get("username").strip(): comment.username = cgi.escape(self.request.get("username").strip()[:255]) else: comment.username = "Anonymous" comment.content = cgi.escape(self.request.get("content").strip()[:2000]) comment.content = re.sub("\r\n|\r|\n", "<br />", comment.content) comment.content = re.sub("\[link\](.*?)\[/link\]", "<a href=\"\g<1>\">\g<1></a>", comment.content) comment.content = re.sub("\[code\](.*?)\[/code\]", self.insert_whitespace, comment.content) comment.content = re.sub("\[code\](.*?)\[/code\]", "<span style=\"font-family:courier new,monospace;font-size:85%;white-space:pre;color:#0000ff;\">\g<1></span>", comment.content) comment.date = comment.date + datetime.timedelta(hours=9) comment.put() self.redirect("/vote?id=%s" % comment.id) def main(): application = webapp.WSGIApplication([("/", MainPage), ("/vote", Vote), ("/vote/add", AddVote), ("/vote/post", PostVote), ("/vote/comment", PostComment), ], debug=False) wsgiref.handlers.CGIHandler().run(application) if __name__ == "__main__": main()
vote.html
<?xml version="1.0" encoding="UTF-8"?> <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1//EN" "http://www.w3.org/TR/xhtml11/DTD/xhtml11.dtd"> <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="ja"> <head> <meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> <link type="text/css" rel="stylesheet" href="/stylesheets/main.css" /> <title>{{ title }}</title> </head> <body> <script type="text/javascript"> function resize_textarea(ev){ var textarea = ev.target || ev.srcElement; var value = textarea.value; var lines = 1; for (var i = 0, l = value.length; i < l; i++) if (value.charAt(i) == '\n') lines++; if (lines > 5) { if (lines > 20) textarea.setAttribute("rows", 20); else textarea.setAttribute("rows", lines); } else textarea.setAttribute("rows", 5); } </script> <div style="font-weight: bold; font-size: large; padding-bottom: 20px;">{{ title }}</div> <div> <p>{{ greeting }}</p> </div> <form action="/vote/add" method="post"> <input type="text" name="content" size="20" maxlength="100"></input> <input type="submit" value="{{ target }}の追加" /> <input type="hidden" name="id" value="{{ id }}"> </form> <form action="/vote/post" method="post"> <div> <table border="0" width="400"> <th>{{ target }}</th><th>{{ love }}</th><th>{{ like }}</th><th>{{ hate }}</th><th>{{ dislike }}</th> {% for vote in votes %} <tr> <td><a href="http://ja.wikipedia.org/wiki/{{ vote.content }}">{{ vote.content }}</a></td> <td align="center"><input type="radio" name="love" value="{{ vote.content }}" /> {% if vote.love %} {{ vote.love }} {% else %} 0 {% endif %}</td> <td align="center"><input type="checkbox" name="like" value="{{ vote.content }}" /> {% if vote.like %} {{ vote.like }} {% else %} 0 {% endif %}</td> <td align="center"><input type="radio" name="hate" value="{{ vote.content }}" /> {% if vote.hate %} {{ vote.hate }} {% else %} 0 {% endif %}</td> <td align="center"><input type="checkbox" name="dislike" value="{{ vote.content }}" /> {% if vote.dislike %} {{ vote.dislike }} {% else %} 0 {% endif %}</td> </tr> {% endfor %} <tr align=center> <td colspan=5><input type="submit" value="投票" style="width: 100px; height: 25px" /></td> </tr> </table> <input type="hidden" name="id" value="{{ id }}"> </div> </form><br /> <div> {% if chart_target2_url %} <p><img src="{{ chart_target2_url }}" /></p> {% endif %} {% if chart_target4_url %} <p><img src="{{ chart_target4_url }}" /></p> {% endif %} </div> コメント記入: <span style="font-size:70%;color:#666666;">リンク: [link]～[/link], ソースコード: [code]～[/code], 1,000文字まで</span><br /> <form action="/vote/comment" method="post"> <div><textarea name="content" rows="5" cols="60" onkeyup="resize_textarea(event)"></textarea></div> 名前: <input type="text" name="username" size="31" maxlength="255" value="{{ username }}"></input> <input type="submit" value="コメント" /> <span style="font-size:70%;color:#666666;">無記名で匿名投稿</span> <input type="hidden" name="id" value="{{ id }}"> </form> <div> {% for comment in comments %} {{ comment.username }} [{{ comment.date.year }}年{{ comment.date.month }}月{{ comment.date.day }}日{{ comment.date.hour }}時{{ comment.date.minute }}分{{ comment.date.second }}秒]: <blockquote>{{ comment.content }}</blockquote> {% endfor %} </div> <!-- ここより下は変更しないでください。不具合やご意見などがありましたら、下記のnoxのリンク先で報告していただけると嬉しいです。 --> <div><img src="http://code.google.com/appengine/images/appengine-silver-120x30.gif" alt="Powered by Google App Engine" /> <span style="font-size:70%;color:#666666;"><a href="http://handasse.blogspot.com/">ソースコードおよびその解説</a></span><br /><span style="font-size:70%;color:#666666;">Copyright © 2011 <a href="http://handasse.blogspot.com/">nox</a>. All rights reserved.</span></div> </body> </html>
2011年10月19日水曜日
コンピュータにラストワンを解かせてみた
先日、地元のまつりがあった。そこの会場となっている区民館では子供用ボーリングやプラ板遊びなどの催しが開かれていたのだが、その中の一つに「ひとりぼっち」というゲームを行なっているのを見つけたので、一緒にいた小4の自分の娘に「このゲーム面白そうだね」と言ったところ、「ラストワンでしょ」と返され、よく知っているゲームとのことだった。遊び方は次の通り。全部で33のマスがあり、真ん中以外の32個のマスにコマを置く。コマは縦・横方向に隣のコマを一つ飛び越えて進むことができる。飛び越えられたコマは取り除かれる。これを繰り返して、最後のコマを真ん中のマスに置けばクリアとなる。32コマを使うステージ以外にもT字型やピラミッド型にコマを配置するなど、さまざまなバリエーションがある。
ネットで調べたところ、1985年にバンダイから同名の電子ゲームが出ていたようだ。リンク先では動画もあるのでそれを見てもらえば遊び方は一目瞭然だろう。普通に解いても面白そうだが、コンピュータを使って解かせたくなったので、早速プログラムを作ってみた。
余談だが、子供のころ、コンピュータは何でも解いてくれる魔法の箱だと思っていた。円周率を何桁まで解いたと聞いて、ワクワクしたものだ。ラストワンのようなゲームをコンピュータに解かせたいという気持ちは、子供のころのコンピュータに対する憧れから来るのかもしれない。
閑話休題。ラストワンは33マスあり、コマが「ある」「ない」の2状態のみを持つので、33ビットを使って盤面の状態を表すことができる。そして、ある状態から取り得る次の状態をすべて求めて、それを優先順位付きキューに入れて探索した。キューが空になるか、ゴールの状態が見つかるまで探索を行う。優先順位付きキューにしたのは探索を短縮するためにヒューリスティックな関数でスコアを決めようと思ったからなんだけど、32コマあるステージも含め、どれも一瞬で解けてしまうので、どうやら高度なヒューリスティック関数を用意しなくても問題なかったらしい。
上図のように初期配置としてT字型にコマを置く場合、以下に示すような入力データを用意して、これをプログラムに渡す。o (オー)でコマのあるマスを表し、. (ドット)で空のマスを表している。
lastone.dat
...
...
.......
..ooo..
...o...
.o.
...
入力データのファイル名が lastone.dat であれば、以下のように実行すればよい。
% ./lastone lastone.dat
実行結果は以下のようになる。
...
...
.......
..ooo..
...o...
.o.
...
...
...
.......
.o..o..
...o...
.o.
...
...
...
.......
.o.oo..
.......
...
...
...
...
.......
.oo....
.......
...
...
...
...
.......
...o...
.......
...
...
最後に作成したコードを示す。
lastone.cpp
// Last One solver by nox, 14 Oct. 2011
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <queue>
#include <algorithm>
#include <utility>
using namespace std;
typedef unsigned long long state_t;
typedef int score_t;
const int N = 33; // number of cells
class LastOne {
private:
vector<vector<pair<int, int> > > D; // table of directions
state_t start, goal;
vector<state_t> ans; // solution
public:
LastOne(const string datafile="");
~LastOne() { }
void read(const string datafile);
int movable(state_t state);
void next(state_t state, vector<state_t>& neighbors);
bool search();
void show(state_t state);
void answer();
};
// initialize and make table
LastOne::LastOne(const string datafile)
{
if (!datafile.empty()) read(datafile);
vector<pair<int, int> > v;
v.push_back(make_pair(1, 2));
v.push_back(make_pair(3, 8));
D.push_back(v);
v.clear();
v.push_back(make_pair(4, 9));
D.push_back(v);
v.clear();
v.push_back(make_pair(1, 0));
v.push_back(make_pair(5, 10));
D.push_back(v);
v.clear();
v.push_back(make_pair(4, 5));
v.push_back(make_pair(8, 15));
D.push_back(v);
v.clear();
v.push_back(make_pair(9, 16));
D.push_back(v);
v.clear();
v.push_back(make_pair(4, 3));
v.push_back(make_pair(10, 17));
D.push_back(v);
v.clear();
v.push_back(make_pair(7, 8));
v.push_back(make_pair(13, 20));
D.push_back(v);
v.clear();
v.push_back(make_pair(8, 9));
v.push_back(make_pair(14, 21));
D.push_back(v);
v.clear();
v.push_back(make_pair(3, 0));
v.push_back(make_pair(7, 6));
v.push_back(make_pair(9, 10));
v.push_back(make_pair(15, 22));
D.push_back(v);
v.clear();
v.push_back(make_pair(4, 1));
v.push_back(make_pair(8, 7));
v.push_back(make_pair(10, 11));
v.push_back(make_pair(16, 23));
D.push_back(v);
v.clear();
v.push_back(make_pair(5, 2));
v.push_back(make_pair(9, 8));
v.push_back(make_pair(11, 12));
v.push_back(make_pair(17, 24));
D.push_back(v);
v.clear();
v.push_back(make_pair(10, 9));
v.push_back(make_pair(18, 25));
D.push_back(v);
v.clear();
v.push_back(make_pair(11, 10));
v.push_back(make_pair(19, 26));
D.push_back(v);
v.clear();
v.push_back(make_pair(14, 15));
D.push_back(v);
v.clear();
v.push_back(make_pair(15, 16));
D.push_back(v);
v.clear();
v.push_back(make_pair(8, 3));
v.push_back(make_pair(14, 13));
v.push_back(make_pair(16, 17));
v.push_back(make_pair(22, 27));
D.push_back(v);
v.clear();
v.push_back(make_pair(9, 4));
v.push_back(make_pair(15, 14));
v.push_back(make_pair(17, 18));
v.push_back(make_pair(23, 28));
D.push_back(v);
v.clear();
v.push_back(make_pair(10, 5));
v.push_back(make_pair(16, 15));
v.push_back(make_pair(18, 19));
v.push_back(make_pair(24, 29));
D.push_back(v);
v.clear();
v.push_back(make_pair(17, 16));
D.push_back(v);
v.clear();
v.push_back(make_pair(18, 17));
D.push_back(v);
v.clear();
v.push_back(make_pair(13, 6));
v.push_back(make_pair(21, 22));
D.push_back(v);
v.clear();
v.push_back(make_pair(14, 7));
v.push_back(make_pair(22, 23));
D.push_back(v);
v.clear();
v.push_back(make_pair(15, 8));
v.push_back(make_pair(21, 20));
v.push_back(make_pair(23, 24));
v.push_back(make_pair(27, 30));
D.push_back(v);
v.clear();
v.push_back(make_pair(16, 9));
v.push_back(make_pair(22, 21));
v.push_back(make_pair(24, 25));
v.push_back(make_pair(28, 31));
D.push_back(v);
v.clear();
v.push_back(make_pair(17, 10));
v.push_back(make_pair(23, 22));
v.push_back(make_pair(25, 26));
v.push_back(make_pair(29, 32));
D.push_back(v);
v.clear();
v.push_back(make_pair(18, 11));
v.push_back(make_pair(24, 23));
D.push_back(v);
v.clear();
v.push_back(make_pair(19, 12));
v.push_back(make_pair(25, 24));
D.push_back(v);
v.clear();
v.push_back(make_pair(22, 15));
v.push_back(make_pair(28, 29));
D.push_back(v);
v.clear();
v.push_back(make_pair(23, 16));
D.push_back(v);
v.clear();
v.push_back(make_pair(24, 17));
v.push_back(make_pair(28, 27));
D.push_back(v);
v.clear();
v.push_back(make_pair(27, 22));
v.push_back(make_pair(31, 32));
D.push_back(v);
v.clear();
v.push_back(make_pair(28, 23));
D.push_back(v);
v.clear();
v.push_back(make_pair(29, 24));
v.push_back(make_pair(31, 20));
D.push_back(v);
}
// read LastOne data
void LastOne::read(const string datafile)
{
fstream fs(datafile.c_str(), ios_base::in);
string line;
goal = 0x10000;
start = 0ULL;
while (getline(fs, line)) {
for (string::iterator p = line.begin(); p != line.end(); ++p) {
if (*p == 'o') { start <<= 1; start |= 1ULL; }
else if (*p == '.') start <<= 1;
}
}
fs.close();
}
// number of movable pieces
int LastOne::movable(state_t state)
{
int s = 0;
for (int i = 0; i < N; i++)
if ((state >> i) & 1)
for (vector<pair<int, int> >::iterator p = D[i].begin(); p != D[i].end(); ++p)
s += ((state >> p->first) & 1) & !((state >> p->second) & 1);
return s;
}
// next states
void LastOne::next(state_t state, vector<state_t>& neighbors)
{
neighbors.clear();
for (int i = 0; i < N; i++)
if ((state >> i) & 1)
for (vector<pair<int, int> >::iterator p = D[i].begin(); p != D[i].end(); ++p)
if (((state >> p->first) & 1) & !((state >> p->second) & 1))
neighbors.push_back(state ^ (1ULL << p->first) ^ (1ULL << p->second) ^ (1ULL << i));
}
// search solution
bool LastOne::search()
{
vector<state_t> checked(1, start);
priority_queue<pair<score_t, vector<state_t> > > queue;
queue.push(make_pair(1, checked));
vector<state_t> neighbors;
while (!queue.empty()) {
pair<score_t, vector<state_t> > q(queue.top());
queue.pop();
state_t last = q.second.back();
if (last == goal) {
ans = q.second;
return true;
}
next(last, neighbors);
for (vector<state_t>::iterator p = neighbors.begin(); p != neighbors.end(); ++p) {
if (binary_search(checked.begin(), checked.end(), *p)) continue;
if (movable(*p) == 0 && *p != goal) continue;
checked.insert(upper_bound(checked.begin(), checked.end(), *p), *p);
vector<state_t> path(q.second);
path.push_back(*p);
queue.push(make_pair(static_cast<score_t>(path.size()), path));
}
}
return false;
}
// show state
void LastOne::show(state_t state)
{
string line;
while (state) {
if (state & 1ULL) line.push_back('o');
else line.push_back('.');
state >>= 1;
}
int n = N - line.size();
for (int i = 0; i < n; i++) line.push_back('.');
reverse(line.begin(), line.end());
cout << " " << line.substr(0, 3) << endl;
cout << " " << line.substr(3, 3) << endl;
cout << line.substr(6, 7) << endl;
cout << line.substr(13, 7) << endl;
cout << line.substr(20, 7) << endl;
cout << " " << line.substr(27, 3) << endl;
cout << " " << line.substr(30, 3) << endl;
}
// show solution
void LastOne::answer()
{
for (vector<state_t>::iterator p = ans.begin(); p != ans.end(); ++p) {
show(*p);
cout << endl;
}
}
int main(int argc, char* argv[])
{
if (argc < 2) {
cerr << "Usage: " << argv[0] << " datafile" << endl;
exit(1);
}
LastOne lastone(argv[1]);
if (lastone.search()) lastone.answer();
else cout << "No solution." << endl;
return 0;
}
