C++の小ネタ。マルチスレッド/マルチプロセスで出力が溢れかえるのをなんとかしたかったので作りました。

print文デバッグは実行中の変数の状態を知るのに手軽な方法としてしばしば使いますが、マルチスレッド/マルチプロセスで動く関数をデバッグしていると、スレッド/プロセスの数だけ一気に出力されてしまって大変見づらくなります。

マルチスレッド化の最中など全スレッドの状態を見たい場合もあるのですが、そうでないときはいろいろと混ぜながら出力されてしまうので何がなんだかわからなくなります。

こんなときはif文を使って一応の解決ができます。どこかに出力フラグisOutのようなものがあったとして、

#include

void someThread(bool isOut){
  if(isOut) std::cout << "hogehoge" << std::endl;
}

みたいに書けば一応解決します。しかしいちいちif文まで書くのはどうも面倒だし、そもそも、C++のストリーム関数はいくらでも好きなものを定義できるのが利点だったはずです。

そんなときのために出力するかどうかを選べるようになっていると便利ですね。

というわけで作りました。

続きを読む →

小ネタ。プログラミングの合間にふと思ったちょっぴり邪悪なアイデアを実装してみました。

C++のbool型変数で、そのままでも否定してもtrueになる変数を用意できるよ、というお話です。
組み込み型に未定義の値を入れたらどうなるんだろうと思ってやってみました。

ゆるC++プログラマーを混乱させることができます。良い子のみんなは真似しないように。

おまけで書きますが、良い子のclangには効きませんでした。

続きを読む →

全体は部分の総和に勝る。 – アリストテレス

マンデルブロー集合、10の200乗倍、グーゴル倍のグーゴル倍に挑戦しました。

Youtubeのほう、最後がぶつ切りになってしまっていますが許してください。
ニコニコ動画はこちら。

(追記)もっと画質のよい動画が欲しい場合はGoogle Driveからダウンロードして下さい。1GB弱あります。

現実世界だと陽子や電子の半径が10の-15乗メートル、観測可能な宇宙の範囲でもたったの(?)10の27乗メートルしかないため、これほどの倍率はまず見かけることはありません。

ニコニコ動画で何度も見返した動画というのはそれほど多くないのですが、その一つに「【超高倍率】マンデルブロー集合【10の100乗倍】 (hoho様)」というものがあります。今回はその動画に感謝しつつ、記録更新を狙いました。
先の動画ではニコニコ最高記録と記されているので、それが正しければ今回がニコニコ新記録です。(logスケールで)一気に2倍まで来たことになります。

続編はこちら:  マンデルブロー集合へ飛び込む

続きを読む →

「サンクトペテルブルグのパラドックス」は、現実的にはあまり儲からなそうな賭けの期待値が、計算上無限大に発散するとされる有名なパラドックスです。

パラドックスを知らない人のために賭けの要点を書いておくと、以下のようになります。

「コインを投げ続け、表が出たらゲーム終了。表が出たのが1回目なら1円獲得、2回目なら2円、3回目なら4円…というように、裏が1回出るたびにもらえる金額が倍になる。いくらの賭け金ならこの賭けに参加するべきか。」

儲け額の期待値Wは以下のように計算されます。1/2で1円、1/4で2円、となるので…

ここまで見てだまされないぞと思った人は、しっかり悩んでみてください。

先日友人との間で話題になって、実際のところおおよそどんな具合で儲かるのかやってみようということでさくっと書いてみたものです。確率についてなんやかんやと議論するのも楽しいですが、ゴリ押しでデータを集めるのもまた面白いものがあります。あと探しても意外とそういう資料が見当たらなかったので。

---

ソースコード

いつもどおりC++です。C++11の機能を使っていますので、例えばgccとかでコンパイルしたい場合は -std=++11 とか -std=c++0x とかのオプションをつけてください。

#include <ctime>
#include <iostream>
#include <random>
#include <functional>

int main(int argc, char** argv){
  if(argc != 2){
    std::cerr << "Wrong argument!" << std::endl;
    return -1;
  }
  long long int loopNum = 1;
  for(int i = 0; i < atoi(argv[1]); i++){
    loopNum *= 2;
  }
  std::cout << "Loop: " << loopNum << std::endl;

  time_t timer;
  time(&timer);
  auto rand = std::bind(std::uniform_int_distribution(0,1),
      std::mt19937(static_cast<unsigned int>(timer)));
  long long int total = 0;
  for(long long int trial = 0; trial < loopNum; trial++){
    long long int earn = 1;
    while(rand()){
      earn*= 2;
    }
    total += earn;
  }
  std::cout << "Earned: "
    << static_cast<double>(total)/static_cast<double>(loopNum)
    << std::endl;
  return 0;
}

特記するようなものもありません。コマンドライン引数から2の何乗回計算するかを決めて、律儀に試行を繰り返して平均獲得金額を求めています。

カウントがlong long int型なので2の50乗回程度までは問題ないと思いますが、そのはるか手前の2の30乗前後(数億程度)で計算時間がきつくなってくると思います。
まあ乱数も毎回1ビットしか使っていないしシングルスレッドだしで、高速化する余地はまだ残っています。

---

結果

大量に試行した結果は以下のようになります。あくまで愚直に。力強く。

横が1セットあたりの試行回数、縦が1回あたりの平均獲得金額です。両方logスケールなのに注意。
(表が出るまでを1回の賭けとし、1セットでn回行うと定義しています。平均獲得金額は儲け額をnで割った値。)

不思議なグラフです。まず、回数を増やすと平均獲得金額は順調に上がっていくことがわかります。1,000,000,000回以上やるつもりなら、1回あたり10円以上儲かることがほぼ確実になります。このままグラフを右に伸ばしていったらどうなるでしょう。

もうひとつ面白いことがあって、それはばらつきが減らないことです。例えば右上にある点、だいたい10,000,000回で平均獲得金額が10,000円を超えていますが、つまりこの1セットで単純にトータル一千億円以上儲けているということです。とんでもない大当たりです。

胴元と挑戦者、どちらの席につきたいですか?

---

今回は横軸がlogスケールで均等になるようにサンプリングしています。離散的な値しかとらないので左側はすかすかに見えますが、実際は点が重なっているだけです。
また、左上から右下へのぼんやりとした斜線が大量に見られます。これは「1回大当たり」を成し遂げたときの線になっています。回数を増やすと同じ大当たりでも「1回あたりの獲得金額」は減っていくので右下がりになるという具合です。その代わり「もっとすごい大当たり」をひく確率も出てきます。

モンテカルロ法などでは試行回数を増やすと値が落ち着くのですが、この問題はおもしろいことにいくら試行回数を増やしてもばらつきが収まらないようになっています。よくできています。

儲け額の分散も出してみようと一度書きかけたのですが、ちょっと考えてからこの問題に対して無駄な試みであることがわかりました。

コンピュータはしばしば想像力を働かせる道具として役に立つのですが、今回のように神通力の効かないヤバそうな問題ではちょっとわけが違ってきます。こんな問題でも扱えてしまう数学の力は偉大です。

結果の意味とこの賭けの本当の期待値についてはいろいろ考えてみて下さい。

前半の補足です。前半で紹介したシミュレータ云々について。

---

この粒子の振る舞いについて

「なんで星の運動を模してるのに楕円軌道じゃないんだ?」と思った方はとても勘のいい方です。

これは2次元の運動をどう考えるかによるものです。
3次元空間の平面上に粒子を配置させて運動させた場合、重力は距離の2乗に反比例しますが、
3次元空間に無限長の棒を(平行に)置いて運動させた場合、重力は距離の1乗に反比例します。

これはまあ真面目に積分してもそうなりますし、電磁気学を勉強したことのある方は雰囲気で理解できるかと思います。後者は宇宙空間では仮定が不自然すぎるだろうと思われるかもしれませんが、万有引力が3次元空間で逆2乗の法則に従っているとするなら、2次元では逆1乗になるのが自然だとも考えられます。

まあ、今回はそこまで哲学的な話ではなく、単に逆1乗のほうが見栄えが良かったからです。

どう見栄えに影響するかですが、これは無限遠点での挙動の違いを考えればだいたいわかります。位置エネルギーは力の積分(の-1倍)で与えられますが、これより力が逆2乗、逆1乗の場合のそれぞれについて位置エネルギーは逆1乗、logとなります。つまり前者ならエネルギーが無限遠である値に収束するのに対して、逆1乗なら際限なく増えていきます。

つまり、逆2乗ではひとたび粒子が第三宇宙速度に達してしまうと二度と戻ってこないのに対し、逆1乗なら少なくともいずれは引き戻されてくるという大きな違いがあります。見た目に粒子がどこかに飛んでいってしまうのは寂しいですね。

よく「この世が2次元だったらいいのに」という言葉を聞きますが、そんなことをしたら、上記の理由から世界中のものがある一点に集まってしまってつまらない世界になるだろうなあと思っています。

なお、粒子があまり近づきすぎると計算上の問題からヤバいことが起きるので、引力を打ち消すだけの反発力を与えてあります。
現実だと衝突しますが、無理にエネルギーを保存させようとして星が完全弾性衝突をするのもあまりリアリティがないですね。

---

以前の3重振り子を使わなかった理由

これは速度ベルレ法の適用範囲に限界があるからです。(下は剛体3重振り子の動画)

シンプレクティック法は、影のハミルトニアンと呼ばれる量を厳密に保存する運動を記述します。この際、位置を計算する際にのときのハミルトニアン(の偏微分)が必要になります。これは普通に考えると右辺と左辺がからまってさくっと解けません。
しかし都合のいいことに、単純な場合ハミルトニアンHは位置qと運動量pの足し算に分離でき、このときハミルトニアンの偏微分の形からqがなくなります。その場合の解法のひとつが速度ベルレ法となるわけです。

以前の3重振り子では、ハミルトニアン(というよりラグランジアン)の計算に3つのθという一般化座標を使っていました。このためハミルトニアンの分離が行えなくなり、(実際にはできるかもしれないけど私はその数学を持ち合わせていないので、)速度ベルレ法をはじめとした陽解法が使えない状態に陥っていました。

こういった場合には陰解法のガウス・ルジャンドル法などを使って解くことができます。

なんとなく当然のような気もしますが、よくよく考えてみるとこの仮定は結構不思議なことです。分子動力学法などではたいていこの仮定が成り立ってしまいます。

---

ソースコード

ちょっとだけかいつまんで書いておきます。C++です。

まずは系の状態の定義。数値積分に関しては位置と速度のペアのリストさえあればいいのだけど、肝心な力の計算ではどの値がどの物理量に対応しているかを知らなければならず、もどかしい。

かといって値を参照して保持しておくほどの大げさなものでもない。そこで、どうせ数値積分では四則演算くらいしか使わないので、粒子のvectorに必要な演算子のみ定義した構造体を持っておきます。

	struct partState{
		double x, y;
		double vx, vy;
		partState()	: x(0.0), y(0.0), vx(0.0), vy(0.0){}
		partState(double _x, double _y, double _vx, double _vy)
			: x(_x), y(_y), vx(_vx), vy(_vy)
		{}
	};
	struct partSystem{
		std::vector<partState> pss;
		partSystem(int size){
			pss = std::vector<partState>(size);
		}
	};
	const partState operator+(const partState& lhs, const partState& rhs);
	const partState operator-(const partState& lhs, const partState& rhs);
	const partSystem operator+(const partSystem& lhs, const partSystem& rhs);
	const partSystem operator-(const partSystem& lhs, const partSystem& rhs);
	const partState operator*(const double lhs, const partState& rhs);
	const partSystem operator*(const double lhs, const partSystem& rhs);

