# Shootout/Partial sums

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

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

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

main = do
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

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

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

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

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

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

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