良いもの。悪いもの。

マチンの公式 で円周率を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)->...

社内で好きなコンピュータ言語と嫌いなコンピュータ言語のアンケートを取ってみた。ことの発端は 自分のツイート なのだが、思いがけず賛同を得ることができたので、実際にアンケートを社内で行うことにしたのだ。各自で社内Wikiを編集してもらってもよかったのだけど、プログラムを作る会社でもあるし、投票システムに興味もあったので、Google App Engine (GAE)で作ってみることにした。以前にGAEで掲示板などを作成していたので、それらを流用して時間をかけずにすぐに仕上げることができた。以下のリンクから実際に投票できる。 好きな・嫌いなコンピュータ言語のアンケート 表示されているコンピュータ言語一覧から、最も好きな言語と最も嫌いな言語をそれぞれ一つずつ、それ以外にも好きな言語や嫌いな言語があれば、それらも選択(複数可)する。投票したい言語が一覧にない場合は「言語の追加」で自由に追加できるようになっている。言語を選択して投票ボタンを押すと、投票数とそのグラフが更新される。表示されるグラフには、「最も好きな・最も嫌いな言語」と「好きな・嫌いな言語」の2つがあり、青のグラフが「好きな言語」、赤のグラフが「嫌いな言語」を表している。「好きな・嫌いな言語」のグラフについては、「最も好きな・最も嫌いな言語」と「好きな・嫌いな言語」を足し合わせた票数となっている。 また、コメントも投稿できるようになっていて、自分の好きな言語を啓蒙するのも良し、嫌いな言語について文句を言うのもアリだ。自由に書き込める。 現時点での社内での票は以下のグラフに示す通りで、幅広く票が入っているように思う。傾向としては、C, C++, C#, ASM, Python, Rubyあたりが人気で、Visual Basic, Java, Perlが不人気みたいだ。 今回作成したコードの簡単な説明を書いておく。上述の「好きな・嫌いなコンピュータ言語」だけでなく、初期変数を設定することで、「好きな・嫌いな動物」「好きな・嫌いな食べ物」などのように、簡単に任意のアンケートを作成できるようになっている。グラフについては Google Chart API を利用しているが、1,000ピクセルを超えることができないので、表示数を1,000ピクセル以内に抑える処理を入れている。また、頻繁にデータベ...

先日、地元のまつりがあった。そこの会場となっている区民館では子供用ボーリングやプラ板遊びなどの催しが開かれていたのだが、その中の一つに「ひとりぼっち」というゲームを行なっているのを見つけたので、一緒にいた小4の自分の娘に「このゲーム面白そうだね」と言ったところ、「ラストワンでしょ」と返され、よく知っているゲームとのことだった。遊び方は次の通り。全部で33のマスがあり、真ん中以外の32個のマスにコマを置く。コマは縦・横方向に隣のコマを一つ飛び越えて進むことができる。飛び越えられたコマは取り除かれる。これを繰り返して、最後のコマを真ん中のマスに置けばクリアとなる。32コマを使うステージ以外にもT字型やピラミッド型にコマを配置するなど、さまざまなバリエーションがある。 ネットで調べたところ、 1985年にバンダイから同名の電子ゲームが出ていた ようだ。リンク先では動画もあるのでそれを見てもらえば遊び方は一目瞭然だろう。普通に解いても面白そうだが、コンピュータを使って解かせたくなったので、早速プログラムを作ってみた。 余談だが、子供のころ、コンピュータは何でも解いてくれる魔法の箱だと思っていた。円周率を何桁まで解いたと聞いて、ワクワクしたものだ。ラストワンのようなゲームをコンピュータに解かせたいという気持ちは、子供のころのコンピュータに対する憧れから来るのかもしれない。 閑話休題。ラストワンは33マスあり、コマが「ある」「ない」の2状態のみを持つので、33ビットを使って盤面の状態を表すことができる。そして、ある状態から取り得る次の状態をすべて求めて、それを優先順位付きキューに入れて探索した。キューが空になるか、ゴールの状態が見つかるまで探索を行う。優先順位付きキューにしたのは探索を短縮するためにヒューリスティックな関数でスコアを決めようと思ったからなんだけど、32コマあるステージも含め、どれも一瞬で解けてしまうので、どうやら高度なヒューリスティック関数を用意しなくても問題なかったらしい。 上図のように初期配置としてT字型にコマを置く場合、以下に示すような入力データを用意して、これをプログラムに渡す。 o (オー)でコマのあるマスを表し、 . (ドット)で空のマスを表している。 lastone.dat ... ... ..........

