source code
--
-- The Computer Language Benchmarks Game
-- https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
--
-- Contributed by Olof Kraigher and Don Stewart.
-- GHC 7.8.1 fix by Ersin Er
--
-- 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
import System.IO.Unsafe
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)
notes, command-line, and program output
NOTES:
64-bit Ubuntu quad core
The Glorious Glasgow Haskell
Compilation System,
version 9.10.1
LLVM version 18.1.3
Thu, 01 Aug 2024 19:28:26 GMT
MAKE:
mv nbody.ghc nbody.hs
~/.ghcup/bin/ghc --make -fllvm -O2 -XBangPatterns -threaded -rtsopts nbody.hs -o nbody.ghc_run
Loaded package environment from /home/dunham/.ghc/x86_64-linux-9.10.1/environments/default
[1 of 2] Compiling Main ( nbody.hs, nbody.o )
nbody.hs:25:29: warning: [GHC-63394] [-Wx-partial]
In the use of ‘head’
(imported from Prelude, but defined in GHC.Internal.List):
"This is a partial function, it throws an error on empty lists. Use pattern matching, 'Data.List.uncons' or 'Data.Maybe.listToMaybe' instead. Consider refactoring to use "Data.List.NonEmpty"."
|
25 | n <- getArgs >>= readIO.head
| ^^^^
[2 of 2] Linking nbody.ghc_run
rm nbody.hs
17.77s to complete and log all make actions
COMMAND LINE:
./nbody.ghc_run +RTS -N4 -RTS 50000000
PROGRAM OUTPUT:
-0.169075164
-0.169059907