水曜日, 8月 19, 2009

パサデナで円周率

 Pythonを用いた科学計算に関する学会に行ってきた[link]。開催地はパサデナにあるCaltech (カリフォルニア工科大学)。バスと電車を乗り継いで行かなくてはならないため、5時起床で6時に家を出たのだが、結局着いたのは10分遅刻の8時40分。最寄り駅から歩けるものと勝手に思っていたら、とんでもない、4km以上はあった。広すぎるぜカリフォルニア。
 専門は計算生物学だが、計算の学会に出たのは実は初めて。プログラミングは独学で学んだ。今日と明日はチュートリアルがあって、僕は基本クラスに参加した(参加者は圧倒的に上級クラスの方が多かった)。午前中はIPythonを使いながら基本的な文法の確認。さすがにこれはやさし過ぎた。午後は一番聞きたかった多次元配列のライブラリNumpyと描画ライブラリMatplotlibの入門。実際に開発に携わっている人が講義してくれたので、マニュアルでは分からないことが聞けて良かった。
 今日は円周率を求める演習をやったので、そのPythonによる実装を示します(Wallisの式による計算(上式))。分数をその都度計算すると端数をばっさり切って値が不正確になるので、分子と分母を先に計算してから割り算する必要がある。これをCで素直にやるとオーバーフローしてしまう(つまり、どんどん掛け算をしていくと、とてつもなく大きな数字になってしまう)。Pythonの場合は、一行目のおまじないのおかげで素直に計算ができる。
 そういえば先日、筑波大のスーパーコンピュータが円周率の計算で世界一になりました[link]。円周率計算はコンピュータ進化のベンチマークなのだそうです。
from __future__ import division

def wallis():

  num = den = my_pi = 1

  for i in range(1, 10000):

    tmp = 4 * i ** 2

    num *= tmp

    den *= tmp - 1

  my_pi = 2 * (num / den)

  return my_pi


if __name__ == "__main__":

  print wallis()

0 件のコメント: