A first stab at the 'nbody' demo, with support for calling C sqrt().

5 files changed, 219 insertions, 1 deletions

diff --git a/Makefile.in b/Makefile.in
index 129b46a81a0..5ddf5b784e7 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -638,7 +638,8 @@ TEST_XFAILS_BOOT := $(TASK_XFAILS) \
                       bad-send.rs \
                       infinite-vec-type-recursion.rs \
                       tail-non-call.rs \
-                      writing-through-read-alias.rs)
+                      writing-through-read-alias.rs) \
+                    $(S)src/test/bench/shootout/nbody.rs
 
 TEST_XFAILS_STAGE0 := $(FLOAT_XFAILS) \
                       $(addprefix $(S)src/test/run-pass/, \
diff --git a/src/rt/rust_builtin.cpp b/src/rt/rust_builtin.cpp
index 983337745b9..2de9d7f47ac 100644
--- a/src/rt/rust_builtin.cpp
+++ b/src/rt/rust_builtin.cpp
@@ -51,6 +51,11 @@ last_os_error(rust_task *task) {
     return st;
 }
 
+extern "C" CDECL
+void squareroot(rust_task *task, double *input, double *output) {
+    *output = sqrt(*input);
+}
+
 extern "C" CDECL size_t
 size_of(rust_task *task, type_desc *t) {
   return t->size;
diff --git a/src/rt/rust_internal.h b/src/rt/rust_internal.h
index 5bd62d2b7a4..a1cc8e514b0 100644
--- a/src/rt/rust_internal.h
+++ b/src/rt/rust_internal.h
@@ -16,6 +16,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <fcntl.h>
+#include <math.h>
 
 #include "rust.h"
 #include "rand.h"
diff --git a/src/rt/rustrt.def.in b/src/rt/rustrt.def.in
index f8ac1ed46ef..df9e26aff9e 100644
--- a/src/rt/rustrt.def.in
+++ b/src/rt/rustrt.def.in
@@ -20,6 +20,7 @@ rust_get_stdout
 rust_run_program
 rust_start
 size_of
+squareroot
 str_alloc
 str_buf
 str_byte_len
@@ -39,6 +40,7 @@ upcall_get_type_desc
 upcall_grow_task
 upcall_join
 upcall_kill
+upcall_log_double
 upcall_log_int
 upcall_log_str
 upcall_malloc
diff --git a/src/test/bench/shootout/nbody.rs b/src/test/bench/shootout/nbody.rs
new file mode 100644
index 00000000000..b0ac2d44697
--- /dev/null
+++ b/src/test/bench/shootout/nbody.rs
@@ -0,0 +1,209 @@
+// based on:
+// http://shootout.alioth.debian.org/u32/benchmark.php?test=nbody&lang=java
+
+fn main() {
+
+    let vec[int] inputs = vec(
+        50000
+        //these segfault :(
+        //500000,
+        //5000000,
+        //50000000
+        );
+
+    let vec[Body.props] bodies = NBodySystem.MakeNBodySystem();
+
+    for (int n in inputs) {
+        // TODO: make #fmt handle floats?
+        log NBodySystem.energy(bodies);
+
+        let int i = 0;
+        while (i < n) {
+            bodies = NBodySystem.advance(bodies, 0.01);
+            i = i+1;
+        }
+        log NBodySystem.energy(bodies);
+    }
+}
+
+// making a native call to sqrt
+native "rust" mod rustrt {
+    fn squareroot(&float input, &mutable float output);
+}
+
+// Body.props is a record of floats, so
+// vec[Body.props] is a vector of records of floats
+
+mod NBodySystem {
+
+    fn MakeNBodySystem() -> vec[Body.props] {
+        // can't iterate over a record?  how about a vector, then?
+        let vec[Body.props] bodies = vec(
+            // these each return a Body.props
+            Body.sun(), 
+            Body.jupiter(), 
+            Body.saturn(), 
+            Body.uranus(), 
+            Body.neptune());
+
+        let float px = 0.0;
+        let float py = 0.0;
+        let float pz = 0.0;
+
+        for (Body.props body in bodies) {
+            px += body.vx * body.mass;
+            py += body.vy * body.mass;
+            pz += body.vz * body.mass;
+        }
+        bodies.(0) = Body.offsetMomentum(bodies.(0), px, py, pz);
+
+        ret bodies;
+    }
+
+    fn advance(vec[Body.props] bodies, float dt) -> vec[Body.props] {
+        for (Body.props ibody in bodies) {
+
+            let Body.props iBody = ibody;
+
+            for (Body.props jbody in bodies) {
+                let float dx = iBody.x - jbody.x;
+                let float dy = iBody.y - jbody.y;
+                let float dz = iBody.z - jbody.z;
+
+                let float dSquared = dx * dx + dy * dy + dz * dz;
+
+                let float distance;
+                rustrt.squareroot(dSquared, distance);
+                let float mag = dt / (dSquared * distance);
+            }
+        }        
+
+        for (Body.props body in bodies) {
+            body.x += dt * body.vx;
+            body.y += dt * body.vy;
+            body.z += dt * body.vz;
+        }
+
+        ret bodies;
+    }
+
+    fn energy(vec[Body.props] bodies) -> float {
+        let float dx;
+        let float dy;
+        let float dz;
+        let float distance;
+        let float e = 0.0;
+
+        for (Body.props ibody in bodies) {
+
+            // do we need this?
+            let Body.props iBody = ibody;
+
+            e += 0.5 * iBody.mass *
+                ( iBody.vx * iBody.vx
+                  + iBody.vy * iBody.vy
+                  + iBody.vz * iBody.vz );
+
+            for (Body.props jbody in bodies) {
+
+                // do we need this?
+                let Body.props jBody = jbody;
+
+                dx = iBody.x - jBody.x;
+                dy = iBody.y - jBody.y;
+                dz = iBody.z - jBody.z;
+
+                rustrt.squareroot(dx*dx + dy*dy + dz*dz, distance);
+                e -= (iBody.mass * jBody.mass) / distance;
+            }
+        }
+        ret e;
+    }
+
+}
+
+mod Body {
+    
+    const float PI = 3.141592;
+    const float SOLAR_MASS = 39.478417; // was 4 * PI * PI originally
+    const float DAYS_PER_YEAR = 365.24;
+
+    type props = rec(float x, 
+                     float y, 
+                     float z, 
+                     float vx, 
+                     float vy,
+                     float vz, 
+                     float mass);
+
+    fn jupiter() -> Body.props {
+        // current limitation of the float lexer: decimal part has to
+        // fit into a 32-bit int.
+        
+        let Body.props p;
+        p.x    =  4.841431e+00;
+        p.y    = -1.160320e+00;
+        p.z    = -1.036220e-01;
+        p.vx   =  1.660076e-03 * DAYS_PER_YEAR;
+        p.vy   =  7.699011e-03 * DAYS_PER_YEAR;
+        p.vz   = -6.904600e-05 * DAYS_PER_YEAR;
+        p.mass =  9.547919e-04 * SOLAR_MASS;
+        ret p;
+    }
+
+    fn saturn() -> Body.props {
+        let Body.props p;
+        p.x    =  8.343366e+00;
+        p.y    =  4.124798e+00;
+        p.z    = -4.035234e-01;
+        p.vx   = -2.767425e-03 * DAYS_PER_YEAR;
+        p.vy   =  4.998528e-03 * DAYS_PER_YEAR;
+        p.vz   =  2.304172e-05 * DAYS_PER_YEAR;
+        p.mass =  2.858859e-04 * SOLAR_MASS;
+        ret p;
+   }
+
+    fn uranus() -> Body.props {
+        let Body.props p;
+        p.x    =  1.289436e+01;
+        p.y    = -1.511115e+01;
+        p.z    = -2.233075e-01;
+        p.vx   =  2.964601e-03 * DAYS_PER_YEAR;
+        p.vy   =  2.378471e-03 * DAYS_PER_YEAR;
+        p.vz   = -2.965895e-05 * DAYS_PER_YEAR;
+        p.mass =  4.366244e-05 * SOLAR_MASS;
+        ret p;
+    }
+
+    fn neptune() -> Body.props {
+        let Body.props p;
+        p.x    =  1.537969e+01;
+        p.y    = -2.591931e+01;
+        p.z    =  1.792587e-01;
+        p.vx   =  2.680677e-03 * DAYS_PER_YEAR;
+        p.vy   =  1.628241e-03 * DAYS_PER_YEAR;
+        p.vz   = -9.515922e-05 * DAYS_PER_YEAR;
+        p.mass =  5.151389e-05 * SOLAR_MASS;
+        ret p;
+   }
+
+   fn sun() -> Body.props {
+        let Body.props p;
+        p.mass = SOLAR_MASS;
+        ret p;
+   }
+
+   fn offsetMomentum(Body.props props,
+                     float px, 
+                     float py, 
+                     float pz) -> Body.props {
+
+       // TODO: should we create a new one or mutate the original?
+       let Body.props p = props;
+       p.vx = -px / SOLAR_MASS;
+       p.vy = -py / SOLAR_MASS;
+       p.vz = -pz / SOLAR_MASS;
+       ret p;
+   }
+
+}