Shootout/Partial sums
< Shootout
Jump to navigation
Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
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
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
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
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