力の計算はもうほんとに全ペアで単純に逆1乗で計算してるだけで、高速化のかけらもありません。
ガチで分子動力学法などをやるときには今回のO(n^2)のようにならない工夫が必要です。

		partSystem calcForce(const partSystem& _ps){
			partSystem ans(size);
			for(int i = 0; i < size; i++){
				for(int j = i+1; j < size; j++){
					double dx = _ps.pss.at(i).x - _ps.pss.at(j).x, dy = _ps.pss.at(i).y - _ps.pss.at(j).y;
					double r2 = dx*dx + dy*dy;
					double r = std::sqrt(r2);
					if(r > r_threshold){
						//double r12_e3 = r2/0.0001;
						//r12_e3 = r12_e3*r12_e3*r12_e3*r12_e3*r12_e3*r12_e3;
						double f = -g*masses.at(i)*masses.at(j)/r;// + g2*masses.at(i)*masses.at(j)/r12_e3;
						double fx = f*dx/r, fy = f*dy/r;
						ans.pss.at(i).vx += fx;
						ans.pss.at(i).vy += fy;
						ans.pss.at(j).vx -= fx;
						ans.pss.at(j).vy -= fy;
					}
				}
				ans.pss.at(i).x = _ps.pss.at(i).vx;
				ans.pss.at(i).y = _ps.pss.at(i).vy;
			}
			return ans;
		}

反発力を与える際、分子動力学チックに高次の多項式の反発力を与えようとした名残があります。結局そこまでする必要もなく削除。かわりに距離にしきい値を導入します。
エネルギーの計算(表示用)でも不整合がないようにしておく必要があります。

あとは今回のテーマだった数値積分。ソースはとっても単純。

		void move_euler(const double dt){
			auto psDiff = calcAccel(ps);
			ps = ps + dt * psDiff;
		}
		void move_heun(const double dt){
			auto k1 = calcAccel(ps);
			auto k2 = calcAccel(ps + dt*k1);
			ps = ps + dt/2.0*(k1+k2);
		}
		void move_runge_kutta(const double dt){
			auto k1 = calcAccel(ps);
			auto k2 = calcAccel(ps + dt/2.0*k1);
			auto k3 = calcAccel(ps + dt/2.0*k2);
			auto k4 = calcAccel(ps + dt*k3);
			ps = ps + dt/6.0*(k1+2*k2+2*k3+k4);
		}
		void move_velocity_velret(const double dt){
			auto k1 = calcAccel(ps);
			auto p2 = ps;
			for(int i = 0; i < size; i++){
				p2.pss.at(i).x += dt*ps.pss.at(i).vx + dt*dt/2.0*k1.pss.at(i).vx;
				p2.pss.at(i).y += dt*ps.pss.at(i).vy + dt*dt/2.0*k1.pss.at(i).vy;
			}
			auto k2 = calcAccel(p2);
			for(int i = 0; i < size; i++){
				p2.pss.at(i).vx += dt/2.0*(k2.pss.at(i).vx + k1.pss.at(i).vx);
				p2.pss.at(i).vy += dt/2.0*(k2.pss.at(i).vy + k1.pss.at(i).vy);
			}
			ps = p2;
		}

速度ベルレ法だけ、位置と速度を明示的に扱っているのがわかります。

前々回の振り子の欄でちょっと触れた運動方程式の数値計算の話の続きです。

運動方程式は立てるだけでなく、実際に方程式を解かないとシミュレーションは行えません。

一番丁寧な運動方程式の解き方は解析的に解くことですが、ほとんどの場合うまい方法がないことがわかっています。例えばN体に万有引力が働いて運動する系ではNが3以上になると解析的に解けなくなることが知られています。太陽と地球と月でもうオーバーです。

もっと一般的に扱いたいときは、コンピュータが強力な武器になります。微分方程式を差分方程式に近似してしまえば計算機で簡単に扱えるので、後は時間の刻み幅を十分に小さくとってあげれば微分方程式に近い解が得られるというわけです。

じゃあそれで話は終わりかというとそうでもなく、この差分方程式にする段階がなかなかうまくいかないわけです。それに伴っていろいろな手法が提案されています。今回はその話です。

---

今回作ったものは互いに引力を及ぼし合う粒子たちのシミュレーションです。

イメージしているのは銀河系的なアレですが、まあ数値積分の例で作っただけなので本気で似せているわけではありません。私は中学校で習う電子の描像のイメージに似ているなと思いました。
IAEAのマークにもあるアレです。
中央に質量の大きめな粒子をひとつ、あとは小さな粒子をばら撒いています。小さな粒子同士も相互作用します。

大事なのはこの系では力学的エネルギー・運動量・角運動量が保存するということです。シミュレーションでエネルギーが保存していなかったら、それは何かおかしいことが起こっているぞ、というわけです。

可視化の部分は振り子の流用ですね。数字キーでシミュレーション手法を変えたり、初期化したり軌跡のオンオフを行ったりします。

後述するオイラー法、ホイン法、ルンゲ・クッタ法、速度ベルレ法を実装しています。

いつものようにプログラム: particle
ソースコード: GitHub

---

微分方程式を解く際、おそらく最も単純な差分化はオイラー法です。

ニュートンの運動方程式は2階の微分方程式ですが、速度v自体も変数と置いてしまうことで1階の微分方程式に置き換えることができます。以下の式で左辺をx、右辺をfと置いてしまえばいいわけですね。

