良いもの。悪いもの。

スピログラフに初めて触れたのは小学校に上がる前だった。当時はそれがスピログラフという名称であることを知らず、輪っかをぐるぐる回して綺麗な模様を描ける道具としてお気に入りの玩具になっていた。自分はマンデルブロ集合やジュリア集合のような数学的に表現できる図形や模様にともて魅力を感じるのだが、幼少のころのスピログラフがその切っ掛けだったのかもしれない。 長方形の定規に大きさの異なる3枚の歯車が付属しており、そして2つの円がくり抜かれている。さらに、矢印、五角形、半円などの型が子供心をくすぐったものだ。因みに、歯車と円に付いている歯数は刻印されており、3つの歯車はそれぞれ36、52、63、円の方は96、105の歯を持っている。 幼少のころ以来、スピログラフを見かけることはなかったのだが、最近になって100円ショップでも売られているという噂を聞き、少し探ってみた。近くの100円ショップでは扱っていなかったが、「デザイン定規」という名称でネットで78円で売っていたので思わず購入してしまった(当然だが送料のほうが高くついた…)。 スピログラフは一般的には内トロコイド(hypotrochoid; 内余擺線)といい、以下の式で表現できる(Wikipediaより)。 定円の半径 r c 、動円の半径 r m 、描画点の半径 r d とし、下図では、それぞれ、青、緑、赤の線で示している。そして、回転角を θとした軌跡がスピログラフの模様となる。 スピログラフでは半径よりも、円周の歯の数で示した方が分かりやすいだろう。円周と半径の比は変わらないのだから、歯の数の比がそのまま半径の比となる。例えば定円の歯が30で動円の歯が10ならば、 r c : r m =3:1ということだ。 さて、幼少のころはこんなに綺麗な模様が数式で表現できるなどと夢にも思わなかったわけだが、今はそれを理解し、コンピュータを使って自由に描くことができる。コンピュータであれば、動円が定円の外側に来る外トコロイド(epitrochoid; 外余擺線)も描くことができるし、描画点を動円の外に出すことも思いのままだ。 せっかくコンピュータを使って簡単に表示できるようになったというのに、プログラムを組むことが難しくては意味がない。しかし、Python+matplotlibであれば簡単に...

世の中のことをもっと知るにはどうしたら良いだろうと思うときがある。世の中の多くの事柄はログやデータに落とされる。Googleなどの検索サイトは良い例だろう。さて、そのログやデータをどうすれば良いのか? 多くの場合、視覚化が有効な手段となる。 まずは身の回りの日常的なデータやログを何とかしたい。ただ、日常のデータを視覚化するのに数十行以上のコードは書きたくない。まるで息をするかのごとく自然に視覚化を行いたいのだ。そのためには1～2行、長くて数行で済ませることが必要だ。そこでPythonとmatplotlibを使う。加えて、IPythonがあればなお良い。IPythonの導入については以前のブログ記事である IPythonの埋め込みプロットが素晴らしい を参考にして欲しい。 まずは事前にnumpyとmatplotlibをインポートしておく。できればscipyも。 >>> from numpy import * >>> from pylab import * 短いコードで視覚化を行うためには、Pythonの内包表記は必須だ。例えば、5, 2, 1, 5, 8をデータとするグラフを書きたいのならIPythonを使って以下の1行で実現できる。 >>> plot([5, 2, 1, 5, 8]) 数値が行ごとにinput.txtというファイルに書かれていた場合は以下の通り。 >>> plot([int(x) for x in file("input.txt")]) map関数を使えばもっとスマートに書ける。ファイル内の数値を文字列として配列で読み込み、それらの文字列をmap関数によりintで整数に変換する。 >>> plot(map(int, file("input.txt"))) さて、ログが2次元で書かれていた場合はどうするのか。例えば以下のデータがinput.txtに書かれていたとする。 2, 2.5 5, 6.2 6, 3.6 7, 6.3 10, 1.9 この場合、以下のようにすればプロットできる。ファイルから行ごとの文字列を読み出して、それを","を区切りとしてリス...

あけましておめでとうございます。 新年を迎えて、ふと、今年の西暦である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モジュールがあることを知ったので、早速インストールして使ってみた。 ところで、大きな素数といえばメルセンヌ素数が有名だが、 マラン・メルセンヌ は、2 n -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),...

