Add the n-body shootout benchmark

.gitignore

@@ -53,6 +53,7 @@ real/veritas/veritas
 
 shootout/binary-trees/binary-trees
 shootout/fannkuch-redux/fannkuch-redux
+shootout/n-body/n-body
 shootout/pidigits/pidigits
 shootout/spectral-norm/spectral-norm
 

shootout/Makefile

 TOP = ..
 include $(TOP)/mk/boilerplate.mk
 
-SUBDIRS = binary-trees fannkuch-redux pidigits spectral-norm
+SUBDIRS = binary-trees fannkuch-redux n-body pidigits spectral-norm
 
 include $(TOP)/mk/target.mk
 

shootout/n-body/Main.hs

+--
+-- The Computer Language Benchmarks Game
+-- http://benchmarksgame.alioth.debian.org/
+--
+-- Contributed by Olof Kraigher and Don Stewart.
+--
+-- To be compiled with:
+--
+--  -O2 -fglasgow-exts -funbox-strict-fields -fbang-patterns -optc-O 
+--
+-- Don't enable -optc-mfpmath=sse -optc-msse2, this triggers a gcc bug on x86
+--
+
+import Foreign
+import Foreign.Storable
+import Foreign.Marshal.Alloc
+import Data.IORef
+import Control.Monad
+import System.Environment
+import Text.Printf
+
+main = do
+    n <- getArgs >>= readIO.head
+    initialize
+    offset_momentum
+    energy 0 planets >>= printf "%.9f\n"
+    replicateM_ n (advance planets)
+    energy 0 planets >>= printf "%.9f\n"
+
+offset_momentum = do
+    m <- foldr (.+.) (Vec 0 0 0)
+             `fmap` (mapM momentum
+                   . take (nbodies - 1)
+                   . iterate next $ next planets)
+
+    setVec (vel planets) $ (-1/solar_mass) *. m
+  where
+    momentum !p = liftM2 (*.) (mass p) (getVec (vel p))
+
+energy :: Double -> Ptr Double -> IO Double
+energy !e !p
+    | p == end = return e
+    | otherwise      = do
+        p1 <- getVec (pos p)
+        v1 <- getVec (vel p)
+        m1 <- mass p
+        e  <- energy2 p1 m1 e p2
+        energy (e + 0.5 * m1 * magnitude2 v1) p2
+    where p2 = next p
+
+energy2 !p1 !m1 !e !p
+    | p  == end = return e
+    | otherwise = do
+        p2 <- getVec (pos p)
+        v2 <- getVec (vel p)
+        m2 <- mass p
+        let distance = sqrt . magnitude2 $ p1 .-. p2
+        energy2 p1 m1 (e - m1 * m2 / distance) (next p)
+
+advance :: Ptr Double -> IO ()
+advance !p1 = when (p1 /= end) $ do
+    pos1 <- getVec $ pos p1
+    m1   <- mass p1
+    let go !p2
+            | p2 /= end = do
+                pos2 <- getVec (pos p2)
+                m2   <- mass p2
+                let vel2       = vel p2
+                    difference = pos1 .-. pos2
+                    distance2  = magnitude2 difference
+                    distance   = sqrt distance2
+                    magnitude  = delta_t / (distance2 * distance)
+                    mass_magn  = magnitude *. difference
+                vel1 -= m2 *. mass_magn
+                vel2 += m1 *. mass_magn
+                go (next p2)
+
+            | otherwise = do
+                v1 <- getVec vel1
+                p1 += delta_t *. v1
+    go p2
+    advance  p2
+  where
+    vel1 = vel p1
+    p2   = next p1
+
+------------------------------------------------------------------------
+
+planets :: Ptr Double
+planets = unsafePerformIO $ mallocBytes (7 * nbodies * 8) -- sizeOf double = 8
+
+nbodies :: Int
+nbodies = 5
+
+solar_mass, delta_t, days_per_year :: Double
+days_per_year = 365.24
+solar_mass    = 4 * pi ** 2;
+delta_t       = 0.01
+
+initialize = mapM_ newPlanet planets
+  where
+   dp = days_per_year
+   planets =
+    [0, 0, 0,
+     0, 0, 0,
+     1 * solar_mass,
+     4.84143144246472090e+00,        (-1.16032004402742839e+00), (-1.03622044471123109e-01),
+     1.66007664274403694e-03*dp,     7.69901118419740425e-03*dp, (-6.90460016972063023e-05)*dp,
+     9.54791938424326609e-04 * solar_mass,
+
+     8.34336671824457987e+00,        4.12479856412430479e+00,    (-4.03523417114321381e-01),
+     (-2.76742510726862411e-03)*dp,  4.99852801234917238e-03*dp, 2.30417297573763929e-05*dp,
+     2.85885980666130812e-04 * solar_mass,
+
+     1.28943695621391310e+01,        (-1.51111514016986312e+01), (-2.23307578892655734e-01),
+     2.96460137564761618e-03*dp,     2.37847173959480950e-03*dp, (-2.96589568540237556e-05)*dp,
+     4.36624404335156298e-05 * solar_mass,
+
+     1.53796971148509165e+01,        (-2.59193146099879641e+01), 1.79258772950371181e-01,
+     2.68067772490389322e-03*dp,     1.62824170038242295e-03*dp, (-9.51592254519715870e-05)*dp,
+     5.15138902046611451e-05 * solar_mass
+    ]
+
+------------------------------------------------------------------------
+-- Support for 3 dimensional mutable vectors
+
+data Vector3 = Vec !Double !Double !Double
+
+end :: Ptr Double
+end = inc planets $ nbodies * 7
+
+next  :: Ptr Double -> Ptr Double
+next p = inc p 7
+
+cursor :: IORef (Ptr Double)
+cursor = unsafePerformIO $ newIORef planets
+
+inc :: Ptr Double -> Int -> Ptr Double
+inc ptr n = plusPtr ptr (n * 8)
+
+newPlanet :: Double -> IO ()
+newPlanet !d = do
+    ptr <- readIORef cursor
+    pokeElemOff ptr 0 d
+    writeIORef cursor (inc ptr 1)
+
+pos :: Ptr Double -> Ptr Double
+pos ptr = ptr
+
+vel :: Ptr Double -> Ptr Double
+vel ptr = inc ptr 3
+
+mass :: Ptr Double -> IO Double
+mass ptr = peekElemOff ptr 6
+
+------------------------------------------------------------------------
+
+(Vec x y z) .+. (Vec u v w) = Vec (x+u) (y+v) (z+w)
+
+(Vec x y z) .-. (Vec u v w) = Vec (x-u) (y-v) (z-w)
+
+k *. (Vec x y z) = Vec (k*x) (k*y) (k*z) -- allocates
+
+magnitude2 (Vec x y z) = x*x + y*y + z*z
+
+------------------------------------------------------------------------
+
+getVec !p = liftM3 Vec (peek p) (f 1) (f 2)
+    where f = peekElemOff p
+
+setVec p (Vec x y z)= do
+    poke        p   x
+    pokeElemOff p 1 y
+    pokeElemOff p 2 z
+
+infix 4 +=
+infix 4 -=
+
+v1 += (Vec u v w) = do
+    x <- peek v1;          poke        v1   (x+u)
+    y <- peekElemOff v1 1; pokeElemOff v1 1 (y+v)
+    z <- peekElemOff v1 2; pokeElemOff v1 2 (z+w)
+
+v1 -= (Vec u v w)  = do
+    x <- peek v1;          poke        v1   (x-u)
+    y <- peekElemOff v1 1; pokeElemOff v1 1 (y-v)
+    z <- peekElemOff v1 2; pokeElemOff v1 2 (z-w)

shootout/n-body/Makefile

+TOP = ../..
+include $(TOP)/mk/boilerplate.mk
+include $(TOP)/mk/target.mk
+
+FAST_OPTS = 500000
+NORM_OPTS = 5000000
+SLOW_OPTS = 50000000  # official shootout setting
+
+# The benchmark game also uses -fllvm, which we can't since it might
+# not be available on the developer's machine.
+HC_OPTS += -O2 -XBangPatterns -fexcess-precision

shootout/n-body/n-body.faststdout

+-0.169075164
+-0.169096567

shootout/n-body/n-body.slowstdout

+-0.169075164
+-0.169059907

shootout/n-body/n-body.stdout

+-0.169075164
+-0.169083134