オイラー法では以下のような差分化を行います。⊿tは十分に短い時間をとります。

見た目にわかりやすいため学習用として頻繁に登場しますが、正確なシミュレーションをしたい場合向きではありません。誤差がとても大きいからです。

誤差が大きい理由は2点の補間を直線、つまり1次式近似にしているからで、それなら2次や3次の項を追加して補間しようというアイデアがあります。実際これはうまく行きます。
よく使われるものに4次多項式で補間するルンゲ・クッタ法というのがあります。

tからt+⊿tまでの4点で傾きを計算して、適当に混ぜることで4次の精度を出しています。式中の諸係数を求めるのにはそこそこの手間がかかります。
次数を上げるほど多項式補間の精度は上がりますが、そうすると計算の安定性の議論から時間刻み⊿tをかなり小さくとらないといけなくなってしまうこと、また計算コストが上がるのでそれなら次数を上げずに時間刻みを減らしたほうが都合がいい場合が多いことなどがあり、結果としてこの4次のルンゲ・クッタ法はよく使われます。

1次、2次、4次の各手法で時間経過とエネルギーのずれをプロットしたものが以下になります。

1次のオイラー法は目に見えて粒子が遠ざかっていってしまいます。2次(ホイン法)と4次(ルンゲ・クッタ法)はシミュレーション結果を見た雰囲気では似ているものの、エネルギーを見るとホイン法はだんだんとずれていくことがわかります。
もちろんオイラー法にも利点はあって、軽くて実装が楽なので、正確さにあまりこだわりがない場合にはいちばん手軽な解法になります。

陰的解法と言われる手法群もあるのですが、一般に複雑なので今回は省略します。

---

速度ベルレ法

ここで、今回メインで書きたかった手法「速度ベルレ法」を紹介したいと思います。2次の精度をもつ手法です。
2次ならルンゲ・クッタ法より精度が悪いのかというとそうでもありません。先のグラフの倍率を上げたものに速度ベルレ法を追加してみます。

ルンゲ・クッタ法に対してぶれは大きいですが、ホイン法のようにエネルギーが変化していくことはありません。しかしまだこれではルンゲ・クッタ法のほうが頼りがいがあるように見えます。
面白いのは長時間の振る舞いです。このまま横軸のスケールを100倍にしてみます。

ルンゲ・クッタ法もこの長時間では計算誤差が積もってエネルギーがずれていくのに対して、速度ベルレ法はほぼ初期の振る舞いを保ったままとなっています。短期間のぶれは大きいのに、長期間でずれていかない不思議なグラフになっています。

速度ベルレ法の式は以下のようになります。位置と速度を分けて扱うのが今までの手法と違います。

なぜこのような振る舞いを見せるのかにはわけがあって、この速度ベルレ法はシンプレクティック性というとても重要な性質を持つように作られています。

ここからちょっと難しい話ですが、ニュートンの運動方程式は時代を追ってラグランジュ形式、ハミルトン形式へと再公式化されていきます。このハミルトン形式で表現したとき、運動は位相空間上で体積を保つ性質があることがわかっています。
速度ベルレ法も含まれるシンプレクティック積分法と呼ばれる手法では、この位相空間での体積を保存する性質を有していて、結果的に長時間でのシミュレーションで(ある種の)誤差が積もっていく現象が起きなくなっています。

ハミルトン形式におけるとある保存量を厳密に保存しており、時間反転対称性(終了状態から反対向きにシミュレーションを行うと開始状態に行くことができる)、エネルギーなどの保存量を保存する、といった今までの手法にはなかった性質があります。

---

速度ベルレ法は、私が研究でよく扱う分子動力学法でもよく使われています。系を長時間観察するときなどにこの性質がとても有効だからです。また分子動力学法では計算コストのネックが力の計算になっており、1回に多段で力の計算を行うルンゲ・クッタ法などは相性があまりよくありません。

式自体は単純なため「分子動力学法で使われる数値計算は単純」という表現をよく耳にします。しかしその構成はよく考えられていて、裏側には先人の知恵がたっぷり詰まっていることは確かです。

---

運動方程式ひとつをとっても様々な数値計算法がありますが、どの積分法が正しいかといった議論はあまり有益ではなく、私のイメージでは「どれも間違っているが、間違え方が違うので場合に応じてどうしても譲れないものに応じて選ぶのが良い」としか言えないかなと思っています。

例えばゲームの物理演算などでは、そこらへんの物が綺麗な放物線を描かなくてもそんなに気になりませんが、プルプル震えたり突然とんでもない勢いで吹っ飛んでいったりすると見た目からしてヤバい、ということになるわけです。それなら計算精度より安定性と計算コストを重視しようとなるわけです。

では、現実のこの世界は、どんな手法で運動方程式が解かれているのでしょう。これはちょっと危険な問いで、そもそも現実世界が運動方程式に厳密に従っているのか、宇宙は連続体なのか、シミュレーション可能なのか、といったことからあれこれ考えなければいけません。
まあでも、私が宇宙をひとつ作っていいと言われたら、たぶんオイラー法を採用します。一番単純だし、世界がそんなに複雑な式でできているとは考えたくないからです。

---

後半では、今回のシミュレーションで作った系について(なぜ前回の振り子のままではだめだったのか、星の運動としてこの動きの胡散臭さは何か、ソースコード云々)の説明をしたいと思います。
速度ベルレ法と言っても万能ではないのです。

Linuxの話です。

続編(?)はこちら:  [フラクタル] マンデルブロー集合 10の200乗倍への挑戦

先日、注目のコンパイラclang/LLVMがC++11に完全対応したというニュースがあったので、せっかくの機会なのでどんな具合なのか試してみました。

