Almkvist-Giullera formula for pi: Difference between revisions

Line 585:

import Text.Printf

⚫

iterations = 52

⚫

-- ~~Factorial~~ for arbitrary-precision integers

⚫

main = do

facInteger :: Integer -> Integer

⚫

printf "N. %44s %4s %s\n"

⚫

facInteger n = if n `leInteger` 1 then 1 else n `timesInteger` facInteger (n `minusInteger` 1)

⚫

"Integral part of Nth term" "×10^" "=Actual value of Nth term"

⚫

forM_ [0..9] $ \n ->

⚫

-- ~~Exponentiation~~ for arbitrary-precision integers

⚫

printf "%d. %44d %4d %s\n" n

powInteger :: Integer -> Integer -> Integer

⚫

(almkvistGiulleraIntegral n)

powInteger 1 _ = 1

(tenExponent n)

powInteger _ 0 = 1

(showCReal 50 (almkvistGiullera n))

⚫

powInteger b 1 = b

⚫

powInteger b e = b `timesInteger` powInteger b (e `minusInteger` 1)

⚫

printf "\nPi after %d iterations:\n" iterations

⚫

putStrLn $ showCReal 70 $ almkvistGiulleraPi iterations

-- The ~~almkvistGiulleraIntegral~~ part of the Nth term in the Almkvist-Giullera series

-- The integral part of the Nth term in the Almkvist-Giullera series

almkvistGiulleraIntegral :: Integer -> Integer

almkvistGiulleraIntegral n =

let polynomial = (532 `timesInteger` n `timesInteger` n) `plusInteger` (126 `timesInteger` n) `plusInteger` 9

Line 605:

Line 607:

-- The exponent for 10 in the Nth term of the series

tenExponent :: Integer -> Integer

tenExponent n = 3 `minusInteger` (6 `timesInteger` (1 `plusInteger` n))

-- The Nth term in the series

-- The Nth term in the series (integral * 10^tenExponent)

almkvistGiullera :: Integer -> CReal

almkvistGiullera n = fromInteger (almkvistGiulleraIntegral n) / fromInteger (powInteger 10 (abs (tenExponent n)))

-- The sum of the first N terms

almkvistGiulleraSum :: Integer -> CReal

almkvistGiulleraSum n = sum $ map almkvistGiullera [0 .. n]

-- The approximation of pi from the first N terms

almkvistGiulleraPi :: Integer -> CReal

almkvistGiulleraPi n = sqrt $ 1 / almkvistGiulleraSum n

⚫

-- Utility: factorial for arbitrary-precision integers

⚫