現在のオフィスでは、エントランスからどの通路を通るのが最も近いのか、未だに話題になることがある。エントランスからオフィスにつながるエレベータまでの道が二通りあり、どちらも似たような距離なので同じオフィスに行く人達なのにそこで別々になってしまうこともしばしば。左の図が問題のフロアなんだけど、エントランスから青丸で示したエレベータまで行く必要がある。だったらコンピュータに解かせてしまえばいいじゃない、ということでプログラムを作って最短経路と正確な距離を調べてみた。定規で測ったほうが早いという意見は聞こえません。 以前、 Python: 画像で与えられた迷路に対し2点間の最短経路を求める というブログ記事で、画像から直接最短距離を求めるプログラムを書いたことがあるので、まずはそれを使って調べてみた。図面の邪魔な文字だけ消して、あとはスタート地点とゴール地点の座標を引数で指定するだけだから簡単なんだけど、これだと問題があることに気がついた。というのはCCL (connected component labeling; 連結成分ラベリング)で移動できる領域を調べて、それをA*で解いているのだけど、隣接するノード(ピクセル)を上下左右斜めの8方向のみで調べているから正しい距離を出すことができないのだ。実際に作成した図は以下のようにギザギザになってしまい最短ルートを通らない。なので通り道の分岐や方向が変化する箇所をノードにしてそれをA*で求めることにした。 全部で10個のノード(画像のピクセル座標で指定)を作って調べたところ、以下のような結果になった。左側を通る通路の距離が324.53、右側を通る通路の距離が322.51となり僅かに右側の通路のほうが近かった。ただし、右側の方は途中でフロア案内と少し重なるし、6つあるエレベータのどれを使うかによっても結果が変わってしまうぐらいなので、実質的には同じ距離といってしまっても差し支えないと思う。 結局、それぞれの距離にはほとんど差がないので、結論としては「好きな方を使え」ということになるだろう。自分は広い通路の方が好みなので右側の通路を利用している。最後にソースコードとその使い方を示しておく。 nodes.datに、一行を"ピクセルX座標 ピクセルY座標 隣接ノード番号リスト"としてノード数の...

しばらく前から goo.gl のAPI (Google URL Shortener API)が使えるようになっていたことは知っていたけど、実際に使ったことなかったので試しに使ってみた。import文を含めなければ3行で短縮URLの作成・復元の自動判別を行うことができた。条件演算子などを使えば1行にまとめられる。簡単。 ただし、正規表現は適当だし、エラー処理はしてないので、ちゃんと使う場合はその辺の処理を行う必要あり。あと、本格的に使うのであれば APIキーを取得 すべき。また、 Google APIs Client Library for Python を 使う方法 もある。 http://handasse.blogspot.com/ を短縮URLに変換: % googl.py http://handasse.blogspot.com/ http://goo.gl/RiIiD http://goo.gl/RiIiD を元のURLに戻す: % googl.py http://goo.gl/RiIiD http://handasse.blogspot.com/ googl.py import sys,re from urllib2 import urlopen as U, Request as R from json import loads as J API,URL="https://www.googleapis.com/urlshortener/v1/url",sys.argv[1] if re.match('http://goo\.gl/.+',URL):print J(U(API+'?shortUrl=%s'%URL).read())['longUrl'] else:print J(U(R(API,'{"longUrl":"%s"}'%URL,{'Content-Type':'application/json'})).read())['id'] 追記 (2015/3/26): 現在、API Keyを利用しないとエラーが返ってくるようなので、以下のように修...

