The Q6600
Benchmarks Game

n-body Haskell GHC #2 program

source code

-- The Computer Language Benchmarks Game
-- https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
--
-- Contributed by Branimir Maksimovic

import Foreign.Ptr
import Foreign.Storable
import Foreign.Marshal.Alloc
import Control.Monad
import System.Environment
import Text.Printf

main = do
    n <- getArgs >>= readIO.head :: IO Int
    pPlanets <- fromList planets
    nbodyInit pPlanets
    energy pPlanets >>= printf "%.9f\n" 
    run n pPlanets
    energy pPlanets >>= printf "%.9f\n" 

run 0 _ = return ()
run i p = do
        advance p
        run (i-1) p
        
data Planet = Planet { x,y,z,vx,vy,vz,mass :: !Double } deriving (Show)
    
offsetMomentum p (px,py,pz) = p {
                                   vx = -px / solar_mass,
                                   vy = -py / solar_mass,
                                   vz = -pz / solar_mass
                                }

nbodyInit pPlanets = do
    let init (px,py,pz) i = do
        if i < length planets
            then do
                p <- peekElemOff pPlanets i
                init (px + vx p * mass p,py + vy p * mass p, pz + vz p * mass p) (i+1)
            else return (px,py,pz)
    s <- init (0,0,0) 0
    p <- peek pPlanets
    poke pPlanets $ offsetMomentum p s

squared x y z = x * x + y * y + z * z
    
energy pPlanets = do
    let
        energy' e i = if i < length planets
                    then do
                        p <- peekElemOff pPlanets i
                        e1 <- energy'' p (i+1) e
                        e2 <- energy' e (i+1)
                        return $ e + 0.5 * mass p * squared (vx p) (vy p) (vz p)+e1+e2
                    else return e
        energy'' p j e = if j < length planets
                        then do
                            pj <- peekElemOff pPlanets j
                            let
                                distance = sqrt $ squared dx dy dz
                                dx = x pj - x p
                                dy = y pj - y p
                                dz = z pj - z p
                            e1 <- energy'' p (j+1) e
                            return $ e - (mass p * mass pj) / distance + e1
                        else return e
    energy' 0.0 0

advance pPlanets = do
    let 
        advance' i = 
            when (i < length planets) $ do
                    let
                        loop j = when (j < length planets) $ do
                                    ii <- peekElemOff pPlanets i
                                    jj <- peekElemOff pPlanets j
                                    let
                                        mag = dt / (dSquared * sqrt dSquared)
                                        dSquared = squared dx dy dz
                                        dx = x ii - x jj
                                        dy = y ii - y jj
                                        dz = z ii - z jj
                                    pokeV pPlanets i ii{
                                        vx = vx ii - dx * mass jj * mag,
                                        vy = vy ii - dy * mass jj * mag,
                                        vz = vz ii - dz * mass jj * mag
                                        }
                                    pokeV pPlanets j jj{
                                        vx = vx jj + dx * mass ii * mag,
                                        vy = vy jj + dy * mass ii * mag,
                                        vz = vz jj + dz * mass ii * mag
                                        }
                                    loop (j+1)
                    loop (i+1)
                    advance' (i+1)
        advance'' i = when (i < length planets) $ do
                            p <- peekElemOff pPlanets i
                            pokeC pPlanets i p { 
                                x = x p + dt * vx p,
                                y = y p + dt * vy p,
                                z = z p + dt * vz p
                                }
                            advance'' (i+1)
    advance' 0
    advance'' 0

planets = [sun, jupiter, saturn, uranus, neptune]
    
sun = Planet {x = 0, y = 0, z = 0,
              vx = 0, vy = 0, vz = 0,
              mass = solar_mass}
              
jupiter = Planet 
    {x = 4.84143144246472090e+00, y = -1.16032004402742839e+00, z= -1.03622044471123109e-01,
     vx = 1.66007664274403694e-03*dp, vy = 7.69901118419740425e-03*dp, vz = -6.90460016972063023e-05*dp,
     mass = 9.54791938424326609e-04 * solar_mass}

