I thought it would be a nice exercise (and a good learning experience) to
try and solve the nbody problem from the debian language shootout.
Unfortunately, my code sucks. There is a massive space leak and performance
is even worse. On the bright side, my implementation is purely functional.
My twist: I manually unrolled a few loops from the C++ version.
I believe that most of my performance problems stem from my abuse of tuple.
The bodies are passed as a tuple of planets, a planet is a tuple of
(position, velocity, mass) and the vectors position and velocity are also
tuples of type double. My lame justification for that is to make it nice and
handy to pass data around.
Any tips would be greatly appreciated.
--ryan
{--
The Great Computer Language Shootout
http://shootout.alioth.debian.org/
N-body problem
C version contributed by Christoph Bauer
converted to C++ and modified by Paul Kitchin
This crappy Haskell version by Ryan Dickie based on the above two
--}
import System
import Text.Printf
(a,b,c) .+ (x,y,z) = (a+x,b+y,c+z)
(a,b,c) .- (x,y,z) = (a-x,b-y,c-z)
x .* (a,b,c) = (x*a,x*b,x*c)
mag2 (x,y,z) = x*x + y*y + z*z
mag = sqrt . mag2
--some getter functions
pVel (_,vel,_) = vel
pPos (pos,_,_) = pos
pMass (_,_,mass) = mass
days_per_year = 365.24::Double
solar_mass = 4 * pi * pi::Double
delta_time = 0.01::Double
planets = (p1,p2,p3,p4,p5) where
p1 = ((0,0,0),(0,0,0),solar_mass)
p2 = (p2pos,p2vel,p2mass)
p3 = (p3pos,p3vel,p3mass)
p4 = (p4pos,p4vel,p4mass)
p5 = (p5pos,p5vel,p5mass)
p2pos = (4.84143144246472090e+00,-1.16032004402742839e+00,-1.03622044471123109e-01)
p2vel = days_per_year .* (1.66007664274403694e-03, 7.69901118419740425e-03 ,-6.90460016972063023e-05)
p2mass = 9.54791938424326609e-04 * solar_mass
p3pos = (8.34336671824457987e+00,4.12479856412430479e+00,-4.03523417114321381e-01)
p3vel = days_per_year .* (-2.76742510726862411e-03, 4.99852801234917238e-03, 2.30417297573763929e-05)
p3mass = 2.85885980666130812e-04 * solar_mass
p4pos = (1.28943695621391310e+01,-1.51111514016986312e+01,-2.23307578892655734e-01)
p4vel = days_per_year .* (2.96460137564761618e-03,2.37847173959480950e-03,-2.96589568540237556e-05)
p4mass = 4.36624404335156298e-05 * solar_mass
p5pos = (1.53796971148509165e+01,-2.59193146099879641e+01,1.79258772950371181e-01)
p5vel = days_per_year .* (2.68067772490389322e-03,1.62824170038242295e-03,-9.51592254519715870e-05)
p5mass = 5.15138902046611451e-05 * solar_mass
offset_momentum (p1,p2,p3,p4,p5) = (pp1,p2,p3,p4,p5) where
pp1 = (pPos p1,ppvel,pMass p1)
ppvel = (-1.0 / solar_mass) .* momentum
momentum = (mul p2) .+ (mul p3) .+ (mul p4) .+ (mul p5)
mul (_,vel,mass) = (mass .* vel)
--trick here is to unroll each loop.. based on the c++ version
advance ps = update $! adv4 $ adv3 $ adv2 $ adv1 ps
update (p1,p2,p3,p4,p5) = (up p1,up p2,up p3,up p4,up p5) where
up (pos,vel,mass) = ((pos .+ (delta_time .* vel)),vel,mass)
adv4 (p1,p2,p3,p4,p5) = (p1,p2,p3,pp4,pp5) where
il45 = innerLoop p4 p5
pp4 = fst il45
pp5 = snd il45
adv3 (p1,p2,p3,p4,p5) = (p1,p2,pp3,pp4,pp5) where
il34 = innerLoop p3 p4
il35 = innerLoop (fst il34) p5
pp3 = fst il35
pp4 = snd il34
pp5 = snd il35
adv2 (p1,p2,p3,p4,p5) = (p1,pp2,pp3,pp4,pp5) where
il23 = innerLoop p2 p3
il24 = innerLoop (fst il23) p4
il25 = innerLoop (fst il24) p5
pp2 = fst il25
pp3 = snd il23
pp4 = snd il24
pp5 = snd il25
adv1 (p1,p2,p3,p4,p5) = (pp1,pp2,pp3,pp4,pp5) where
il12 = innerLoop p1 p2
il13 = innerLoop (fst il12) p3
il14 = innerLoop (fst il13) p4
il15 = innerLoop (fst il14) p5
pp1 = fst il15
pp2 = snd il12
pp3 = snd il13
pp4 = snd il14
pp5 = snd il15
innerLoop p1 p2 = (pp1,pp2) where
difference = (pPos p1) .- (pPos p2)
distance_squared = mag2 difference
distance = sqrt distance_squared
magnitude = delta_time / (distance * distance_squared)
planet2_mass_magnitude = (pMass p2) * magnitude
planet1_mass_magnitude = (pMass p1) * magnitude
pp1 = (pPos p1, (pVel p1) .- (planet2_mass_magnitude .* difference), pMass p1)
pp2 = (pPos p2, (pVel p2) .+ (planet1_mass_magnitude .* difference), pMass p2)
energy (p1,p2,p3,p4,p5) = sum2 where
sum2 = loop5 + loop4 + loop3 + loop2 + loop1
loop0 (pos,vel,mass) = (0.5) * mass * (mag2 vel)
loop5 = (loop0 p5)
loop4 = (loop0 p4) - (delE p4 p5)
loop3 = (loop0 p3) - (delE p3 p5) - (delE p3 p4)
loop2 = (loop0 p2) - (delE p2 p5) - (delE p2 p4) - (delE p2 p3)
loop1 = (loop0 p1) - (delE p1 p5) - (delE p1 p4) - (delE p1 p3) - (delE p1 p2)
delE (pos,_,mass) (pos2,_,mass2) = (mass * mass2) / (mag (pos .- pos2))
runIt 0 args = args
runIt cnt args = runIt (cnt-1) $! advance args
--runIt cnt args = advance $ runIt (cnt-1) args
main :: IO()
main = do
n <- getArgs >>= readIO.head
let ps = offset_momentum planets
let results = runIt n ps
printf "%.9f\n" (energy ps)
printf "%.9f\n" (energy results)
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
