モンテカルロ法による円周率の計算 (ベンチマーク用)

test-monte-carlo-pi.el

;;; -*- lexical-binding: t; -*-

(defvar rand-next 1)

(defsubst myrand ()
  "C 言語版や Common Lisp 版等と計算結果を合わせたい為に 32 bit の整形合同法を実装してます…"
  (setq rand-next (logand (+ (logand (* rand-next 1103515245) #x7fffffff) 12345) #x7fffffff))
  (logand (ash rand-next -16) #x7fff))

(defsubst float-rand ()
  (/ (float (myrand)) 32768))

(defun calc-pi ()
  (let ((count 0)
        (max 10000000)
        x y z)
    (dotimes (_ max)
      (setq x (float-rand)
            y (float-rand)
            z (+ (* x x) (* y y)))
      (when (< z 1.0)
        (setq count (1+ count))))
    (* (/ (float count) max) 4)))

(let ((time (current-time))
      my-pi)
  (setq my-pi (calc-pi))
  (print (float-time (time-since time)))
  (print my-pi))