先日、「 『フカシギの数え方』 おねえさんといっしょ！ みんなで数えてみよう！ 」という動画を見た。格子状のマスの左上から右下までの経路が何通りあるのかを調べて、格子が多くなればなるほど組み合わせの数が爆発的に増えることを教えてくれる動画だ。これは 自己回避歩行(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スクリプトで統合している。これは、クヌー...

5年ほど前、このブログで 虚数のテトレーション という記事を書いたことがある。 テトレーション(tetration) とは、自らのべき乗を指定された回数反復する演算のことで、 n a と表現する。 3 5 の場合、5 5 5 = 1.911×10 2185 となる。Pythonの関数で表現すれば以下のようになる。 以前の記事では、 ∞ i が 0.43828+0.36059 i に収束することを見つけたのだが、今回はそれに関連したフラクタルについて紹介したい。 テトレーション n a の 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++による実装: % ...

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

先日、地元のまつりがあった。そこの会場となっている区民館では子供用ボーリングやプラ板遊びなどの催しが開かれていたのだが、その中の一つに「ひとりぼっち」というゲームを行なっているのを見つけたので、一緒にいた小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座標 隣接ノード番号リスト"としてノード数の...

今回は興味を持ってもらえる人が限定される記事だと思うけど、以前作ったコードを腐らせてしまうのも勿体無いので敢えて晒してみる。 生体内のタンパク質構造を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 ...

問題 : 1, 5, 10, 50, 100, 500グラムの6種類の錘があります。それぞれの錘はいくつでも使うことができます。ある重さを量るときに使う錘の組み合わせは何通りになりますか。それを求めるプログラムを書いてください。ここで量る重さはグラム単位の整数とし、最大で100キログラムとします。10グラムの重さを量る場合、10グラムの錘が1個、5グラムの錘が2個、5グラムの錘1個と1グラムの錘が5個、1グラムの錘が10個の組み合わせになり、全部で4通りとなります。 以下、解答例。 weights.cpp // 指定した重さを量るのに使う錘の組み合わせの数を求めるプログラム. // 錘は 1g, 5g, 10g, 50g, 100g, 500g の6種類. // 重さ 700,000g 程度までなら桁あふれせずに計算できる. #include <iostream> #include <map> #include <algorithm> #include <cstdlib> using namespace std; const int weights[] = { 1, 5, 10, 50, 100, 500 }; // 錘(昇順). const int N = sizeof(weights) / sizeof(int); // 錘の種類. // 求めた錘の組み合わせの数のメモ. // 下位3ビットは使用する錘(solveのidx)を表す. // 残りのビットは重さ(solveのv)を表す. map<int, unsigned long long> memo; // 指定した重さに対する錘の組み合わせの数を求める. // v : 重さ // idx : 使用する錘 // N-1の場合はweights[0]からweights[N-1]までの錘を使用する(すべての錘). // 2の場合はweights[0]からweights[2]までの錘を使用する(1, 5, 10の錘). // 戻り値: 組み合わせの数 unsigned long long solve(int v, int idx=N-1) { // 既に...

Haskellによる エラトステネスの篩(sieve of Eratosthenes) の美しい実装を見て、初めてアルゴリズムの素晴らしさに気付かせてくれたのがエラトステネスの篩だったことを思い出した。たしか中学生の頃だ。そこで、当時を懐かしみながら簡単な実装を書き留めておくことにした。とりあえず、Haskell, C++, Pythonの実装を以下に示す。コードは比較的短いが、実行効率を優先させているわけではない。 Haskell これはHaskellのサンプルコードでよく出てくる無限リストと遅延評価による実装だが、とても分かりやすいし、美しいコードだと思う。 primes = sieve [2..] sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p /= 0] 以下のように関数 take を使って必要な分だけ素数を取り出すことができる。ただし、実行効率はあまり良くない。 take 10 primes [2,3,5,7,11,13,17,19,23,29] C++ 次にC++を使って実装してみた。ここではC++0xを使っている。そのほうが簡潔で分かりやすく書けるし、これからはC++0xがより使われ、普及して欲しいという意味もある。因みにgccではバージョン4.5以降、Intelコンパイラではバージョン11.0以降、Visual C++では2010(16.0)以降でコンパイルできる。 std::vector<int> primes; for (int i = 3; i < 100 ; i += 2) primes.push_back(i); auto end = primes.end(); for (auto x = primes.begin(); *x * *x <= *(end-1); ++x) end = std::remove_if(x + 1, end, [&x](int p){ return p % *x == 0; }); primes.erase(end, primes.end()); for (auto p = primes.begin(); p != primes.end(); ++p) std::cout <...

