Haskell implementation
Integration Scheme: 2nd order predictor corrector leapfrog
Compiler: g++4.6.1 with -O4 , ghc 7.0.3 with --make -O3
Operating System: Arch Linux
Hardware: Intel(R) Core(TM) i7 CPU 975 @ 3.33GHz
Initial Conditions: input/input128 and input/input256
Time step: constant and shared 1.e-3 N-body unit
tEnd = 1.0, file = input/input128
source tCPU(s) dE/E
main.cpp 0.303 1.7273e-6
MainUnboxed.hs 6.780 2.1276e-6
Main.hs 12.196 2.1276e-6
MainWithList.hs 55.336 2.1276e-6
tEnd = 1.0, file = input/input256
source tCPU(s) dE/E
main.cpp 1.206 -6.6877e-4
MainUnboxed.hs 31.508 -6.6877e-4
Main.hs 49.187 -6.6877e-4
MainWithList.hs 494.682 -6.6877e-4
Note:
Main.hs uses 'vector', an efficient array library of Haskell.
MainUnboxed.hs uses unboxed version of Vector.
MainWithList.hs uses Haskell's list, hope it's more intuitive.
You may need
> cabal install vector
to compile those codes.
import Control.Monad (when)
import Data.List (transpose)
import Data.Vector (Vector, (!))
import qualified Data.Vector as V
import Prelude hiding (Real)
type Real = Double
type Cluster = Vector Star
dimension = 3
zero = V.replicate dimension 0
dt = 1e-3
innerProduct xs ys = V.sum $ V.zipWith (*) xs ys
data Star
= Star
{ mass :: Real,
posi :: Vector Real,
velo :: Vector Real
}
deriving (Eq, Show)
energies stars = [total, kinetic, potential] where
kinetic = V.sum $ V.map kin stars
kin (Star m _ v) = 0.5 * m * innerProduct v v
potential =
V.sum $
flip V.imap stars $
\i si ->
V.sum $
flip V.imap (V.take i stars) $
\j sj -> pot si sj
pot si sj =
let rij = V.zipWith (-) (posi si) (posi sj)
in -mass si * mass sj / sqrt(innerProduct rij rij)
total = kinetic + potential
update :: Cluster -> Cluster
update stars = stars'
where
stars' = V.zipWith3 Star ms rs' vs'
ms = V.map mass stars
rs = V.map posi stars
vs = V.map velo stars
as = calcAcc ms rs
rs' = flip V.imap rs $
\i r -> V.generate dimension $
\ex -> rs!i!ex + dt * vs!i!ex + 0.5 * as!i!ex * dt*dt
as' = calcAcc ms rs'
vs' = flip V.imap vs $
\i v -> V.generate dimension $
\ex -> vs!i!ex + 0.5*dt *(as!i!ex + as'!i!ex)
calcAcc ms rs = let n = V.length ms in V.generate n $
\i -> V.generate dimension $
\ex -> V.sum $ V.generate n $
\j -> if i==j
then 0
else let
rij = V.zipWith (-) (rs!i) (rs!j)
r2 = innerProduct rij rij
apre = 1/sqrt(r2*r2*r2)
in - ms!j * apre * rij!ex
main = do
str <- getContents
let
input :: Vector (Vector Double)
input = V.fromList $ map V.fromList $
map (map read) $
filter ((>=8) . length) $
map words $ lines str
stars = V.generate (V.length input) $
\i -> Star
{ mass = input!i!1,
posi = V.fromList [input!i!j | j <- [2,3,4]],
velo = V.fromList [input!i!j | j <- [5,6,7]]
}
e0 = head $ energies stars
loop 0 0 e0 stars
loop t k e0 stars = do
let stars' = update stars
es = energies stars'
dE = (head es - e0) / e0
when (mod k 10 == 0) $ putStrLn $
"t= " ++ show t ++ " E= " ++ unwords (map show es) ++ " dE = " ++ show dE
when (k<1000) $ loop (t+dt) (k+1) e0 stars'
import Control.Monad (when)
import Data.List (transpose)
import Prelude hiding (Real)
type Real = Double
type Cluster = [Star]
dimension = 3
dt = 1e-3
innerProduct xs ys = sum $ zipWith (*) xs ys
mapi xs f = map (uncurry f) $ zip [0..] xs
zip3With xs ys zs f = zipWith3 f xs ys zs
generate n f = map f [0..n-1]
data Star
= Star
{ mass :: Real,
posi :: [Real],
velo :: [Real]
}
deriving (Eq, Show)
energies :: [Star] -> [Double]
energies stars = [total, kinetic, potential] where
kinetic = sum $ map kin stars
kin (Star m _ v) = 0.5 * m * innerProduct v v
potential =
sum $ mapi stars $
\i si ->
sum $
mapi (take i stars) $
\j sj -> pot si sj
pot si sj =
let rij = zipWith (-) (posi si) (posi sj)
in -mass si * mass sj / sqrt(innerProduct rij rij)
total = kinetic + potential
update :: Cluster -> Cluster
update stars = stars'
where
stars' = zipWith3 Star ms rs' vs'
ms = map mass stars
rs = map posi stars
vs = map velo stars
as = calcAcc ms rs
rs' = zip3With rs vs as $
\r v a -> generate dimension $
\ex -> r!!ex + dt * v!!ex + 0.5 * a!!ex * dt*dt
as' = calcAcc ms rs'
vs' = zip3With vs as as' $
\v a a' -> generate dimension $
\ex -> v!!ex + 0.5*dt *(a!!ex + a'!!ex)
calcAcc ms rs = let n = length ms in generate n $
\i -> generate dimension $
\ex -> sum $ generate n $
\j -> if i==j
then 0
else let
rij = zipWith (-) (rs!!i) (rs!!j)
r2 = innerProduct rij rij
apre = 1/sqrt(r2*r2*r2)
in - ms!!j * apre * rij!!ex
main = do
str <- getContents
let
input :: [[Double]]
input = map (map read) $
filter ((>=8) . length) $
map words $
lines str
stars = flip map input $
\xs -> Star
{ mass = xs!!1,
posi = [xs!!j | j <- [2,3,4]],
velo = [xs!!j | j <- [5,6,7]]
}
e0 = head $ energies stars
loop 0 0 e0 stars
loop t k e0 stars = do
let stars' = update stars
es = energies stars'
dE = (head es - e0) / e0
when (mod k 10 == 0) $ putStrLn $
"t= " ++ show t ++ " E= " ++ unwords (map show es) ++ " dE = " ++ show dE
when (k<1000) $ loop (t+dt) (k+1) e0 stars'
{-# LANGUAGE FlexibleInstances #-}
import Control.Monad (when)
import Data.List (transpose)
import Data.Vector.Unboxed (Vector, (!))
import qualified Data.Vector.Unboxed as V
import Prelude hiding (Real)
type Real = Double
type R3 = (Real, Real, Real)
dimension = 3
dt = 1e-3
instance Num a => Num (a,a,a) where
(a,b,c) + (d,e,f) = (a+d, b+e, c+f)
(a,b,c) - (d,e,f) = (a-d, b-e, c-f)
negate (d,e,f) = (-d, -e, -f)
fromInteger x = (fromInteger x, fromInteger x, fromInteger x)
innerProduct (a,b,c) (d,e,f) = a*d + b*e + c*f
scalarProduct a (d,e,f) = (a*d, a*e, a*f)
data Cluster
= Cluster
{ size :: Int,
mass :: Vector Real,
posi :: Vector R3,
velo :: Vector R3
}
deriving (Eq, Show)
calcAcc :: Cluster -> Vector R3
calcAcc (Cluster n ms rs _) =
V.generate n $
\i -> sum $ [force i j |j <- [0..n-1], i/=j]
where
force i j =
let rji = rs!j - rs!i
in (ms!j / sqrt(innerProduct rji rji^(3::Int))) `scalarProduct` rji
update :: Cluster -> Cluster
update stars@(Cluster n ms rs vs) = Cluster n ms rs' vs'
where
as = calcAcc stars
as' = calcAcc stars { posi = rs' }
rs' = flip V.imap rs $
\i r -> r +
scalarProduct dt (vs!i) +
scalarProduct (0.5*dt*dt) (as!i)
vs' = flip V.imap vs $
\i v -> v +
scalarProduct (0.5*dt) (as!i + as'!i)
energies :: Cluster -> [Real]
energies (Cluster n ms rs vs) = [kinetic+potential, kinetic, potential]
where
kinetic =
sum $
[0.5 * ms!i * innerProduct (vs!i) (vs!i)| i<-[0..n-1]]
potential =
sum $
[ let rij = rs!i - rs!j in
-ms!i * ms!j / sqrt (innerProduct rij rij)
| i <- [0..n-2], j <- [i+1..n-1] ]
main = do
str <- getContents
let
input :: [[Double]]
input = map (map read) $
filter ((>=8) . length) $
map words $ lines str
stars
= Cluster
{ size = length input,
mass = V.fromList $ map (!!1) input,
posi = V.fromList $
map (\l -> (l!!2, l!!3, l!!4)) input,
velo = V.fromList $
map (\l -> (l!!5, l!!6, l!!7)) input
}
e0 = head $ energies stars
loop 0 0 e0 stars
loop t k e0 stars = do
let stars' = update stars
es = energies stars'
dE = (head es - e0) / e0
when (mod k 10 == 0) $ putStrLn $
"t= " ++ show t ++ " E= " ++ unwords (map show es) ++ " dE = " ++ show dE
when (k<1000) $ loop (t+dt) (k+1) e0 stars'