I've been dabbling with implementing the "n-body" benchmark over at the "Great Computer Language Shootout"... http://shootout.alioth.debian.org/benchmark.php?test=nbody ...and I've stumbled on to a problem with space leaks. The code below works great for low numbers of iterations, but crashes and burns with large numbers. You can temporarily fix it up by specifying larger stack and heap sizes (with the +RTS -K & -H options), but this only gets you so far. I haven't been able to glean much information from the profiling reports. You can see the reports I generated (like the biography report) over at http://sleepingsquirrel.org/nbody. I've tried several different variations on the "advance" function, as well as trying "seq" in numerous spots, but I haven't stumbled on to the right combination yet. Any pointers? Thanks, Greg Buchholz
-- n-body simulation -- compile: ghc -O2 -o nbody nbody.hs -- run: nbody 100000
import System(getArgs) import Numeric
data Vec = V !Double !Double !Double deriving Show instance Num Vec where (V a b c) + (V x y z) = (V (a+x) (b+y) (c+z)) (V a b c) - (V x y z) = (V (a-x) (b-y) (c-z)) fromInteger 0 = V 0.0 0.0 0.0 -- for sum function instance Eq Vec where (V a b c) == (V x y z) = (a==x) && (b==y) && (c==z)
dot (V a b c) (V x y z) = a*x + b*y + c*z scale (V a b c) n = V (n*a) (n*b) (n*c) mag (V x y z) = sqrt (x*x + y*y + z*z)
data Planet = Planet Vec Vec Double deriving Show --Position Velocity Mass dist (Planet p1 _ _) (Planet p2 _ _) = mag $ p1 - p2 mass (Planet _ _ m) = m vel (Planet _ v _) = v pos (Planet p _ _) = p
main = do [arg] <- getArgs let iter = read arg let n = 5::Int let bodies = offset_momentum n [sun, jupiter, saturn, neptune, uranus] let begin = energy n bodies putStrLn $ showFFloat (Just 9) begin "" let final = (iterate (advance n 0.01) bodies) let end = energy n (final !! iter) putStrLn $ showFFloat (Just 9) end ""
days_per_year = 365.24 solar_mass = 4.0 * pi * pi
advance :: Int -> Double -> [Planet] -> [Planet] advance n dt (p:ps) = take n ((Planet (pos p + delta_x) (new_v) (mass p) ) : advance n dt (ps++[p])) where delta_v = sum (map (\q -> (pos p - pos q) `scale` ((mass q)*dt/(dist p q)^3)) ps) new_v = (vel p) - delta_v delta_x = new_v `scale` dt
energy:: Int -> [Planet] -> Double energy n ps = kinetic - potential where kinetic = 0.5 * (sum (map (\q->(mass q)*((vel q) `dot` (vel q))) ps)) potential = sum [(mass (ps!!i))*(mass (ps!!j))/(dist (ps!!i) (ps!!j)) | i<-[0..n-1], j<-[i+1..n-1]]
offset_momentum n ((Planet p v m):ps) = (Planet p new_v m):ps where new_v = (sum (map (\n->(vel n) `scale` (mass n)) ps)) `scale` ((-1.0)/solar_mass)
{- rotations 0 _ = [] rotations n xss@(x:xs) = flipped:rotations (n-1) flipped where flipped = xs++[x]
advance :: Int -> Double -> [Planet] -> [Planet] advance n dt pss = map (adv dt) (rotations n pss)
adv dt (p:ps) = Planet (pos p + delta_x) (new_v) (mass p) where delta_v = sum (map (\q -> (pos p - pos q) `scale` ((mass q)*dt/(dist p q)^3)) ps) new_v = (vel p) - delta_v delta_x = new_v `scale` dt --} jupiter = (Planet (V 4.84143144246472090e+00 (-1.16032004402742839e+00) (-1.03622044471123109e-01)) (V ( 1.66007664274403694e-03 * days_per_year) ( 7.69901118419740425e-03 * days_per_year) ((-6.90460016972063023e-05) * days_per_year)) (9.54791938424326609e-04 * solar_mass))
saturn = (Planet (V 8.34336671824457987e+00 4.12479856412430479e+00 (-4.03523417114321381e-01)) (V (-2.76742510726862411e-03 * days_per_year) (4.99852801234917238e-03 * days_per_year) (2.30417297573763929e-05 * days_per_year)) (2.85885980666130812e-04 * solar_mass))
uranus = (Planet (V 1.28943695621391310e+01 (-1.51111514016986312e+01) (-2.23307578892655734e-01)) (V (2.96460137564761618e-03 * days_per_year) (2.37847173959480950e-03 * days_per_year) (-2.96589568540237556e-05 * days_per_year)) (4.36624404335156298e-05 * solar_mass))
neptune = (Planet (V 1.53796971148509165e+01 (-2.59193146099879641e+01) 1.79258772950371181e-01) (V (2.68067772490389322e-03 * days_per_year) (1.62824170038242295e-03 * days_per_year) (-9.51592254519715870e-05 * days_per_year)) (5.15138902046611451e-05 * solar_mass))
sun = (Planet (V 0.0 0.0 0.0) (V 0.0 0.0 0.0) solar_mass)