{-# OPTIONS_GHC -funbox-strict-fields -fglasgow-exts -fbang-patterns -fexcess-precision -O3 -optc-O2 -optc-ffast-math -optc-mfpmath=sse -optc-msse2 -optc-march=pentium4 #-}


{-  
    The Computer Language Shootout
       http://shootout.alioth.debian.org/
       
    Contributed by Olof Kraigher, with help from Don Stewart.
    Munged by Lennart Augustsson.
    
    Compile with:
        -funfolding-keeness-factor=10 -funfolding-use-threshold=2
-}

{-
coxorange:~/Proj] lennart% time ./a.out 2000000
-0.169075164
-0.169026286
2.049u 0.009s 0:02.06 99.0%     0+0k 0+0io 0pf+0w
-}
 
import Control.Monad
import System
import Text.Printf
import Control.Monad.ST
import Data.STRef
import Data.Array.MArray
import Data.Array.Base(unsafeWrite, unsafeRead)
import Data.Array.ST

data Vector3 = Vec !Double !Double !Double
 
nr_of_planets :: Int
nr_of_planets = 5
 
delta_t :: Double
delta_t = 0.01
 
days_per_year :: Double
days_per_year = 365.24
 
solar_mass :: Double
solar_mass = 4 * pi ** 2

type Bodies s = STUArray s Int Double
 
end :: Int
end = nr_of_planets * 7
 
nextBody ptr = ptr + 7
 
{- # INLINE getVec # -}
getVec planets ptr =
    let p o = unsafeRead planets (ptr+o) in
    liftM3 Vec (p 0) (p 1) (p 2)
 
{- # INLINE setVec # -}
setVec :: Bodies s -> Int -> Vector3 -> ST s ()
setVec planets ptr (Vec x y z) = do
    unsafeWrite planets (ptr + 0) x
    unsafeWrite planets (ptr + 1) y
    unsafeWrite planets (ptr + 2) z
 
mass planets ptr = unsafeRead planets (ptr + 6)
 
pos ptr = ptr + 0
 
vel ptr = ptr + 3
 
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)
 
magnitude2 (Vec x y z) = x*x + y*y + z*z

bodyData = 
    let dp = days_per_year in
        [   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
        ]

initialize :: ST s (Bodies s)
initialize = newListArray (0, length bodyData - 1) bodyData
 
loopOverBodies start body =
    let loop idx
            | idx /= end = do body idx; loop (nextBody idx)
	    | otherwise  = return ()
    in  loop start

advance planets = do
    loopOverBodies 0 $ \ planet1 -> do
        p1 <- getVec planets $ pos planet1
        m1  <- mass planets planet1
        let vel1 = vel planet1
        
	    v1 += v2 = do 
	        v1' <- getVec planets v1
		setVec planets v1 $ v1' .+. v2

	    v1 -= v2 = do 
	        v1' <- getVec planets v1
		setVec planets v1 $ v1' .-. v2

        loopOverBodies (nextBody planet1) $ \ planet2 -> do
            p2 <- getVec planets $ pos planet2
            m2 <- mass planets planet2
            let vel2 = vel planet2
                difference = p1 .-. p2
                distance2 = magnitude2 difference
                distance = sqrt distance2
                magnitude = delta_t / (distance2 * distance)
                mass_magn = magnitude *. difference
            vel1 -= (m2 *. mass_magn)
            vel2 += (m1 *. mass_magn)
        v1 <- getVec planets vel1
        planet1 += (delta_t *. v1)
 
energy planets = do
    e <- newSTRef 0
    
    loopOverBodies 0 $ \ planet1 -> do
        p1 <- getVec planets $ pos planet1
        v1 <- getVec planets $ vel planet1
        m1 <- mass planets planet1

	loopOverBodies (nextBody planet1) $ \ planet2 -> do
            p2 <- getVec planets $ pos planet2
            v2 <- getVec planets $ vel planet2
            m2 <- mass planets planet2
            let distance = sqrt $ magnitude2 $ p1 .-. p2
            modifySTRef e (\x -> x - m1 * m2 / distance)
	modifySTRef e (+ 0.5 * m1 * magnitude2 v1)

    readSTRef e
 
offset_momentum planets = do
    let momentum' planet = do
        v <- getVec planets $ vel planet
        m <- mass planets planet
        return $ m *. v
    
    momentum <- return . foldl (.+.) (Vec 0 0 0) =<< 
        (mapM momentum' $ take (nr_of_planets - 1) $ tail $ iterate nextBody 0)
 
    setVec planets (vel 0) $ (-1/solar_mass) *. momentum
 
computeST n = do
    planets <- initialize
    offset_momentum planets
    e1 <- energy planets
    replicateM_ n (advance planets)
    e2 <- energy planets
    return (e1, e2)

compute n = runST (computeST n)

main = do
    n <- getArgs >>= readIO.head

    let (e1, e2) = compute n
    printf "%.9f\n" e1
    printf "%.9f\n" e2