saturn = Planet
    { x = 8.34336671824457987e+00, y = 4.12479856412430479e+00, z = -4.03523417114321381e-01,
     vx = -2.76742510726862411e-03*dp,  vy = 4.99852801234917238e-03*dp, vz = 2.30417297573763929e-05*dp,
     mass = 2.85885980666130812e-04 * solar_mass}

uranus = Planet
    {x = 1.28943695621391310e+01,y = -1.51111514016986312e+01,z = -2.23307578892655734e-01,
     vx = 2.96460137564761618e-03*dp,vy = 2.37847173959480950e-03*dp, vz = -2.96589568540237556e-05*dp,
     mass = 4.36624404335156298e-05 * solar_mass}

neptune = Planet
    {x = 1.53796971148509165e+01,y = -2.59193146099879641e+01,z = 1.79258772950371181e-01,
     vx = 2.68067772490389322e-03*dp,vy = 1.62824170038242295e-03*dp, vz = -9.51592254519715870e-05*dp,
     mass = 5.15138902046611451e-05 * solar_mass}
     
days_per_year = 365.24
solar_mass    = 4 * pi ^ 2
dp = days_per_year
dt = 0.01

instance Storable Planet where
    sizeOf _ = 8 * dblSz
    alignment _ = dblSz
    peekElemOff p i = peek (plusPtr p (i * sizeOf (undefined::Planet)))
    pokeElemOff p i e = poke (plusPtr p (i * sizeOf e)) e
    peek p = do
        x <- peek (offset 0)
        y <- peek (offset 1)
        z <- peek (offset 2)
        vx <- peek (offset 3)
        vy <- peek (offset 4)
        vz <- peek (offset 5)
        mass <- peek (offset 6)
        return $ Planet {x=x,y=y,z=z,vx=vx,vy=vy,vz=vz,mass=mass}
            where
                offset i = plusPtr (castPtr p::Ptr Double) (i*8)
    poke p e = do
        poke (offset 0) $ x e
        poke (offset 1) $ y e
        poke (offset 2) $ z e
        poke (offset 3) $ vx e
        poke (offset 4) $ vy e
        poke (offset 5) $ vz e
        poke (offset 6) $ mass e
            where
                offset i = plusPtr (castPtr p::Ptr Double) (i*8)

dblSz = sizeOf (undefined::Double)

pokeC p i e = do
        poke (offset 0) $ x e
        poke (offset 1) $ y e
        poke (offset 2) $ z e
        where
            offset o = (plusPtr (castPtr p::Ptr Double)(o*8+i*64))

pokeV p i e = do
        poke (offset 3) $ vx e
        poke (offset 4) $ vy e
        poke (offset 5) $ vz e
        where
            offset o = (plusPtr (castPtr p::Ptr Double)(o*8+i*64))

fromList :: [Planet]->IO (Ptr Planet)
fromList l = do
    let len = length l
    pa <- mallocBytes (len * sizeOf (undefined::Planet))
    let 
        loop [] _ = return ()
        loop (x:xs) i = do
                            poke (pa `plusPtr` (i * sizeOf(undefined::Planet))) x
                            loop xs (i+1)
    loop l 0
    return pa
    

notes, command-line, and program output

NOTES:
64-bit Ubuntu quad core
The Glorious Glasgow Haskell Compilation System,
version 8.8.1


Tue, 05 May 2020 22:54:49 GMT

MAKE:
mv nbody.ghc-2.ghc nbody.ghc-2.hs
/opt/src/ghc-8.8.1/bin/ghc --make -fllvm -O2 -XBangPatterns -threaded -rtsopts  nbody.ghc-2.hs -o nbody.ghc-2.ghc_run
Loaded package environment from /home/dunham/.ghc/x86_64-linux-8.8.1/environments/default
[1 of 1] Compiling Main             ( nbody.ghc-2.hs, nbody.ghc-2.o )
You are using an unsupported version of LLVM!
Currently only 7 is supported.
We will try though...
Linking nbody.ghc-2.ghc_run ...
rm nbody.ghc-2.hs

23.98s to complete and log all make actions

COMMAND LINE:
./nbody.ghc-2.ghc_run +RTS -N4 -RTS 50000000

PROGRAM OUTPUT:
-0.169075164
-0.169059907