πの計算

ライプニッツの公式

---

	import Data.Ratio
	a = cycle [1,-1]
	b = [1,3..]
	c = zipWith (\x y -> x * 1 % y) a b
	main = do
	let a = take 1000 c
	let b = sum a
	print $ fromRational $ b * 4

view raw l.hs hosted with ❤ by GitHub

---

実行結果
＄　runghc ｌ.hs
3.140592653839793

マチンの公式によるπ

---

	n=10^10^4
	m0 = 5
	m1 = 239
	g n m = let a = scanl f (n `div` m) [5,9..]
	b = scanl f ((n `div` m ^ 3) `div` 3) [7,11..]
	in sum (takeWhile (>0) a) - sum (takeWhile (>0) b)
	where f x y = (x `div` m^4)*(y-4) `div` y
	main = do
	let a = g n m0
	let b = g n m1
	let r =(a * 4 - b) * 4
	print r

view raw m.hs hosted with ❤ by GitHub

参考にしたページ
http://handasse.blogspot.com/2011/12/1.html

---

Parallel版

	import Control.Parallel
	import Control.Parallel.Strategies (rpar, Strategy, using)
	import Text.Printf
	n=10^18^4
	m0 = 5
	m1 = 239
	g n m = let a = scanl f (n `div` m) [5,9..]
	b = scanl f ((n `div` m ^ 3) `div` 3) [7,11..]
	in sum (takeWhile (>0) a) - sum (takeWhile (>0) b)
	where f x y = (x `div` m^4)*(y-4) `div` y
	main = print $ (a * 4 - b) * 4
	where (a,b) = (g n m0, g n m1) `using` parPair
	parPair :: Strategy (a,b)
	parPair (a,b) = do
	a' <- rpar a
	b' <- rpar b
	return (a',b')

view raw pm.hs hosted with ❤ by GitHub

1万桁では差がでないので
１０万桁ほど計算します。

コンパイルします。

ghc  --make -O2 pm.hs -rtsopts -threaded -eventlog

CPUはCore2Duo。デュアルコアなので　-N2 オプションで実行します。
$./pm +RTS -N2 -l
threadscope で動作を確認します。
$ threadscope pm.eventlog

僅かですが速くなっています。

並列計算を効率よく利用するには最初から並列化を意識してコードを書く必要があるようです。

---

spigot アルゴリズムをLoopで書いたコード
kk62526　さんが１０進Basicで書いたコードをそのまま移植しています。

	import Control.Monad.State
	loop = do
	(q,r,t,i,u,y) <- get
	liftIO $ putStr $ show y
	put(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1,u,y)
	(q,r,t,i,u,y) <- get
	put (q,r,t,i,3 * (3 * i + 1) * (3 * i + 2),y)
	(q,r,t,i,u,y) <- get
	put (q,r,t,i,u,(q * (27 * i - 1) + 5 * r) `div` 5 `div` t)
	loop
	main = do
	let (q,r,t,i) = (1,180,60,2)
	let u = 3 * (3 * i + 1) * (3 * i + 2)
	let y = (q * (27 * i - 12) + 5 * r) `div` 5 `div` t
	runStateT loop (q,r,t,i,u,y)

view raw spi.hs hosted with ❤ by GitHub

グローバル変数　「i(注）」に現在何桁目を計算しているかが保持されてますので（ただし＋２された数字）
(注：　本当は無名ですが、こう捉えたほうが分かりやすいです。）
if (i-2) `mod` 100 == 0　then ・・・・・・
のようなコードを挿めば１００桁ごとに区切ったりできます。

	import Control.Monad.State
	f = do
	(q,r,t,i,u,y) <- get
	put(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1,u,y)
	(q,r,t,i,u,y) <- get
	put (q,r,t,i,3 * (3 * i + 1) * (3 * i + 2),y)
	(q,r,t,i,u,y) <- get
	put (q,r,t,i,u,(q * (27 * i - 1) + 5 * r) `div` 5 `div` t)
	loop0 = liftIO (print 3) >> f >> loop
	loop = do
	(q,r,t,i,u,y) <- get
	liftIO $ putStr $ show y
	if (i-2) `mod` 100 == 0
	then do liftIO (print ())
	liftIO $ putStrLn $ show $ i - 2
	f
	loop
	else f >> loop
	main = do
	let (q,r,t,i) = (1,180,60,2)
	let u = 3 * (3 * i + 1) * (3 * i + 2)
	let y = (q * (27 * i - 12) + 5 * r) `div` 5 `div` t
	runStateT loop0 (q,r,t,i,u,y)

view raw p.hs hosted with ❤ by GitHub

実行結果
 $ runghc p.hs > a
$ less a
3 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679()
100 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196()
200 4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273()
300
~
~
~
（略）
ちなみに１００００桁を計算し終えるのに２秒少々でした。
ギボンズ博士が spigot アルゴリズムに対し
"not competitive with state-of-the-art arithmetic-geometric mean algorithms"
と言うとおり決して高速とはいえませんが、いきなり１桁目から表示を始めますので、
体感的には超高速です。
(注：１００００桁の計算では８万数千桁のIntegerが３つ、作られます。）

ギボンズ博士のspigot アルゴリズムは
https://www.cs.ox.ac.uk/jeremy.gibbons/publications/spigot.pdf
コードだけ紹介しますと

	piG3 = g(1,180,60,2)
	where
	g(q,r,t,i) = let (u,y)=(3*(3*i+1)*(3*i+2),div(q*(27*i-12)+5*r)(5*t))
	in y : g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
	main = print piG3

view raw sppi.hs hosted with ❤ by GitHub