module Main where
 
import Data.Array (Array, accumArray, elems, Ix, inRange)
import System.Random (RandomGen, newStdGen, random)
import Control.Monad.State (State, state, evalState)
 
histogram :: Ix i => (i, i) -> [i] -> Array i Int
histogram bounds bins = accumArray (+) 0 bounds [(x, 1) | x <- bins, inRange bounds x]
 
bernoulli :: RandomGen gen => Double -> State gen Int
bernoulli p =
  do x <- state random
     return $ if x < p then 1 else 0
 
binomial :: RandomGen gen => Int -> Double -> State gen Int
binomial count p =
  do ns <- sequence $ replicate count $ bernoulli p
     return $ sum ns
 
simulate :: RandomGen gen => gen -> Int -> State gen a -> [a]
simulate gen tries m = evalState (sequence $ replicate tries m) gen
 
test :: Double -> Int -> Int -> IO ()
test bias nails tries =
  do gen <- newStdGen
     let bins = simulate gen tries $ binomial nails bias
     let h = histogram (0, nails) bins
     print $ elems h
 
main :: IO ()
main = test 0.5 100 10000
 

bW9kdWxlIE1haW4gd2hlcmUKCmltcG9ydCBEYXRhLkFycmF5IChBcnJheSwgYWNjdW1BcnJheSwgZWxlbXMsIEl4LCBpblJhbmdlKQppbXBvcnQgU3lzdGVtLlJhbmRvbSAoUmFuZG9tR2VuLCBuZXdTdGRHZW4sIHJhbmRvbSkKaW1wb3J0IENvbnRyb2wuTW9uYWQuU3RhdGUgKFN0YXRlLCBzdGF0ZSwgZXZhbFN0YXRlKQoKaGlzdG9ncmFtIDo6IEl4IGkgPT4gKGksIGkpIC0+IFtpXSAtPiBBcnJheSBpIEludApoaXN0b2dyYW0gYm91bmRzIGJpbnMgPSBhY2N1bUFycmF5ICgrKSAwIGJvdW5kcyBbKHgsIDEpIHwgeCA8LSBiaW5zLCBpblJhbmdlIGJvdW5kcyB4XQoKYmVybm91bGxpIDo6IFJhbmRvbUdlbiBnZW4gPT4gRG91YmxlIC0+IFN0YXRlIGdlbiBJbnQKYmVybm91bGxpIHAgPQogIGRvIHggPC0gc3RhdGUgcmFuZG9tCiAgICAgcmV0dXJuICQgaWYgeCA8IHAgdGhlbiAxIGVsc2UgMAoKYmlub21pYWwgOjogUmFuZG9tR2VuIGdlbiA9PiBJbnQgLT4gRG91YmxlIC0+IFN0YXRlIGdlbiBJbnQKYmlub21pYWwgY291bnQgcCA9CiAgZG8gbnMgPC0gc2VxdWVuY2UgJCByZXBsaWNhdGUgY291bnQgJCBiZXJub3VsbGkgcAogICAgIHJldHVybiAkIHN1bSBucwoKc2ltdWxhdGUgOjogUmFuZG9tR2VuIGdlbiA9PiBnZW4gLT4gSW50IC0+IFN0YXRlIGdlbiBhIC0+IFthXQpzaW11bGF0ZSBnZW4gdHJpZXMgbSA9IGV2YWxTdGF0ZSAoc2VxdWVuY2UgJCByZXBsaWNhdGUgdHJpZXMgbSkgZ2VuCgp0ZXN0IDo6IERvdWJsZSAtPiBJbnQgLT4gSW50IC0+IElPICgpCnRlc3QgYmlhcyBuYWlscyB0cmllcyA9CiAgZG8gZ2VuIDwtIG5ld1N0ZEdlbgogICAgIGxldCBiaW5zID0gc2ltdWxhdGUgZ2VuIHRyaWVzICQgYmlub21pYWwgbmFpbHMgYmlhcwogICAgIGxldCBoID0gaGlzdG9ncmFtICgwLCBuYWlscykgYmlucwogICAgIHByaW50ICQgZWxlbXMgaAoKbWFpbiA6OiBJTyAoKQptYWluID0gdGVzdCAwLjUgMTAwIDEwMDAwCiA=