GIGAZINE に 数学のエキスパートが3ヶ月かけて作成した「世界一難しい数独」 なるものがあったので、即席でPythonワンライナー(一行プログラム)を作ったら1秒かからず解けた。問題を作るのは大変でも解析プログラムで一瞬とは儚い。 python -c "import sys;L=[];S=lambda D:(0in D)and[L.append(D.index(0)),[(D.__setitem__(L[-1],i),S(D),D.__setitem__(L.pop(),0))for i in set(range(1,10))-set(D[L[-1]/9*9:L[-1]/9*9+9]+D[L[-1]%9:81:9]+[d for n in(0,1,2)for d in D[L[-1]/27*27+L[-1]%9/3*3+n*9:L[-1]/27*27+L[-1]%9/3*3+n*9+3]])]]or([sys.stdout.write('%d'%d+('\n'if i%9==8 else' '))for i,d in enumerate(D)],sys.exit());S([int(c)if c!='.'else 0for c in'..53.....8......2..7..1.5..4....53...1..7...6..32...8..6.5....9..4....3......97..'])" コードの最後のリスト内の文字列が9行9列の数独の問題となっていて、上の行の数字から順に入っている。未知の部分はドット(.)で表している。上記のコードではGIGAZINEの問題となっている。 因みにPythonによるもっと効率の良いアルゴリズムについては Solving Every Sudoku Puzzle が参考になるかも。 以下、ワンライナーによる出力結果。 1 4 5 3 2 7 6 9 8 8 3 9 6 5 4 1 2 7 6 7 2 9 1 8 5 4 3 4 9 6 1 8 5 3 7 2 2 1 8 4 7 3 9 5 6 7 5 3 2 9 6 4 8 1 3 6 7 5 4 2 8 1 9 9 8 4 7 6 1 2...

小学3年生の娘の日曜授業参観に行ってきた。算数の授業だ。授業の後半、以下のような問題が出された。 問題 : ゴマダラチョウとトノサマバッタがあるゲームをしている。0から9までの数字の書かれたカードがそれぞれ一枚ずつ全部で10枚あって、それを使って3桁の数を2つ作り、その差をできるだけ小さくした方が勝ちとなるゲームだ。ただし百の位は0にできない。できるだけ小さくするにはどのようにカードを選べばよいだろうか。 これを約30人の生徒に考えさせていた。解き方を先に教えるということはしない。生徒が問題を考えている間、先生は生徒たちを見回り、質問などに答える。しばらくするといろいろと答えが挙がってきた。自分の娘は以下のように考えたようだ。 百の位は差が1であればどの数字でも良いので、まず十の位を最小にする数字を考えてみると、最小の数字0から最大の数字9を引いた場合が最も小さくなる。同様に一の位では0と9以外の最小・最大の数字を選ぶ。つまり1から8を引いた場合が最小になる。百の位は残りの数字カードから差が1となるものを選ぶ。解答例は「501-498」や「701-698」となり、差は3になる。 これはまさにアルゴリズムだ。Pythonであれば以下のコードと同じだろう。 digits = range(10) ten = (digits.pop(), digits.pop(0)) one = (digits.pop(), digits.pop(0)) idx = random.randrange(len(digits) - 1) hund = (digits.pop(idx), digits.pop(idx)) a = reduce(lambda x, y: x * 10 + y, map(array, (hund, ten, one))) print "%d - %d = %d" % (a[1], a[0], a[1] - a[0]) 出力例: 601 - 598 = 3 小学3年生の算数だと思ってそれほど期待せずに観に行ったのだが、生徒たちはみな楽しそうに問題に取り組んでいたし、なかなか面白い授業をしているようで何だか安心した。因みに、今回のようなクラス全員で考える問題は普段からよく行っているとのことだった。

先日、 構造体を格納したコンテナに対してSTLのアルゴリズムを適用する方法 について記したが、今回はそれを C++0x で書き直してみた。C++0xの規格はまだきちんと定まっていないが、草案ではかなりのところまでつめられてきており、また対応するC++コンパイラも出てきていることから今のうちにC++0xに慣れておいた方が良いと思う。今回のプログラムで行う内容については前回と全く同じだが、C++0xを使うことによりかなりすっきり書くことができた。因みに現時点の最新の規格(草案)は N3000 (PDF)になる。日本語での資料は C++0xの言語拡張まとめ (Faith and Brave - C++で遊ぼう) や 本の虫 が詳しい。余談だが、2010年になっても0xから1xに変わることはないようだ。C++の生みの親であるBjarne Stroustrupによれば、 0xのxは16進数と考えて欲しい とのこと。 C++0xでは多くの便利な機能が追加されているが、今回はラムダ式、auto型指定子、初期化子リストを使用した。使用したコンパイラは Visual Studio 2010 beta2 (VS2010) 、 GCC 4.5 (snapshot 20100114) 、 Intel Compiler 11.1 だが、残念ながらVS2010とIntel Compiler 11.1では初期化子リストの機能が実装されていないので、それらを使用する場合、今回示したコードのうち素数を格納する静的メンバコンテナの初期化部分は変更しなくてはならない。また、GCCについては最新リリース版である4.4.2で ラムダ式が使えない ので、ここでは開発版の4.5を使用している。その他のコンパイラも含めてC++0xへの対応は C++0xCompilerSupport で確認できる。 初期化子リストを使用しない: static int init_primes[] = { 2, 3, 5, 7 }; vector<int> IsPrime::primes(init_primes, init_primes + 4); 初期化子リストを使用する: vector<int> IsPrime::primes = { 2, 3, 5, 7 }; 今回の...

