Almkvist-Giullera formula for pi: Difference between revisions
Content added Content deleted
(→{{header|Haskell}}: Add implementation) |
|||
Line 571: | Line 571: | ||
3.1415926535897932384626433832795028841971693993751058209749445923078164 |
3.1415926535897932384626433832795028841971693993751058209749445923078164 |
||
</pre> |
</pre> |
||
=={{header|Haskell}}== |
|||
{{libheader|numbers}} |
|||
{{trans|Common Lisp}} |
|||
<lang haskell>import GHC.Integer |
|||
import Data.Number.CReal |
|||
import Text.Printf |
|||
-- Factorial for arbitrary-precision integers |
|||
facInteger :: Integer -> Integer |
|||
facInteger n = if n `leInteger` 1 then 1 else n `timesInteger` facInteger (n `minusInteger` 1) |
|||
-- Exponentiation for arbitrary-precision integers |
|||
powInteger :: Integer -> Integer -> Integer |
|||
powInteger 1 _ = 1 |
|||
powInteger _ 0 = 1 |
|||
powInteger b 1 = b |
|||
powInteger b e = b `timesInteger` powInteger b (e `minusInteger` 1) |
|||
-- The integral part of the Nth term in the Almkvist-Giullera series |
|||
integral :: Integer -> Integer |
|||
integral n = |
|||
let polynomial = (532 `timesInteger` n `timesInteger` n) `plusInteger` (126 `timesInteger` n) `plusInteger` 9 |
|||
numerator = 32 `timesInteger` (facInteger (6 `timesInteger` n)) `timesInteger` polynomial |
|||
denominator = 3 `timesInteger` (powInteger (facInteger n) 6) |
|||
in numerator `divInteger` denominator |
|||
-- 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 |
|||
almkvist'giullera :: Integer -> CReal |
|||
almkvist'giullera n = fromInteger (integral n) / fromInteger (powInteger 10 (abs (tenExponent n))) |
|||
-- The sum of the first N terms |
|||
a'g_sum :: Integer -> CReal |
|||
a'g_sum n = sum $ map almkvist'giullera [0 .. n] |
|||
-- The approximation of pi from the first N terms |
|||
a'g_pi :: Integer -> CReal |
|||
a'g_pi n = sqrt $ 1 / a'g_sum n |
|||
iterations = 52 |
|||
main :: IO () |
|||
main = do |
|||
(printf "N. %44s %4s %s\n" "Integral part of Nth term" "×10^" "=Actual value of Nth term") |
|||
mapM_ (\n -> printf "%d. %44d %4d %s\n" n (integral n) (tenExponent n) (showCReal 50 (almkvist'giullera n))) [0..9] |
|||
printf "\nPi after %d iterations:\n" iterations |
|||
putStrLn $ showCReal 70 $ a'g_pi iterations</lang> |
|||
{{Out}} |
|||
<pre>N. Integral part of Nth term ×10^ =Actual value of Nth term |
|||
0. 96 -3 0.096 |
|||
1. 5122560 -9 0.00512256 |
|||
2. 190722470400 -15 0.0001907224704 |
|||
3. 7574824857600000 -21 0.0000075748248576 |
|||
4. 312546150372456000000 -27 0.000000312546150372456 |
|||
5. 13207874703225491420651520 -33 0.00000001320787470322549142065152 |
|||
6. 567273919793089083292259942400 -39 0.0000000005672739197930890832922599424 |
|||
7. 24650600248172987140112763715584000 -45 0.000000000024650600248172987140112763715584 |
|||
8. 1080657854354639453670407474439566400000 -51 0.0000000000010806578543546394536704074744395664 |
|||
9. 47701779391594966287470570490839978880000000 -57 0.00000000000004770177939159496628747057049083997888 |
|||
Pi after 52 iterations: |
|||
3.1415926535897932384626433832795028841971693993751058209749445923078164</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |