[add 'nbody' shootout entry to regression test, integrate Kleisli failing case to KindInference test
John Meacham <john@repetae.net>**20090224225206
 Ignore-this: fc0c2798f91d61813f089fed609365c0
] addfile ./regress/tests/3_shootout/nbody.args
addfile ./regress/tests/3_shootout/nbody.expected.stdout
addfile ./regress/tests/3_shootout/nbody.hs
hunk ./regress/tests/2_language/KindInference.hs 12
-main = return ()
+main = runKleisli (arr id) ()
hunk ./regress/tests/3_shootout/nbody.args 1
+129333
hunk ./regress/tests/3_shootout/nbody.expected.stdout 1
+-0.169075164
+-0.169058208
hunk ./regress/tests/3_shootout/nbody.hs 1
+{-# OPTIONS_GHC -fexcess-precision #-}
+--
+-- The Computer Language Shootout
+-- http://shootout.alioth.debian.org/
+--
+-- Contributed by Olof Kraigher, with help from Don Stewart.
+--
+-- Compile with:
+--
+--  -funbox-strict-fields -fglasgow-exts -fbang-patterns -O3
+--      -optc-O3 -optc-mfpmath=sse -optc-msse2 -optc-march=pentium4
+--
+
+import Foreign
+import Foreign.Storable
+import Foreign.Marshal.Alloc
+import Data.IORef
+import Control.Monad
+import System
+import Text.Printf
+
+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"
+    return ()
+
+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 = p `seq` liftM2 (*.) (mass p) (getVec (vel p))
+
+energy :: Double -> Ptr Double -> IO Double
+energy e p
+    | e `seq` p `seq` False = undefined
+    | 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
+    | p1 `seq` m1 `seq` e `seq` p `seq` False = undefined
+    | 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   | p1 `seq` False = undefined
+advance p1 = when (p1 /= end) $ do
+    pos1 <- getVec $ pos p1
+    m1   <- mass p1
+    let go p2
+            | p2 `seq` False = undefined
+            | 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
+
+cursor :: IORef (Ptr Double)
+cursor = unsafePerformIO $ newIORef planets
+
+end :: Ptr Double
+end  = inc planets (nbodies * 7)
+
+next  :: Ptr Double -> Ptr Double
+next = flip inc 7
+
+inc :: Ptr Double -> Int -> Ptr Double
+inc ptr n = ptr `seq` n `seq` plusPtr ptr (n * 8)
+
+newPlanet :: Double -> IO ()
+newPlanet d = d `seq` 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 = flip peekElemOff 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 = p `seq` 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)