Diagramsを使って円周率をモンテカルロ法で計算する

プログラミング

Using queries
Diagramsのqueryを使うと、点が図形の内部にあるか外部にあるかが判別できる。

> {-# LANGUAGE NoMonomorphismRestriction #-}
> 
> import System.Random
> import Data.Monoid
> import Control.Applicative
> import Diagrams.Prelude
> import Diagrams.Backend.Cairo.CmdLine
> 
> main = putStrLn . show $ calcPi points

`gens`は乱数のタプルのリストを作っている。何かもっといい方法はない？

> tmap f (a,b) = (f a, f b)
> gens = uncurry zip $ tmap randoms $ split . mkStdGen $ 100

`c`は半径1の円。

> c :: Diagram Cairo R2
> c = circle 1

`points`は乱数から生成された点のリスト。

> points = map p2 $ take 10000 $ gens

`sample c pt`で点`pt`が図形`c`の内部にあるかどうかが分かる。

> count p = length . filter (== Any True) $ sample c <$> p

`calcPi`で円周率の計算。(内部の点の個数) / (全部の点の個数) * 4

> calcPi p = (fromIntegral . count $ p) / (fromIntegral . length $ p) * 4

実行結果:

monte$ ghc monte.lhs
[1 of 1] Compiling Main             ( monte.lhs, monte.o )
Linking monte ...
monte$ ./monte 
3.146
monte$

(Haskell力が足りない……)