円周率計算方法

無理数の一つである円周率の計算方法には昔からさまざまな方法が知られている．最も有名なものはモンテカルロ法だろうか．
これは一辺の長さが2rの正方形の中に，半径rの縁を描き，$(x,y), x,y \in [-2,2]$となる乱数座標(x,y)を大量にプロットし，全点数に対する円内部にプロットされた数の比率から円周率を求める方法である．そのほかにも，ビュフォンの針などさまざまな方法が知られている

剛体運動から導出される円周率

そして，2005年に物理現象から円周率を求める論文が出ていたことに最近になって気がついた．解説記事はQiita上にも複数投稿されているので詳細はそちらに譲る．
https://www.maths.tcd.ie/~lebed/Galperin.%20Playing%20pool%20with%20pi.pdf

質量m, Mの二つの物体があり，$M=m\cdot 100^n$とする．これらの物体を，壁・質量m物体・質量M物体の順に並べ，M物体が壁方向にある速度で運動することを仮定する．
この時，質量m物体に衝突した際，運動量の一部がm物体に移動し，m,M両物体が壁に向かって運動する．m物体の方が速度が速いため先に壁に衝突してM物体方向へ跳ね返ってくる．以後，これらの運動が繰り返される．
この一連の運動で生じる壁および　m，M物体間の衝突回数を合算すると，円周率がn-1桁まで求まるという内容だ．

簡単な運動で円周率がもとまるとは面白いので，haskellで簡単に書いてみた．
浮動小数点周りの計算がいけておらず，まだ小数点以下2桁までしか計算できないが，確かに衝突回数から円周率がもとまっている

以下にコードを掲載する．

import Debug.Trace

data Body = Body {pos :: Pos, mass :: Float, vel :: Float} deriving (Eq)
type Pos = Float

instance Show Body where
  show b = show (pos b) ++ " [m] " ++ show (vel b) ++ " [m/s]"

nextPos :: Body -> Body
nextPos x = x {pos = pos x + vel x * dt}

reflect :: Body -> Body
reflect b = b {pos = pos b + v' * dt, vel = v'} where
  v' = negate (vel b)

collision :: (Body, Body) -> (Body, Body)
collision (b1,b2) = (b1 {pos = pos b1 + v1'*dt, vel = v1'}, b2 {pos = pos b2 + v2'*dt, vel = v2'}) where
  v1' = ((mass b1 - mass b2) * vel b1 + 2 * (mass b2) * (vel b2)) / (mass b1 + mass b2)
  v2' = (2 * mass b1 * vel b1 + (mass b2 - mass b1) * (vel b2)) / (mass b1 + mass b2)

init_b1, init_b2 :: Body
init_b1 = Body {pos = -100, mass = 10^4, vel = 1}
init_b2 = Body {pos = 0, mass = 1, vel = 0}
wall = 100  -- x = 100の位置に壁
dt = 0.0001 -- 1ms

count :: (Body, Body) -> Int
count (b1,b2)
  | pos b1 > wall = trace ("stop" ++ show' b1 b2) 1
  | vel b1 <= 0 && vel b2 <= 0 && vel b1 <= vel b2 = trace ("finish" ++ show' b1 b2) 0
  | pos b2   >= wall   = trace ("wall " ++ show' b1 b2) $ 1 + count ( b1, reflect b2 )
  | pos b1   >= pos b2 = trace ("collision " ++ show' b1 b2) $ 1 + count ( collision (b1, b2) )
  | otherwise = count (nextPos b1, nextPos b2)

show' b1 b2 = show b1 ++ ", " ++ show b2


main = do
  putStrLn "衝突回数"
  putStrLn . show $ count (init_b1, init_b2)

年末の隙間時間に遊ぶのにちょうどよかった．
より大きい桁が計算できるよう改善していきたい