Re: [Caml-list] [Benchmark] NBody

[Christoph Bauer <ich@christoph-bauer.net>]

[-- Attachment #1: Type: text/plain, Size: 788 bytes --]

Hi,

>
> The C translation is attached for your amusement.

and here is a mlton version. The test runs on a Pentium-M.

compile flags:

gcc -march=pentium-m -O2 -o nbody-gcc nbody.c
mlton -cc-opt -ffast-math -output nbody-mlton nbody.sml
ocamlopt -ccopt -ffast-math -inline 3 -unsafe    -o nbody-ocaml nbody.ml 

nbody         gcc         ocaml         mlton
1000	      0.000       0.000         0.000
10000         0.014       0.019         0.024
100000        0.070       0.301         0.239
200000        0.231       0.621         0.478
500000        0.552       1.556         1.189
1000000       1.123       3.112         2.377
2000000       2.246       6.229         4.704
3000000       3.367       9.286         8.865
4000000       4.457       12.501        9.687

Christoph Bauer


[-- Attachment #2: nbody.sml --]
[-- Type: application/octet-stream, Size: 4628 bytes --]

val pi = 3.141592653589793
val solar_mass = 4.0 * pi * pi
val days_per_year = 365.24

type planet = {
  x : real ref,
  y : real ref,
  z : real ref,
  vx: real ref,  
  vy: real ref,  
  vz: real ref,
  mass : real
}


fun for (start, stop, f) =
   let
      fun loop i =
         if i > stop
            then ()
         else (f i; loop (i + 1))
   in
      loop start
   end

fun advance (bodies : planet array) dt =
  let 
     val n = Array.length bodies - 1
  in 
   for (0, n, fn i =>
   let 
      val b = Array.sub (bodies, i) 
   in
     for (i+1, n, fn j =>
         let 
           val b' = Array.sub (bodies, j)
           val dx = !(#x b) - !(#x b')  
           val dy = !(#y b) - !(#y b')
           val dz = !(#z b) - !(#z b')
           val distance = Math.sqrt(dx * dx + dy * dy + dz * dz) 
           val mag = dt / (distance * distance * distance) 
         in
           #vx b := !(#vx b) - dx * (#mass b') * mag;
           #vy b := !(#vy b) - dy * (#mass b') * mag;
           #vz b := !(#vz b) - dz * (#mass b') * mag;
	   
           #vx b' := !(#vx b') + dx * (#mass b) * mag;
           #vy b' := !(#vy b') + dy * (#mass b) * mag;
           #vz b' := !(#vz b') + dz * (#mass b) * mag
        end)
  end);

  for( 0, n, fn i =>
    let 
       val b = Array.sub (bodies, i) 
    in
	#x b := !(#x b) + dt * !(#vx b);
	#y b := !(#y b) + dt * !(#vy b);
	#z b := !(#z b) + dt * !(#vz b)
    end)
end

fun energy (bodies : planet array) =
  let 
     val e = ref 0.0 
  in
    for(0, Array.length bodies - 1, fn i =>
    let 
        val b = Array.sub (bodies, i) 
    in
       e := !e + 0.5 * (#mass b) * 
		 (!(#vx b) * !(#vx b) + !(#vy b) 
		  * !(#vy b) + !(#vz b) * !(#vz b));
       for (i+1, Array.length bodies - 1, fn j => 
         let 
            val b' = Array.sub (bodies, j)
            val dx = !(#x b) - !(#x b') 
            val dy = !(#y b) - !(#y b')  
            val dz = !(#z b) - !(#z b')
            val distance = Math.sqrt(dx * dx + dy * dy + dz * dz) 
         in
             e := !e - ((#mass b) * (#mass b')) / distance
         end)
   end);
   !e
  end



fun offset_momentum (bodies : planet array) =
  let 
      val px = ref 0.0
      val py = ref 0.0 
      val pz = ref 0.0
   in
     for (0, Array.length bodies - 1, 
	  fn i =>
	     (px := !px + !(#vx (Array.sub (bodies, i))) * (#mass (Array.sub (bodies, i)));
	      py := !py + !(#vy (Array.sub (bodies, i))) * (#mass (Array.sub (bodies, i)));
	      pz := !pz + !(#vz (Array.sub (bodies, i))) * (#mass (Array.sub (bodies, i)))));
     #vx (Array.sub (bodies, 0)) := ~ (!px / solar_mass);
     #vy (Array.sub (bodies, 0)) := ~ (!py / solar_mass);
     #vz (Array.sub (bodies, 0)) := ~ (!pz / solar_mass)
  end

val jupiter = {
  x = ref 4.84143144246472090,
  y = ref ~1.16032004402742839,
  z = ref ~1.03622044471123109e~1,
  vx = ref (1.66007664274403694e~3 * days_per_year),
  vy = ref (7.69901118419740425e~3 * days_per_year),
  vz = ref (~6.90460016972063023e~5 * days_per_year),
  mass = 9.54791938424326609e~4 * solar_mass
}

val saturn = {
  x = ref 8.34336671824457987,
  y = ref 4.12479856412430479e00,
  z = ref ~4.03523417114321381e~01,
  vx = ref (~2.76742510726862411e~03 * days_per_year),
  vy = ref (4.99852801234917238e~03 * days_per_year),
  vz = ref (2.30417297573763929e~05 * days_per_year),
  mass = 2.85885980666130812e~04 * solar_mass
}

val uranus = {
  x = ref 1.28943695621391310e01,
  y = ref ~1.51111514016986312e01,
  z = ref ~2.23307578892655734e~01,
  vx = ref (2.96460137564761618e~03 * days_per_year),
  vy = ref (2.37847173959480950e~03 * days_per_year),
  vz = ref (~2.96589568540237556e~05 * days_per_year),
  mass = 4.36624404335156298e~05 * solar_mass
}

val neptune = {
  x = ref 1.53796971148509165e01,
  y = ref ~2.59193146099879641e01,
  z = ref 1.79258772950371181e~01,
  vx = ref (2.68067772490389322e~03 * days_per_year),
  vy = ref (1.62824170038242295e~03 * days_per_year),
  vz = ref (~9.51592254519715870e~05 * days_per_year),
  mass = 5.15138902046611451e~05 * solar_mass
}

val sun = {
  x = ref 0.0,  y = ref 0.0,  z = ref 0.0,  vx = ref 0.0,  vy = ref 0.0, vz = ref 0.0,
  mass= solar_mass
}



fun printr r = 
  let 
    val s = Real.fmt (StringCvt.FIX (SOME 9)) r
  in (print s; print "\n") end

val bodies = Array.fromList [ sun, jupiter, saturn, uranus, neptune ]


fun main args =
    let 
	val n = 
	    case Int.fromString (List.hd args) of
		SOME i => i
	      | NONE => 0
			
    in
	offset_momentum bodies;
	printr (energy bodies);
	for (1, n, fn _ => advance bodies 0.01);
	printr (energy bodies)
    end

val _ = main (CommandLine.arguments ())

[-- Attachment #3: Type: text/plain, Size: 172 bytes --]




-- 
let () = let rec f a w i j = Printf.printf "%.20f\r" a; let a1 = a *. i /. j in
if w then f a1 false (i +. 2.0) j else f a1 true i (j +. 2.0) in f 2.0 false 2.0 1.0