ビュフォンの針で交差数を求めるアルゴリズム

Pythonでプログラミングを学習中に，ビュフォン(Buffon)の針のシミュレータを実装したのですが，そのときのアルゴリズムについてメモ。

# ビュフォン(Buffon)の針とは
等間隔に並んだ直線の上に針をたくさん落とし，直線と交わっている針の割合をもとに円周率$\pi$を推定しようというものです。

今，直線の間隔$d$，針の長さを$l$とします。
<img width="337" alt="キャプチャ.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/280743/0faf3166-f912-2837-7cd5-c31cd5be46f9.png">
一般に$l \le d$のとき，針を落とした時に直線と交わる確率は$\frac{2l}{\pi d}$になることが知られています。よって針を$n$本落としたとき，直線と$m$本が交わったとすると，
$$\frac{m}{n} = \frac{2l}{\pi d}$$
と近似でき，このことから
$$ \pi = \frac{2nl}{md}$$と円周率が計算できます。

詳しくは[Wikipedia](https://ja.wikipedia.org/wiki/ビュフォンの針)等を参照してください。

# 準備
さて，ビュフォンの針を実装していきます。`BuffonCal`というクラスにまとめます。
### 初期設定

```python

import numpy as np

class BuffonCal:
     def __init__(self, length, interval, width_x=1000, width_y=1000):
         self.length = length
         self.interval = interval
         self.width_x = width_x
         self.width_y = width_y
```

`length`が針の長さ$l$，`interval`で直線の間隔$d$を表すことにします。
`width_x`と`width_y`は針を落とすキャンパスを記しています。

### 直線のy座標を返すメソッド

```python

    def base_lines(self):
        self.lines = (i * self.interval
                      for i in range(self.width_y//self.interval + 1))
        return self.lines
```

直線のy座標をイテレータの形で返します。

### 針を落とすメソッド

```python
     def drop(self, number):
         self.number = number
         self.drops_x1 = self.width_x * np.random.rand(self.number)
         self.drops_y1 = self.width_y * np.random.rand(self.number)
         theta = np.pi * np.random.rand(self.number)
         self.drops_x2 = self.drops_x1 + self.length * np.cos(theta)
         self.drops_y2 = self.drops_y1 + self.length * np.sin(theta)
```

針の端点の座標を$(x_1,y_1)$，$(x_2, y_2)$として，まず$(x_1,y_1)$を乱数で生成し，さらに偏角$\theta$を乱数で生成して$(x_2, y_2)$を決定しています。必ず$y_1 \le y_2 $になるようにしてあります。
各座標ごとに`numpy.ndarray`型として格納しています。

# 落とした針のうち，直線と交わっている本数を数える

さて，落とした針のうち，直線と交わっている針の本数を数えていきましょう。最も素直なのは以下のプログラムでしょうか。

```python

    def count(self):
        self.counts = 0
        for y1, y2 in zip(self.drops_y1, self.drops_y2):
            for line_y in self.base_lines():
                if y1 <= line_y <= y2 :
                    count += 1
                    break
        return self.counts
```

いわゆる**総当たりプログラム**で，針が線をまたいでいるとき，$y_1 \le (線のy座標) \le y_2$となる(線のy座標)が存在するので，そのようなものがあるかを`for`文を用いて全探索しています(このプログラムを採用する場合はNumpyではなくリストを使ってください)。 

それに対し，今回使ったプログラムは以下になります。


```python

    def count(self):
        drops_y1_int = self.drops_y1 // self.interval
        drops_y2_int = self.drops_y2 // self.interval
        self.counts = len(np.where(drops_y1_int != drops_y2_int)[0])
        return self.counts
```

2~3行目で，各$y1,y2$を直線の間隔$d$で割った商を求めています。
**もし，針が直線と交差している場合，商は異なるはずだし，交差していない場合は同じ商になるはず**なので，4行目で商が異なるものの数を計算しています。`for`文を回すより速いと思います。

ちなみに，Pythonの「//」は，割る数が小数であっても計算可能です。以下に例を挙げます。


```

>>> 5 // 2
2
>>> 5 // 2.1
2.0
>>> 7.4 // 0.9
8.0
```


# 残りの処理

$\pi$の推定値を返すコードも記載しておきます。

```python

    def pi_cal(self):
        if self.counts > 0:
            self.pi = (2*self.length*self.number) / (self.counts*self.interval)
        else:
            self.pi = float("inf")
        return self.pi

    def difference(self):
        return '{:+}'.format(self.pi - np.pi)
```

# 終わりに
このプログラムは最終的にtkinterを使って，簡単なGUIアプリケーションにまとめました。コードは[こちら](https://github.com/mo-mo-666/Buffon)。