整数の分割(integer partitions) を列挙するプログラムを書いてみた。ここでは再帰を用いた比較的単純な実装を示すが、より効率的なアルゴリズムを考えるとなると奥が深い。 整数の分割とは、ある整数を自然数の和で表したものだ。例えば分割する整数を5とすると求める自然数は以下のようになる。 5 4 1 3 2 3 1 1 2 2 1 2 1 1 1 1 1 1 1 1 今回のC++で書いたコードは単純で理解しやすいと思うが、高速化などは考えていない。もしよりよいアルゴリズムに興味があるのなら Antoine Zoghbiu and Ivan Stojmenovic. Fast Algorithms for Generating Integer Partitions. Intern. J. Computer Math. , 70 , 319-332 (1998) (PDF)あたりが参考になるかも。 使い方は以下の通り。 partitions 分割したい整数 [ 分割した要素の最大値 ] 分割したい整数が5であれば、 partitions 5 とする。出力は次のようになる。 5 (max 5): 7 5 4 1 3 2 3 1 1 2 2 1 2 1 1 1 1 1 1 1 1 また、分割したい整数が5で、分割した要素の最大値を3とする場合は、 partitions 5 3 で、以下のような出力を得る。 5 (max 3): 5 3 2 3 1 1 2 2 1 2 1 1 1 1 1 1 1 1 partitions.cpp #include <iostream> #include <vector> #include <algorithm> #include <cstdlib> using namespace std; void partitions(int n, int max_n, vector<int>& v, vector<vector<int> >& vv) { if (n == 0) { vv.push_back(v); return; ...

「 騎士の巡歴(Knight's Tour) 」を解くプログラムをC++で書いてみた。騎士の巡歴とはチェスボード上の駒「ナイト」を移動させてすべてのマスを一回ずつ通過させる問題だ。特に最初の駒の位置に戻ってくる解を「騎士の周遊(Closed Knight's Tour)」と呼ぶ。 一応、C++らしいコードを心掛けてみたけど、思ったより長くなってしまった。騎士の巡歴では Warnsdorffのアルゴリズム 、騎士の周遊においてはそれに加えて Schwenkの定理 を用いている。また、 Puzzle DE Programming で書かれている方法も参考にさせてもらった。 まず、ボード上のマスをNode構造体のポインタで表現し、 *node とする。その node からナイトが動けるマスを vector<Node*> next; として保持する。また、一度通った node は node->visited にそれを記録しておく。 next を順に巡っていくが、Warndorffのアルゴリズムにより、移動先となる next 内の node のうち、動けるマスの最も少ない node を移動先とする。すべてのマスを通れば探索終了だ。 また、騎士の周遊については、最初にSchwenkの定理から解があるかどうかを判断する。解があれば、騎士の巡歴と同様の手順でまず解を見つける。ここで、見つけた解の最初のマスからナイトが動けるマスをA、最後のマスからナイトが動けるマスをBとし、BからAの順となるマスがあるかどうか調べる。もしあれば、そのAのマスから最後のマスまでの順路を逆にすればそれが騎士の周遊の解となる。もし、BからAの順となるマスがなければ、動けるマスが見つかるまでバックトラックにより探索する。この手順については上記の Puzzle DE Programming が詳しい。 プログラムの使い方は以下の通り。 knight 行数 列数 [駒のスタート行=0 駒のスタート列=0 巡歴か周遊か?] スタート位置はデフォルトで(0, 0)で、何も指定がなければ騎士の巡歴を求める。騎士の周遊の場合は、コマンドラインの最後に c を付ける。 例えば、16x16マス上で騎士の周遊を解く場合は以下のようにする。 knight 16...