円周率ベンチはだいぶ下です。結論を言ってしまうと、速度はそこまで大きく変わりませんでした。<

まずはclang/LLVMのインストール。更新も早いし、基本的には本家ページのガイドラインに従うのがよいと思います。C++11に対応したのはつい最近なので、ソースコードをリポジトリから入手します。

% mkdir llvm
% svn co http://llvm.org/svn/llvm-project/llvm/trunk llvm
% cd llvm/tools
% svn co http://llvm.org/svn/llvm-project/cfe/trunk clang
% cd ../projects
% svn co http://llvm.org/svn/llvm-project/compiler-rt/trunk compiler-rt

コンパイルにはコンパイラやgmakeなどいくつか基本的なソフトが必要になります。本家ページにも書いてありますが、だいたいのLinuxディストリビューションにはあらかた対応しているかと思います。
ただしGCC4.3以前やARMでの4.6以前だとコンパイルにいろいろ問題があるようです。

気分と環境が整ったらコンパイルとインストールをします。

% ./configure --enable-optimized --enable-jit --enable-targets=x86_64 --prefix=${HOME}/どこどこ
% gmake -j コア数
% make check
% make install

オプションは適当に。–prefixオプションはなくても構いませんが、今回とは別にclangを既にインストールしている、root権限がないのでホームディレクトリ以下にインストールしたい、などの場合には指定しておきましょう。
その場合はインストールディレクトリにパスを通しておくことをおすすめします。
gmake -j の後は、CPUのコア数に応じて数字を指定しておくとそれなりに早くなります。
手元の環境では10分もあれば終わりました。

GCCのインストールがあれこれ依存してて泣きたくなるのに比べたら、LLVMのインストールは非常にシンプルです。この辺りはさすが後発。

test1.cppにHello Worldを書いて確認。しかもあえて間違えて。

#include <iostream>

int main(int argc, char** argv){
  cout << "Hello, World!" << std::endl;
return 0;
}

そしてコンパイル。

% clang++ test1.cpp -o test1_clang
test1.cpp:4:2: error: use of undeclared identifier 'cout'; did you mean 'std::cout'?
        cout << "Hello, World!" << std::endl;
        ^~~~
        std::cout
/***/include/c++/4.7.2/iostream:62:18: note: 'std::cout' declared here
  extern ostream cout;          /// Linked to standard output
                 ^

なかなか賑やかでフレンドリーなエラーが出力されました。(といってもこの程度ならg++でも比較的穏やかなエラーを返してくれますが。)

coutをstd::coutに書きなおして実行。無事Hello Worldです。

---

インストールだけで満足してもつまらないので、円周率の計算でベンチマークをとってみます。円周率の計算ソフトはいろいろありますが、今回は多倍長計算でGMP (The GNU Multiple Precision Arithmetic Library)を導入してGMP謹製の円周率計算ベンチを回してみます。
(なんでBSDライクなライセンスのclang/LLVMを準備したのにGNUプロジェクトのソフトを持ってくるのかという感もありますが、まあGMPはLGPLですし許してください。)

なお、遊びでなく本気のベンチマークが見たい方はPhoronixのページなどを参考にして下さい。
ベンチマークに円周率を使ったのは気分です。

GMPのサイトに行って最新バージョンを持ってきます。現時点では5.1.1が最新。GCCとclangでコンパイルします。

% mkdir gmp
% cd gmp
% wget http://ftp.jaist.ac.jp/pub/GNU/gmp/gmp-5.1.1.tar.bz2
% tar xvfj gmp-5.1.1.tar.bz2
% cd gmp-5.1.1
% ./configure --prefix=${HOME}/どこどこ/gcc-4.7 --enable-cxx CC=gcc CXX=g++
% make
% make check
% make install
% ./configure --prefix=${HOME}/どこどこ/clang-3.3 --enable-cxx CC=clang CXX=clang++
% make
% make check
% make install

GMPの円周率計算ベンチのページからソースコードを持ってきて、適当にpi.cなどの名前で保存します。
そしてコンパイル。-Iと-Lのパスは適当に読み替えてください。

% gcc pi.c -I../gmp/gcc-4.7/include -L../gmp/gcc-4.7/lib -lgmp -lm -O3 -o pi-gcc
% clang pi.c -I../gmp/clang-3.3/include -L../gmp/clang-3.3/lib -lgmp -lm -O3 -o pi-clang3.3
/tmp/pi-fDvkj6.o: In function `bs':
pi.c:(.text+0xc90): undefined reference to `fac_mul2'
pi.c:(.text+0xd89): undefined reference to `fac_mul2'
pi.c:(.text+0xdf8): undefined reference to `fac_mul_bp'
pi.c:(.text+0xeb7): undefined reference to `fac_mul_bp'
pi.c:(.text+0xed6): undefined reference to `fac_mul_bp'
clang: error: linker command failed with exit code 1 (use -v to see invocation)

あれ?clangでリンカエラーが出るぞ?

ソースを見てみるとinline指定されている関数が見つからないことになっているもよう。プリプロセッサの問題かと思って-Eオプションつけてみても無事残っている。
…と調べてみたら本家に説明がありました。曰く、「C99の仕様ではinline指定をした関数定義はインライン展開されたときにのみ有効になるため、inline指定のない定義は別に必要になります。」とのこと。つまり、inlineと書いてあってもインライン展開されるかどうかはコンパイラが自由に決めることであり(clangでは「ゆるい提案」扱いになるとのこと)、インライン展開しないと決めたときに関数未定義になるらしい。

知らない。知らなすぎる。
clangで起きるエラーは実は手元のソースコードの仕様違反であった(がコンパイラの拡張やバグで潜在化していた)ということがよくあるようですが、今回もそうでした。コンパイラを変えてみるというのもそれだけで勉強になります。
まあ、今まで自分が知っていたと思っていた世界は実はC/C++でなくgcc/g++であった、というのが改めてよくわかります。

