Shootout/Partial sums
From HaskellWiki
< Shootout(Difference between revisions)
(→Benchmarks) |
|||
| Line 21: | Line 21: | ||
0.693127181 Alternating Harmonic | 0.693127181 Alternating Harmonic | ||
0.785388163 Gregory | 0.785388163 Gregory | ||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
| - | |||
== Proposed entry == | == Proposed entry == | ||
Current revision
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 nave 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
Contents |
1 Proposed entry
Int comparision, cache 2/3 on the stack
{-# OPTIONS -O2 -fexcess-precision #-} -- -- The Great Computer Language Shootout -- http://shootout.alioth.debian.org/ -- -- Translation of the Clean by Don Stewart -- import System import Numeric import Monad main = do n <- getArgs >>= readIO . head mapM_ draw $ sigma (n::Int) (1::Int) 1 (2/3) 0 0 0 0 0 0 0 0 0 draw (s,t) = putStrLn $ (showFFloat (Just 9) (s::Double) []) ++ "\t" ++ t sigma !n !i !alt !tt !a1 !a2 !a3 !a4 !a5 !a6 !a7 !a8 !a9 | i <= n = sigma n (i+1) (-alt) tt (a1 + tt ** (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
2 Old 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
3 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
4 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
5 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
6 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))
7 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
7.1 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