ちょっと出先で手持ち無沙汰になったとき、なにかコードが書きたくなることってあると思う。その時に書くプログラムは何でもよくて、言語も選ばない。でも、駅で電車を待っている間とか、注文したメニューが来るまでとか、ノートPCを引っ張り出してまでやりたくはない。そんなとき、Android端末を片手に SL4A でコードを書くっていうのは丁度よい。 というわけで、ちょっと外出したときの暇な時間にSL4AのPythonで書いたコードが以下のコード。本当にくだらないプログラムなわけだけど、暇潰しにはなった。何をするプログラムか興味のある方は読み解いて欲しい。暇潰しぐらいにしかならないけど。まあ、勘のよい人ならファイル名だけでわかってしまうかも。 Android端末を持っている方なら、SL4Aのメニューを開いてAddを選択し、Scan Barcodeで冒頭のQRコードを読み込むことでソースコードを取得できる。SL4Aで書いたプログラムだけど、Android端末以外でもPythonの動作する環境なら大抵は実行できると思う。 Zm9yIGkgaW4gcmFuZ2UoMSwxMDEpOnByaW50J0ZpenonKigxLWklMykrJ0J1enonKigxLWklNSlvciBp.py exec(r"""%s'''gzga*%%%%,hmkl*ocr*nco`fc"a8ajp*mpf*a+\\0+.p%%%s'kormpv"mq.`cqg469gzga*`cqg46,fgamfgqvpkle*mq,rcvj,`cqglcog*]]dkng]]+Y8/1_++')))%%+++''')))"""%(lambda _:(_,_))("exec(''.join(map(lambda _:chr(ord(_)^2),"))

最近、本格的にScalaを使用するつもりで環境を整えた。ビルドツールとして Simple Build Tool (sbt) を利用し、エディタは ensime を入れたEmacsにしている。利用するに当たって .emacsに必要な記述を加えておく こと。初めてのプロジェクトの場合、適当に作成したディレクトリ内でsbtを起動してプロジェクトの雛形を作った後、Emacsを立ち上げて M-x ensime-config-gen で下準備をする。ここまでは初回のみの作業となる。その後は M-x ensime を実行して、src/main/scala/内でコードを書くだけ。これで、コーディング中にタブで補完してくれるし、文法などが間違っている場合に赤でハイライトもしてくれる。コンパイルや実行などは、C-c C-v s とすればsbtがEmacs上で立ち上がるので compile や run するだけ。簡単。 Scalaと言えばアクターでしょ、やっぱり。ということで今回はアクターモデルを使用したコードを書いてみた。自分はScalaをいじるよりも前にErlangをいじっていたので余計にそう考えてしまうのかもしれない。それでもマルチコアやらメニーコアが溢れているこの時代、並行処理の重要性は誰もが認識しているはず。 書いたコードはマンデルブロ集合の計算だ。毎回飽きもせずにまたそれかと思うかもしれないが、並列処理させるのが簡単だし、いつも同じ内容にしておけば比較も容易。なによりもマンデルブロ集合の画像はともて魅力的だと自分では思っているので、今回もマンデルブロ集合だ。冒頭に作成した画像を載せたけど、やっぱり美しいと思う。 今回のアクターモデルでは、scala.actors.Actorから派生させずに、scala.actors.Actor.actorを使用した。クラスを作るほど大行ではないし、その方がスッキリするからだ。そして、計算のループ部分は末尾再帰とパターンマッチを使っている。関数型言語として使うならやっぱりその方がしっくりくる。ただ、ちゃんと末尾再帰させるように、final修飾子をつけるなど、気をつけること。それと@tailrecアノテーション(scala.annotation.tailrec)は常に付けたほうが良い。@tailrecアノテーションを末尾再帰させたい...

以前、 複数のGPUでマンデルブロ集合を並列計算 というブログ記事を書いた。このときはTeslaが1枚で7.45秒、Teslaが2枚で4.26秒で計算できたと報告したが、マルチGPU化は多少面倒だったし、ある程度のオーバーヘッドが掛かっていた。 しかし、 CUDA 4.0 がリリースされてから状況が一変した。自由にcudaSetDeviceが呼び出せるようになり、またUVA (Unified Virtual Addressing)を利用することにより、それはそれは簡単にマルチGPU化や、ホスト・複数GPU間での一元的なメモリの参照を実現できることになった。 というわけで、CUDA 4.0を使ったマンデルブロ集合のプログラムを以下に示す。シングルGPUからマルチGPUへ変更した部分については赤字で示してある。ただし、簡便のために、GPUデバイスは2枚で決め打ちとし、マンデルブロ集合の幅は2の倍数のピクセルとなるようにしてある。また、UVAは使わなくても良かったのだが、使い方を示すために利用している。 CUDA 4.0によるマルチGPU化によりマンデルブロ集合計算の実行時間(計算条件は 前回 と同じ)は 7.45秒から3.73秒 となった。ほぼ2倍の速度である。シングルGPUからの修正箇所はほんの僅か。カーネル関数に至っては一文字たりとも変更していない。良い時代になったものだ。 mandelbrot_multigpu.cu // madelbrot using multi GPU devices by nox, 2011.06.06 // nvcc -lcutil_x86_64 -arch sm_13 -use_fast_math -prec-sqrt=false -keep -L ~/NVIDIA_GPU_Computing_SDK/C/lib -I ~/NVIDIA_GPU_Computing_SDK/C/common/inc -g mandelbrot_multigpu.cu -o mandelbrot_multigpu #include <iostream> #include <fstream> #include "cutil_inline.h" using namespace std; ...

覚書。Base64エンコード・デコードを行う関数とURLエンコードを行う関数をどんな環境でもすぐに参照できるように書いておく。ただし、短いコードのままにしておきたかったので、まったくもって安全ではない。入力前にチェックを入れるとか、関数内で制限を掛けるとか、引数に最大文字数を渡すとか、適当な処理を入れること。よく分からない場合は使わないのが吉。 #include <iostream> #include <iomanip> #include <sstream> #include <algorithm> #include <cstring> using namespace std; int base64encode(const char* p, char* buf) { const char b64[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/="; const int in_size = strlen(p); if (in_size == 0 || p == buf) return 0; int cnt = 0; for (int i = 0; i < (in_size - 1) / 3 + 1; i++) { const unsigned char* c = reinterpret_cast<const unsigned char*>(&p[i*3]); for (int j = 0; j < 4; j++) if (in_size >= i * 3 + j) buf[cnt++] = b64[(c[j-1]<<(6-j*2) | c[j]>>(j*2+2)) & 0x3f]; else buf[cnt++] = '='; } buf[cnt] = '\0'; return cnt; } int base64decode(const c...

今回は興味を持ってもらえる人が限定される記事だと思うけど、以前作ったコードを腐らせてしまうのも勿体無いので敢えて晒してみる。 生体内のタンパク質構造をX線回折やNMRで解析したデータの多くは PDB (Protein Data Bank)データ として登録される。以下に示すように原子タグ、原子の通し番号、原子名、アミノ酸名、アミノ酸の通し番号、原子の座標などがテキストで書かれている。そして、機能の解明などのために 分子動力学(MD; molecular dynamics) 計算プログラムなどでそれらのPDBデータは利用される。MD計算プログラムは世に沢山あるが、自分がよく使っていたのは Amber というソフトウェアだ。 PDBデータの一部 ATOM 1439 N ILE 212 52.408 -3.243 59.013 1.00 40.03 N ATOM 1440 CA ILE 212 52.511 -3.151 57.556 1.00 33.72 C ATOM 1441 C ILE 212 51.464 -2.243 56.911 1.00 34.65 C ATOM 1442 O ILE 212 51.733 -1.619 55.886 1.00 38.03 O ATOM 1443 CB ILE 212 52.482 -4.572 56.921 1.00 38.88 C ATOM 1444 CG1 ILE 212 53.615 -4.703 55.919 1.00 40.75 C ATOM 1445 CG2 ILE 212 51.129 -4.901 56.281 1.00 39.55 C ATOM 1446 CD1 ILE 212 54.931 -4.446 56.550 1.00 46.31 C ATOM 1447 N HIS ...

東日本大震災(東北関東大震災)により被害を受けられた方々には心よりお見舞い申し上げます。 以前、 日本全国コンビニ店舗分布地図 を Processing で作ったが、今回の地震の震源とその大きさについて視覚的に表した地図を作成した。少しでも何かの役に立てば良いのだが。 震源とその大きさを 日本地図 上に円として表現している。円の大きさはマグニチュードと震源の深さに依存している(これは大まかな目安であり、揺れた地域や震度・エネルギーに正確に対応しているわけではない)。また、震源の深さによって円の色が変わり、浅ければ赤に、深ければ緑に近づく。海の色は時刻を現しており、正午が最も明るい青で表現され、深夜であれば真っ黒となる。データについては、 日本気象協会 が提供している3月9日から20日までの 地震情報 を利用させて頂いた。 操作方法は、スペースキーでポーズ、'N'キーで3時間進め、'B'キーで3時間戻る。'+'キーで時間の進みが速くなり、'-'キーで遅くなる。マウスの左クリックでドラッグすれば地図を動かすことができ、右クリックで上下にドラッグすれば拡大・縮小となる。 This browser does not have a Java Plug-in. Get the latest Java Plug-in here. 可視化地図から、3月11日14時46分のマグニチュード9.0の地震の起こる二日ほど前から、三陸沖で中規模の地震が頻発していることが分かる。ところが、M9.0の地震の数時間前からは揺れが全て途絶え、嵐の前の静けさの様相となる。そして、M9.0の地震の直後からはひっきりなしに揺れが続き、数日してやや沈静化してきている。ただ、M6程度の地震は未だに起こっているので、安心するにはまだ早いかもしれない。最近の地震では福島県沖や茨城県沖に多く、ややと南下している様子が分かる。 以下、ソースコードを示す。 earthquake.pde /* earthquake.pde by nox, 2011.3.22 */ String convData ...

昨年は開発エンジニアが転職するという話題が多かったように思う。それに触発されたというわけではないが、自分も昨年11月に開発エンジニアとして転職した。そしてつい先日、試用期間の3ヶ月が過ぎた。これで何かが変わったというわけではないが、正式に採用されたということで、転職して思うことなどをブログに書くことにした。転職先については敢えて名前を出さないが、主にマルチコアやGPGPUなどの並列処理を扱うベンチャー企業といえばピンとくる人もいるかも知れない。 自分はこれまでずっとアカデミックな研究所に所属していたのだけど、少々居座りすぎてしまったことは否めない。基本的に研究者は流動的であることが推奨される。そして、所属していた研究所は所長でさえも任期制で、研究員は普段から職探しをするのが当たり前になっている。しかし、アカデミックなポストは空きも少なく、枠が一人のところに100人200人と応募してくることはザラで、終身雇用の職を得るのはかなり難しい。 そんな中、自分はほとんど職探しをしていなかった。第一の理由として、自分が専攻した薬学の分野でバリバリとコードを書くことと研究が両立できるか疑問だったことがある。ただでさえ狭き枠にそのような特殊な環境を求めるのが非常に厳しいものに思えたからだ。第二の理由として、所属していた研究所ではメインでHPCを扱っており薬学に関連した研究を遂行するためのプログラムを自由に書くことができ、それが求められていたことがある。相手の要求と自分の希望が合致していた。 そうは言っても、ここ最近の政治的判断による不透明さや年齢に伴うマネージメントの必要性など、当たり前だが要求は厳しくなってくる。また、地理的な面でも大きく変わることになったし、こちらの事情で現在の住処を離れたくないということもあった。研究所で一緒に仕事をしていた人たちには大変お世話になっていたので少々迷うところもあったのだが、最終的に転職することにした。 自分の仕事選びの基準は単純で「楽しめるかどうか」だけだ。楽しめればどんなにきつい仕事でも苦にならないし、逆にきつければきついほどチャレンジングであり、やりがいを感じられる。特に開発エンジニアの方なら、自分の技術で困難な問題を克服する達成感、創造する喜び、知識を得る楽しみなどについて分かってもらえるのではないかと思う。それに、企業と...

範囲(-0.005, -0.005, 0.005, 0.005)のマンデルブロ集合を描画サイズ4,800x4,800、繰り返しの上限10,000として計算させてみた。C++でコードを書いてXeon X5650 1コアで走らせたところ942.97秒かかった(SIMD最適化はしていない)。TBBを利用した12論理コアの並列計算では88.68秒だった。次いでCUDAを使ってTeslaで走らせてみたが7.45秒まで短縮された。Teslaパない。そして贅沢にも2枚のTeslaをスレッドで並列化して使ってみたら4.26秒で計算できた。実にXeon 1コアの220倍の速度が出たわけだ。カリカリチューンをしなくてもこの程度の高速化ができるということが重要だ。 因みにマンデルブロ集合の上記の範囲は数値が発散しない集合部分であり、描画させてみても真っ黒なので面白くはない。これは並列化した際に計算が偏らないようにしたかったのでこのような範囲を指定している。プログラムで出力させたファイルは色の生データなのでそのままでは表示できない。一応、画像表示するためのPythonコードも付け加えておく。 以下に、GPU2枚を利用してマンデルブロ集合を計算するCUDAコードを示しておく。 mandelbrot_thread.cu // Mandelbrot set using GPGPU by nox, 2011.02.12 // nvcc -lcutil_x86_64 -arch sm_13 -use_fast_math -prec-sqrt=false -keep -L ~/NVIDIA_GPU_Computing_SDK/C/lib -I ~/NVIDIA_GPU_Computing_SDK/C/common/inc -g mandelbrot_thread.cu -o mandelbrot_thread #include <iostream> #include <fstream> #include <cutil_inline.h> #include <multithreading.h> using namespace std; const int BLOCK_SIZE_X = 16; const int BLOCK_SIZ...