気を取り直して再コンパイル。inline 関数をstaticにするとかオプションに-std=gnu89をつけるとかいろいろありますが、単純にinlineを外すのが簡単です。どうせinlineつけてもclang最適化のほうが偉いので。

---

計測タイム。環境はXeon E5-2690 x2 (Sandy Bridge, 2.9/3.8GHz 2CPU/16コア/32スレッド)、メモリ16Gbyte、Red Hat Linux 64bit。 Xeonでコアが多いけど今回はマルチコア未対応なので関係なし。あえて言うならメモリ転送速度が4チャネルx2でやたら早い。

あまり変わらない。桁をもっと増やしてみる。

やっぱりほとんど変わらない。まあ、GMPが相当最適化されていてコンパイラの入る余地も少ないのでしょう。
個人的にはGMPをコンパイルできただけで私はだいぶ満足です。あと、円周率10億桁なんてとんでもないものが即座に計算できる時代になったことを再認識して喜んでいます。(clang関係ないけど)

---

それにしてもclangの進化はすごいですね。オープンソースでC/C++のコンパイラを作るなんて既にデファクトスタンダードのGCCもあるしで当分微妙なままなんじゃないかと思っていましたが、大間違いでした。セルフホスティングになってからあれよあれよという間にC++11完全対応まで来てしまいました。

clangはコンパイラとしてだけではなく、統合開発環境との連携も目玉のひとつです。今後Eclipseなど既存のIDEと深く関わっていけばコーディングがより気楽になるのでしょうか。

私も手元の遊びで書くC++はclangを使ってもいいかなと思っています。何よりエラーが見やすくて助かります。

以前ここで書きかけた剛体振り子のプログラムですが、いろいろと紛失してしまったので新しく作りました。

ソフトのダウンロードはこちらから pendulums

ソースコードはgithubからどうぞ https://github.com/So-Takamoto/triple_pendulums/
手元の環境はVisual Studio 2010です。

(以前の記事: [昔話] 振り子のシミュレーションをした話)

例によってC++, OpenGL+freeglutで作成しています。
キーの1,2,3で微分方程式を解く方法を変更でき、iで初期化、lで軌跡の表示切り替え、qで終了です。

質量や重力加速度、初期条件、各振り子の長さを変えたい場合はソースコードをいじってみてください。
64bit版でコンパイルを行いたい場合はfreeglutのライブラリを64bit版に変更してください。今回は軽いのであまり必要はないかと思いますが。

例によって動画を用意しておきます。

なんでこんなシミュレーションをしようと思ったかの話は以前書いた(早い話がこういうおもちゃが大学に置いてあった)ので、今回は原理とか内容の紹介メインにしたいと思います。

今回扱う系は平面上の棒3個なので、そのままで言えば3自由度×3個=9自由度の系なのですが、各ジョイントが2自由度分の拘束を持つため結局本来の自由度は3になります。

ひとつのやり方は、これらのジョイントにかかる力を陽に扱う方法です。利点は自由度を実座標と対応付けられるため意味がわかりやすいことと、比較的機械的な処理で立式できることです。
これらジョイントにかかる力は幾何学的な関係から導出することもできますし、今回のような単純なリンク機構であればエンドエフェクタ(一番先の振り子)から芋づる式に数値解析の範囲で導出することもできます。ロボットハンドなどの制御ではこの辺りの話が動力学計算としてたくさん出てきます。

一方で今回のような系ならそこまで面倒な話をしなくても、解析的に運動方程式を求めることができます。こういう系ではラグランジュ形式を出発点にすると便利です。ジョイントに働く力を陽に扱わずに、自由度の数だけ運動方程式を考えればいいからです。

---

ラグランジュの運動方程式の出発点は、ラグランジアンL(=T-U)の表現です。運動エネルギーTは重心の並進運動と回転運動から、位置エネルギーUは重心の位置から求まります。

3つの剛体を根本から順に0,1,2と置いて(0から始めるのはC言語使いの癖です)、それぞれ棒の角度を、長さを、質量を、根元側のジョイントからの重心の相対位置をと置きます。3自由度はここでは3つのθに対応します。

このとき、相対位置と相対速度について、

また、絶対位置について、

速度も上式を微分するだけです。欲しいのはなので、 展開して内積を計算すればそれぞれ導出できます。

回転のエネルギーについてはと慣性モーメントから計算できます。長さlの棒の慣性モーメントIはなので、回転のエネルギーについてはと求まります。

なんやかんやで計算すると以下のようになります。

汚いようですがθに関しては対称的な形をしています。

同様に位置エネルギーに関しても

ラグランジアンLは先に書いたようにL=T-Uとして求まります。あとは ラグランジュの運動方程式

を解けばについて求まるという算段です。ただし、上の式を解こうとするとについての3元連立方程式になるので、微分方程式を解く前に適当な方法で変数を独立させておく必要があります。3変数なら一応解析的にも書き下すことができますが、大変面倒なので今回はガウス・ジョルダンの消去法を使いました。

ラグランジュ形式は便利ですが、最後に出てくる式がいまいち意味付けしにくいという不思議なところがあります。もっとも世の中の物理というのは別に人間にわかってもらうために運動しているわけではないのですが。
そもそもラグランジアンL自体なんの値なのかと考えだすと悩みどころはどんどん増えていきます。教科書には「意味のある値ではない」ときっぱり書いてあったりします。清々しいですね。

まあそれでもハミルトン形式まで持っていけばハミルトニアンHに意味を持たせることができるのですが。

---

運動方程式が解けたら次は微分方程式の数値計算の出番です。

微分方程式を数値計算で解く方法は数多くあります。今回のプログラムでは3種類から選べるようにしています。すべて陽解法で、それぞれ1次、2次、4次の精度を持っています。

おそらく最もシンプルでわかりやすい方法がオイラー法で、離散化した点の間を直線で補間する方法です。直線なので1次の精度がある、と表現されます。
オイラー法は精度が悪いため一般には使われません。左上のエネルギーの値のぶれが大きいことからもわかると思います。

より高い次数まで精度を持たせた方法があり、例えば4次のルンゲ・クッタ法が有名です。今回は2次のホイン法も用意しました。エネルギーの値を見れば違いがわかると思います。

ただし、たとえ4次精度を持っていたとしても長期間のシミュレーションを行うと次第にエネルギーが変化していってしまいます。長期間のシミュレーションでこういう性質があると結構問題で、例えば太陽系のシミュレーションをしていたら計算誤差で惑星が吹っ飛んでいったとか、分子の運動をシミュレーションしていたら計算誤差で勝手に沸騰した、なんていうよくわからないことが起きてしまいます。こういう場合にはエネルギーが長期間で安定するよう工夫した数値計算を行う必要があります。シンプレクティック性を持つ計算などと呼ばれますが、いちいち書いていたらとんでもない量になるのでここでは割愛します。

私の研究分野でよく使われる速度ベルレ法などの計算手法はうまい具合にこのシンプレクティック性を持っており、2次精度の手法ですが、計算が簡便でしかも陽解法でありながら長期間でエネルギーが変化しにくいという便利な性質を持っています。

---

系としてはとっても単純でありながら、運動はとても不規則です。時には大きくゆったり回ったり、時には先の振り子だけぐるぐる回ったり、それらを予告なく切り替え続けます。

このような運動はカオスと呼ばれますが、単純な法則から複雑な運動が生まれる例は他にもたくさんあって、例えば前回の反応拡散系も支配している方程式は単純ながら模様が複雑になったりもします。

それなら逆に、世の中の一見複雑な現象も実は単純な法則によって成り立っているのではないか?と考えるものの見方があります。では複雑な運動から単純な法則は見つけることができるのか?振り子ではなく、例えば世の中の気候の変化、株価、人の行動のような複雑な現象は実は単純な法則から成り立っているのではないか?そうだとしたらその法則は何か?という疑問に答えるのが、「複雑系」と呼ばれる学問のひとつの目標なのではないか、と私は勝手に考えています。

 

まあ、私は複雑系の研究者ではないので、好き勝手に妄想しているところが大きいのですが。

我々には少し先までしか見えないが、そこになすべきことをたっぷりと見付けることができる。
―アラン・チューリング

以前，セル・オートマトンに関係するプログラミングで作った反応拡散方程式シミュレーションのお話です．
眠らせておくのももったいないので，ここで紹介しておきます．

反応拡散方程式というのは，ある空間上での何かの物質の濃度の変化を表現した微分方程式です．もともとは化学の分野で反応をモデル化する際に考えられたものですが，その解が複雑な振る舞いをすることから様々な分野で興味の対象となり，実世界の模様との関連性が指摘されたりもしています．

今回はこの反応拡散方程式のシミュレータの紹介です．まずは動画をどうぞ．

この模様はチューリングパターンと呼ばれ，動物の模様や植生の分布などに類似の模様が見られることがあります．

今回扱う反応拡散方程式は1980年代に提案されたGray-Scottモデルと呼ばれるもので，以下の形に書き下されます．

 

uとvという2種類の物質(?)についての式で，左辺が変化の度合い，右辺第1項が拡散，第2項が反応，第3項が生成と消滅に対応しています．Fとかkとかはパラメータです．

解析的に解けなさそうな式ですが，差分方程式にしてしまえばコンピュータで解くのはそれほど困難ではありません．得られるパターンが独特かつパラメータによって多様に変化するので，数値解析の例題としてはなかなか楽しい部類に入るのではないかと思います．

私がこの反応拡散方程式にはじめて出会ったのは，まともにプログラミングを学んだ少し後の頃でしたが，当時紹介されていた動画などを見てすっかり気に入ってしまった記憶があります．

1つは「ぬるぬる」している点です．特にはじめて触った頃のプログラミングというのはテキストや数式を扱うことが多いため，どうにも真四角でお固い世界だと錯覚しがちです．こんな「ぬるぬる」は当時の私にとってあまりコンピュータっぽくなかったというのがあります．
もう1つは，式が簡単であることです．すごく複雑なものが複雑な式から出てくるならまあわかりますが，この式自体はそれほど複雑ではありません．
各部分では単に隣接した点の情報のみから変化しているだけなのですが，それでも全体としてパターンが現れるというのはなかなか面白いと思います．

プログラムを置いておきます．手元の 64bit Windows 用（7で動作確認済み）なので，実行できない方は一番下のソースコードからコンパイルしてみて下さい．
マウスでドラッグした部分が消せます．キーボードの左2列（1,2,q,w,a,s,z,x)でパラメータの変更，iキーで分布を初期化できます．

rd_system_2 [追記： 録画用でウィンドウが小さいままだったので大きくしました．]

rd_system_32bit [さらに追記： 32bit版も用意しました]

---

私の大学の機械工学科には「シミュレーション & ビジュアライゼーション演習」という授業があります．この授業では何らかのシミュレーションの可視化に挑戦します．

実際に扱う対象は連続体というよりはどちらかというと剛体系の人が多いかと思いますが，何にしても可視化してみるというのは楽しいもので，あれこれ計算した結果が絵になって出てくる達成感がある気がします．
私も昨年度TAで参加したときに楽しませてもらいました．

