Personal tools

Shootout/Partial sums

From HaskellWiki

< Shootout(Difference between revisions)
Jump to: navigation, search
(Benchmarks)
Line 49: Line 49:
   
 
== Proposed entry ==
 
== Proposed entry ==
  +
  +
Int comparision, cache 2/3 on the stack
  +
  +
<haskell>
  +
{-# 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
  +
</haskell>
  +
  +
== Old entry ==
   
 
Ported to GHC 6.6. Performance looks good.
 
Ported to GHC 6.6. Performance looks good.

Revision as of 04:01, 17 July 2007

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 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

1.1 Todo

This entry could be improved now, that the rules have changed, to allow separate loops and certain caching of values. See the C++ example.

1.2 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

1.3 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

1.4 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

2 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

3 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

4 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))

5 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

5.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