Shootout/Partial sums
This is a Shootout Entry for partial-sums.
The description of the program:
diff program output N = 25000 with this output file to check your program is correct before contributing.
Each program should use the same na�ve iterative double-precision algorithms to calculate partial-sums of the series.
For more information see "General Series." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/topics/GeneralSeries.html
Correct output
3.000000000 (2/3)^k 314.770573775 k^-0.5 0.999960002 1/k(k+1) 30.314520404 Flint Hills 42.994899946 Cookson Hills 10.703866769 Harmonic 1.644894068 Riemann Zeta 0.693127181 Alternating Harmonic 0.785388163 Gregory
= Benchmarks
Timing on Linux/P4, N=2500000
||Author|| Time || Flags|| ||dons-fast3|| 2.207 || -O2 -optc-O3 -optc-ffast-math -fexcess-precision || ||dons-fast2|| 2.578 || -O2 -optc-O3 -optc-ffast-math -fexcess-precision || ||dons-fast|| 2.746 || -O2 -optc-O3 -optc-ffast-math -fexcess-precision || ||dons || 3.790 || -O2 -optc-O3 -optc-ffast-math -fexcess-precision || ||dons || 3.817 || -O2 -optc-O3 -fexcess-precision || ||dons || 3.906 || -O2 -fasm -fexcess-precision || ||chris || 3.998 || -O2 -optc-O3 -optc-ffast-math -fexcess-precision || ||chris || 4.043 || -O2 -optc-O3 -fexcess-precision || ||dons || 4.140 || -O2 -optc-O3 -optc-ffast-math || ||chris || 4.371 || -O2 -fasm -fexcess-precision || ||chris || 4.660 || -O2 -optc-O3 -optc-ffast-math || ||dons (eye-candy) || 5.098 || -O2 -optc-O3 -optc-ffast-math -fexcess-precision ||
I agree with this chart (running on powerboook G4) -- ChrisKuklewicz
Proposed entry
Ported to GHC 6.6. Performance looks good. Submitted.
Compile with ghc -O -fglasgow-exts -fbang-patterns -optc-O3 -optc-march=pentium4.
{-# OPTIONS -fexcess-precision #-}
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Translation of the Clean by Don Stewart
--
-- Compile with
-- ghc -O -fglasgow-exts -fbang-patterns -optc-O3 -optc-march=pentium4
--
import System
import Numeric
import Monad
main = do
n <- getArgs >>= readIO . head
mapM_ draw $ go n (1::Int) 1 0 0 0 0 0 0 0 0 0
draw (s,t) = putStrLn $ (showFFloat (Just 9) (s::Double) []) ++ "\t" ++ t
go !n !i !alt !a1 !a2 !a3 !a4 !a5 !a6 !a7 !a8 !a9
| k <= n = go n (i+1) (-alt)
(a1 + (2/3) ** (k-1))
(a2 + sqrt dk)
(a3 + 1 / (k * (k + 1)))
(a4 + 1 / (k3 * sk * sk))
(a5 + 1 / (k3 * ck * ck))
(a6 + dk)
(a7 + dk * dk)
(a8 + alt * dk)
(a9 + alt / (2 * k - 1))
| otherwise = [(a1, "(2/3)^k"), (a2, "k^-0.5"), (a3, "1/k(k+1)")
,(a4, "Flint Hills"), (a5, "Cookson Hills"), (a6, "Harmonic")
,(a7, "Riemann Zeta"),(a8, "Alternating Harmonic"), (a9, "Gregory")]
where k = fromIntegral i
k3 = k2*k; k2 = k*k; dk = 1/k; sk = sin k; ck = cos k
Old entry
Use {{{1 / sqrt k}}} for {{{k ** (-0.5)}}} as the other entries do.
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Haskell version of Isaac Gouy's Clean version, translated by Don Stewart
--
import System; import Numeric
main = do n <- getArgs >>= readIO . head
let sums = loop (1::Int) n 1 0 0 0 0 0 0 0 0 0
fn (s,t) = putStrLn $ (showFFloat (Just 9) s []) ++ "\t" ++ t
mapM_ (fn :: (Double, String) -> IO ()) (zip sums names)
names = ["(2/3)^k", "k^-0.5", "1/k(k+1)", "Flint Hills", "Cookson Hills"
, "Harmonic", "Riemann Zeta", "Alternating Harmonic", "Gregory"]
loop i n alt a1 a2 a3 a4 a5 a6 a7 a8 a9
| i !n !alt !a1 !a2 !a3 !a4 !a5 !a6 !a7 !a8 !a9 !False = undefined -- strict
| k > n = [ a1, a2, a3, a4, a5, a6, a7, a8, a9 ]
| otherwise = loop (i+1) n (-alt)
(a1 + (2/3) ** (k-1))
(a2 + 1 / sqrt k)
(a3 + 1 / (k * (k + 1)))
(a4 + 1 / (k3 * sk * sk))
(a5 + 1 / (k3 * ck * ck))
(a6 + dk)
(a7 + 1 / k2)
(a8 + alt * dk)
(a9 + alt / (2 * k - 1))
where k3 = k2*k; k2 = k*k; dk = 1/k; k = fromIntegral i; sk = sin k; ck = cos k; x!y = x`seq`y
Faster still
Takes an idea to pass an int counter to the loop from the C entry. I walked the C output generated by GHC, and the same operations are generated (though the current C entry rearranges the computations in a more optimal order -- at least I think that's what they're doing).
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Haskell version of Isaac Gouy's Clean version, translated by Don Stewart
--
import System; import Numeric
main = do n <- getArgs >>= readIO . head
let sums = loop (1::Int) n 1 0 0 0 0 0 0 0 0 0
fn (s,t) = putStrLn $ (showFFloat (Just 9) s []) ++ "\t" ++ t
mapM_ (fn :: (Double, String) -> IO ()) (zip sums names)
names = ["(2/3)^k", "k^-0.5", "1/k(k+1)", "Flint Hills", "Cookson Hills"
, "Harmonic", "Riemann Zeta", "Alternating Harmonic", "Gregory"]
loop i n alt a1 a2 a3 a4 a5 a6 a7 a8 a9
| i !n !alt !a1 !a2 !a3 !a4 !a5 !a6 !a7 !a8 !a9 !False = undefined -- strict
| k > n = [ a1, a2, a3, a4, a5, a6, a7, a8, a9 ]
| otherwise = loop (i+1) n (-alt)
(a1 + (2/3) ** (k-1))
(a2 + k ** (-0.5))
(a3 + 1 / (k * (k + 1)))
(a4 + 1 / (k3 * sk * sk))
(a5 + 1 / (k3 * ck * ck))
(a6 + dk)
(a7 + 1 / k2)
(a8 + alt * dk)
(a9 + alt / (2 * k - 1))
where k3 = k2*k; k2 = k*k; dk = 1/k; k = fromIntegral i; sk = sin k; ck = cos k; x!y = x`seq`y
Current Entry
Cleverer idea, from Isaac Gouy. -O2 -optc-O3 -fexcess-precision -optc-ffast-math
Once I saw this, I checked the other entries, and the winning C gcc #4 entry uses combined loops as well. So they will accept code like this. -- ChrisKuklewicz
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Haskell version of Isaac Gouy's Clean version, translated by Don Stewart
--
import System; import Numeric
main = do n <- getArgs >>= readIO . head
let sums = loop 1 n 1 0 0 0 0 0 0 0 0 0
fn (s,t) = putStrLn $ (showFFloat (Just 9) s []) ++ "\t" ++ t
mapM_ (fn :: (Double, String) -> IO ()) (zip sums names)
names = ["(2/3)^k", "k^-0.5", "1/k(k+1)", "Flint Hills", "Cookson Hills"
, "Harmonic", "Riemann Zeta", "Alternating Harmonic", "Gregory"]
loop k n alt a1 a2 a3 a4 a5 a6 a7 a8 a9
| k !n !alt !a1 !a2 !a3 !a4 !a5 !a6 !a7 !a8 !a9 !False = undefined -- strict
| k > n = [ a1, a2, a3, a4, a5, a6, a7, a8, a9 ]
| otherwise = loop (k+1) n (-alt)
(a1 + (2/3) ** (k-1))
(a2 + k ** (-0.5))
(a3 + 1 / (k * (k + 1)))
(a4 + 1 / (k3 * sk * sk))
(a5 + 1 / (k3 * ck * ck))
(a6 + 1 / k)
(a7 + 1 / k2)
(a8 + alt / k)
(a9 + alt / (2 * k - 1))
where k3 = k2*k; k2 = k*k; sk = sin k; ck = cos k; x ! y = x `seq` y
Current Entry
Combines some ideas from Chris', with careful attention to how loops are unboxed. Compile with: -O2 -optc-O3 -optc-ffast-math -fexcess-precision
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Contributed by Chris Kuklewicz and Don Stewart
--
import System; import Numeric
main = do n <- getArgs >>= readIO . head
let run :: (String, Int -> Double -> Double) -> IO ()
run (s,f) = putStrLn $ (showFFloat (Just 9) (f n 0) []) ++ "\t" ++ s
mapM_ run [("(2/3)^k", twth), ("k^-0.5" , k05)
,("1/k(k+1)", kk1), ("Flint Hills", flhl)
,("Cookson Hills",cook), ("Harmonic", harm)
,("Riemann Zeta", rmzt), ("Alternating Harmonic",alth)
,("Gregory", greg)]
twth k s = if k == -1 then s else twth (k-1) (s+((2/3)**fromIntegral k))
k05 k s = if k == 0 then s else k05 (k-1) (s+(fromIntegral k**(-0.5)))
kk1 k s = if k == 0 then s else let j = fromIntegral k in kk1 (k-1) (s+1/(j*(j+1)))
harm k s = if k == 0 then s else harm (k-1) (s+1/fromIntegral k)
rmzt k s = if k == 0 then s else rmzt (k-1) (s+1/(fromIntegral k**2))
flhl k s = if k == 0 then s else let j = fromIntegral k in flhl (k-1) (s+1/((j**3)*(sin j**2)))
cook n s = loop 1 s
where loop k s = if k == n then s
else let j = fromIntegral k in loop (k+1) (s+1/((j**3)*((cos j)**2)))
alth k s = loop k (-1) s
where loop k a s = if k == 0 then s else loop (k-1) (-a) (s+a/fromIntegral k)
greg k s = loop k (-1) s
where loop k a s = if k == 0 then s else loop (k-1) (-a) (s+a/(2*fromIntegral k-1))
Chris' Entry
This seems to run very very fast for me.
--
-- The Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- By Chris Kuklewicz
-- adpated from Isaac Gouy's Clean entry
--
-- compile with: ghc -O3 -optc-O3 -optc-ffast-math -fexcess-precision partial-sums.hs -o partial-sums.ghc_run
-- run with: ./partial-sums.ghc_run %A
--
import Control.Monad(mapM_)
import System(getArgs)
import Text.Printf
todo = [("(2/3)^k",twoThirds),("k^-0.5",k05),("1/k(k+1)",kk1),
("Flint Hills",flintHills),("Cookson Hills",cooksonHills),
("Harmonic",harmonic),("Riemann Zeta",riemannZeta),
("Alternating Harmonic",altHarmonic),("Gregory",gregory)]
run n = mapM_ (\(name,fun)->printf "%.9f\t%s\n" (fun n) name) todo
main = do args <- getArgs
run $ if null args then (25000::Double) else read (head args)
twoThirds = succ . compute (\k-> (2/3)**k )
k05 = compute (\k-> k**(-0.5) )
kk1 = compute (\k-> 1/(k*(k+1)) )
flintHills = compute (\k-> 1/((k**3)*((sin k)**2)) )
cooksonHills n = cooksonHills' 1 n 0
where cooksonHills' k n sum | k == n = sum
| otherwise = cooksonHills'(k+1) n $! (sum + 1/((k**3)*((cos k)**2)))
harmonic = compute (\k-> 1/k )
riemannZeta = compute (\k-> 1/k**2 )
altHarmonic = altCompute (\k-> 1/k )
gregory = altCompute (\k -> 1/(2*k-1) )
compute term k = compute' k 0
where compute' k sum | k > 0 = compute' (k-1) $! (sum + term k)
| otherwise = sum
altCompute term k = compute' k (-1) 0
where compute' k s sum | k > 0 = compute' (k-1) (-s) $! (sum + s * (term k))
| otherwise = sum
An eye candy idea
An idea for a `most beautiful' point-free entry.
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Contributed by Don Stewart
--
import System
import Text.Printf
twoThirds = ((2/3) **)
k05 = (** (-0.5))
harmonic = (1 /)
riemannZeta = (1 /) . (** 2)
kk1 = (1 /) . (\k -> (k + 1) * k)
flintHills = (1 /) . (\k -> (k ** 3) * (sin k ** 2))
cooksonHills = (1 /) . (\k -> (k ** 3) * (cos k ** 2))
altHarmonic a = (a /)
gregory a = (a /) . (subtract 1 . (2 *))
main = do n <- getArgs >>= readIO . head
let run (s,f) = printf "%.9f\t%s\n" (compute f n :: Double) s
mapM_ run [("(2/3)^k", C twoThirds)
,("k^-0.5" , A k05)
,("1/k(k+1)", A kk1)
,("Flint Hills", A flintHills)
,("Cookson Hills", D cooksonHills)
,("Harmonic", A harmonic)
,("Riemann Zeta", A riemannZeta)
,("Alternating Harmonic", B altHarmonic)
,("Gregory", B gregory)]
data A = A (Double -> Double)
| C (Double -> Double)
| D (Double -> Double)
| B (Double -> Double -> Double)
-- todo, generalise all these loops:
compute (A f) k = loop k 0
where loop k s | k == 0 = s
| otherwise = loop (k-1) $! s + f k
compute (C f) k = loop k 0
where loop k s | k == -1 = s
| otherwise = loop (k-1) $! s + f k
compute (D f) n = loop 1 0
where loop k s | k == n = s
| otherwise = loop (k+1) $! s + f k
compute (B f) k = loop (-1) k 0
where loop a k s | k == 0 = s
| otherwise = loop (-a) (k-1) $! s + f a k