プログラミングでちょっと遊んでみたいという人は，マイコンのような動くものにつなげてみるか，あるいはこういう可視化のできるもので何か作ってみると手ごたえを感じやすいと思います．

---

実装はC++/OpenGL+freeglut，方程式は1次の陽解法です．私はこういうものを作るときはたいていC# & 速度の必要な部分のみC++で実装してしまうのですが(実はこれもそうでした)，今回は公開する前に移植しやすいようにC++で書き直してみました．

ソースはgithubにあげておきます．
githubを使うのははじめてなので，何か不具合があったら申し訳ございません．

ソース

知性は、方法や道具に対しては鋭い鑑識眼を持っていますが、目的や価値については盲目です。 – アインシュタイン

大学院の夏学期の授業で自由課題が出たのですが，その提出物を授業で発表することになりました．

課題で愉快な結果が出たわけでもないのですが，せっかくなのでそのとき作ったプログラムとその結果を書き留めておきます．
テーマは「モンテカルロ法」．何か身近な問題についてシミュレーションを行ってこいとのことでした．

モンテカルロ法は乱択アルゴリズムの一種です．ざっくり言ってしまえばたくさんやってみて統計をとるという方法ですね．
コンピュータのくせに「答えは厳密ではないが，たぶんだいたい合ってる」という結果の出方が大好きです．

さて，たくさんのものをぼんやり見る…となるとやはり気になるのは人の群れのシミュレーションです．人って群れると不思議な振る舞いをしますよね．
研究となると面倒な分野ですが，授業の課題なら気にせず変なことができます．

それで扱う現象ですが，人って誰かの隣には座らないんですよね．電車などに乗っているとよく見る景色です．
これは映画館とか大学の講義でも同じで，隣が開いてる席があったらそっちを使うというのが人間というものの性質のようです．

男子トイレでも同じことが起きますね．大きなトイレでがら空きなら誰かの隣の便器というのは普通使いません．もしがらがらのトイレで隣に人がやってきたら，何か他の目的があるのではという考えが頭をよぎります．

で，あのふるまいのシミュレーションをしてみようと思ったわけです．

---

トイレのモデルは以下のようになります．

便器の数をn個として，そこに人が次々と到着します．各個人の便器使用時間はガウス分布，到着時間はポアソン到着としています．

各個人は両隣が空いている便器があればそれを，なければそれ以外の便器を使います．全部埋まっていたら待ち行列に並びます．待ち行列の人は便器が空き次第詰めていきます．
同程度の魅力の便器があればランダムに選びます．

パラメータは便器数n, 便器使用時間の平均μおよび標準偏差σ，それと平均到着時間間隔λとなります．
便器使用時間についてはどこかにデータがあればよかったのですが，見つからなかったので自分で1週間計測しました．
平均は20秒，標準偏差は3秒となりました．今回は両側に3秒追加して26±3秒としてシミュレーションを行いました．

---

n=5として平均到着時間間隔λを変化させた結果が以下になります．右の番号はトイレの便器ですね．端から順に1から5です．

横軸は到着頻度（λ）になります．左端は人が来ない状況で，右に行くにつれてトイレが賑わっていきます．右端は常にトイレが使われているという相当やばい状態ですね．現実にこんなトイレがあったら早いところ増設したほうが事故を未然に防げていいと思います．

人がとてもまばらな時は，どの便器も同様に使われます．また相当混んでいる時も同様に使われます．
まあ他人がいなければ区別なく使いますし，混んでる時は選ぶ余地もなく全部使われてるので当然とも言えます．

問題はこの中間のときで，ちょうど限界の半分くらいの人が押し寄せるときに便器の使用頻度に差が生じます．
端のトイレが比較的よく使われる一方で，端から2番目のトイレ（2,4）の使用頻度は落ちています．おそらく端っこは隣が片方しかないため，そこそこ混んでいても入りやすいのだと思います．それでその隣は使われないのでしょうね．

以下が便器の数n=31とした結果です．大きなトイレですね．

多くの便器達はおおよそグラフの中央に集まっています．これは便器ごとに使用率が変わらないことを意味しています．
一方で上に大きく曲がった線が2本あります．これが両端の便器です．使用率は他のものに対して有意に高くなっています．
下の2本は端から2番目の便器になります．こちらは使われにくくなっていますね．

現実問題を表現できていますかね?
今度駅とかのトイレに行ったときは便器の使用状況に気をつけてみるのもいいかもしれません．

---

この結果から何がわかるか?というと，まあせいぜい端の便器は汚れてるかもしれないよ，ということくらいです．
まあトイレに限らずなんでもそうかもしれません．大学の教室の椅子でも5席あったら1,5番目，詰まったらその後3番目から埋まりますね．

これ自体は非常に単純なシミュレーションですが，簡単な仮定をおくだけでこのようなデータが得られるのがシミュレーションの面白いところです．
気軽にコードが書けるとこういった遊びの幅も広がります．

コメントも書いてないものですが一応ソースコードを置いておきます．Visual Studio 2010で作ったものですが単純なC++のコードなのでだいたいの環境で動くと思います．

restroom

発表に使った資料も置いておきます．後半は違う話になっていますがいずれここで書けたら面白いかなと思っています．

マルチスケール課題1

シミュレーションの内容に補足をすると，今回はキネティックモンテカルロ法という手法になっています．これは扱う時間のスケールに対して稀な現象を含めて扱うときに使われる手法のひとつで，系で起きるイベントをリストアップし，次のイベントまでにかかる時間を計算してイベント単位で時間を進めていく方法です．

また，このコードにはトイレの待ち行列が伸びすぎるとクラッシュする機能があります．シミュレーションとはいえあんまり我慢させるのはかわいそうですから．

ちなみにトイレシミュレータの原案は昨年の研究室の先輩方でした．先輩方ありがとうございました．