Got the nbody demo working, and hopefully faster, too!

1 files changed, 112 insertions, 106 deletions

diff --git a/src/test/bench/shootout/nbody.rs b/src/test/bench/shootout/nbody.rs
index bfa77674..9df439f2 100644
--- a/src/test/bench/shootout/nbody.rs
+++ b/src/test/bench/shootout/nbody.rs
@@ -4,8 +4,6 @@
 
 fn main() {
 
-    log "foo";
-
     let vec[int] inputs = vec(
         50000
         //these segfault :(
@@ -17,13 +15,12 @@ fn main() {
     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;
+            NBodySystem.advance(bodies, 0.01);
+            i += 1;
         }
         log NBodySystem.energy(bodies);
     }
@@ -40,7 +37,6 @@ native "rust" mod rustrt {
 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(), 
@@ -53,10 +49,13 @@ mod NBodySystem {
         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;
+        let int i = 0;
+        while (i < 5) {
+            px += bodies.(i).vx * bodies.(i).mass;
+            py += bodies.(i).vy * bodies.(i).mass;
+            pz += bodies.(i).vz * bodies.(i).mass;
+
+            i += 1;
         }
 
         // side-effecting
@@ -65,15 +64,15 @@ mod NBodySystem {
         ret bodies;
     }
 
-    fn advance(vec[Body.props] bodies, float dt) -> vec[Body.props] {
-        for (Body.props ibody in bodies) {
+    fn advance(vec[Body.props] bodies, float dt) -> () {
 
-            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 int i = 0;
+        while (i < 5) {
+            let int j = i+1;
+            while (j < 5) {
+                let float dx = bodies.(i).x - bodies.(j).x;
+                let float dy = bodies.(i).y - bodies.(j).y;
+                let float dz = bodies.(i).z - bodies.(j).z;
 
                 let float dSquared = dx * dx + dy * dy + dz * dz;
 
@@ -81,23 +80,30 @@ mod NBodySystem {
                 rustrt.squareroot(dSquared, distance);
                 let float mag = dt / (dSquared * distance);
 
-                iBody.vx -= dx * jbody.mass * mag;
-                iBody.vy -= dy * jbody.mass * mag;
-                iBody.vz -= dz * jbody.mass * mag;
+                bodies.(i).vx -= dx * bodies.(j).mass * mag;
+                bodies.(i).vy -= dy * bodies.(j).mass * mag;
+                bodies.(i).vz -= dz * bodies.(j).mass * mag;
+
+                bodies.(j).vx += dx * bodies.(i).mass * mag;
+                bodies.(j).vy += dy * bodies.(i).mass * mag;
+                bodies.(j).vz += dz * bodies.(i).mass * mag;
+
+                j += 1;
 
-                jbody.vx += dx * iBody.mass * mag;
-                jbody.vy += dy * iBody.mass * mag;
-                jbody.vz += dz * iBody.mass * mag;
             }
-        }        
 
-        for (Body.props body in bodies) {
-            body.x += dt * body.vx;
-            body.y += dt * body.vy;
-            body.z += dt * body.vz;
+            i += 1;
         }
 
-        ret bodies;
+        i = 0;
+        while (i < 5) {
+
+            bodies.(i).x += dt * bodies.(i).vx;
+            bodies.(i).y += dt * bodies.(i).vy;
+            bodies.(i).z += dt * bodies.(i).vz;
+
+            i += 1;
+        }
     }
 
     fn energy(vec[Body.props] bodies) -> float {
@@ -107,110 +113,110 @@ mod NBodySystem {
         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;
+        let int i = 0;
+        while (i < 5) {
+            e += 0.5 * bodies.(i).mass *
+                (  bodies.(i).vx * bodies.(i).vx
+                 + bodies.(i).vy * bodies.(i).vy
+                 + bodies.(i).vz * bodies.(i).vz );
+
+            let int j = i+1;
+            while (j < 5) {
+                dx = bodies.(i).x - bodies.(j).x;
+                dy = bodies.(i).y - bodies.(j).y;
+                dz = bodies.(i).z - bodies.(j).z;
 
                 rustrt.squareroot(dx*dx + dy*dy + dz*dz, distance);
-                e -= (iBody.mass * jBody.mass) / distance;
+                e -= (bodies.(i).mass * bodies.(j).mass) / distance;
+                
+                j += 1;
             }
+
+            i += 1;
         }
         ret e;
-    }
 
+    }
 }
 
 mod Body {
-    
-    const float PI = 3.14;
-    const float SOLAR_MASS = 39.47; // was 4 * PI * PI originally
+
+    const float PI = 3.141592653589793;
+    const float SOLAR_MASS = 39.478417604357432; // 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, 
+    type props = rec(mutable float x, 
+                     mutable float y, 
+                     mutable float z, 
+                     mutable float vx, 
+                     mutable float vy,
+                     mutable 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.84e+00;
-        p.y    = -1.16e+00;
-        p.z    = -1.03e-01;
-        p.vx   =  1.66e-03 * DAYS_PER_YEAR;
-        p.vy   =  7.69e-03 * DAYS_PER_YEAR;
-        p.vz   = -6.90e-05 * DAYS_PER_YEAR;
-        p.mass =  9.54e-04 * SOLAR_MASS;
-        ret p;
+        ret rec(
+            mutable x  =  4.84143144246472090e+00,
+            mutable y  = -1.16032004402742839e+00,
+            mutable z  = -1.03622044471123109e-01,
+            mutable vx =  1.66007664274403694e-03 * DAYS_PER_YEAR,
+            mutable vy =  7.69901118419740425e-03 * DAYS_PER_YEAR,
+            mutable vz = -6.90460016972063023e-05 * DAYS_PER_YEAR,
+            mass       =  9.54791938424326609e-04 * SOLAR_MASS
+            );
     }
 
     fn saturn() -> Body.props {
-        let Body.props p;
-        p.x    =  8.34e+00;
-        p.y    =  4.12e+00;
-        p.z    = -4.03e-01;
-        p.vx   = -2.76e-03 * DAYS_PER_YEAR;
-        p.vy   =  4.99e-03 * DAYS_PER_YEAR;
-        p.vz   =  2.30e-05 * DAYS_PER_YEAR;
-        p.mass =  2.85e-04 * SOLAR_MASS;
-        ret p;
+        ret rec(
+            mutable x  =  8.34336671824457987e+00,
+            mutable y  =  4.12479856412430479e+00,
+            mutable z  = -4.03523417114321381e-01,
+            mutable vx = -2.76742510726862411e-03 * DAYS_PER_YEAR,
+            mutable vy =  4.99852801234917238e-03 * DAYS_PER_YEAR,
+            mutable vz =  2.30417297573763929e-05 * DAYS_PER_YEAR,
+            mass       =  2.85885980666130812e-04 * SOLAR_MASS
+            );
    }
 
     fn uranus() -> Body.props {
-        let Body.props p;
-        p.x    =  1.28e+01;
-        p.y    = -1.51e+01;
-        p.z    = -2.23e-01;
-        p.vx   =  2.96e-03 * DAYS_PER_YEAR;
-        p.vy   =  2.37e-03 * DAYS_PER_YEAR;
-        p.vz   = -2.96e-05 * DAYS_PER_YEAR;
-        p.mass =  4.36e-05 * SOLAR_MASS;
-        ret p;
+        ret rec(
+            mutable x  =  1.28943695621391310e+01,
+            mutable y  = -1.51111514016986312e+01,
+            mutable z  = -2.23307578892655734e-01,
+            mutable vx =  2.96460137564761618e-03 * DAYS_PER_YEAR,
+            mutable vy =  2.37847173959480950e-03 * DAYS_PER_YEAR,
+            mutable vz = -2.96589568540237556e-05 * DAYS_PER_YEAR,
+            mass       =  4.36624404335156298e-05 * SOLAR_MASS
+            );
     }
 
     fn neptune() -> Body.props {
-        let Body.props p;
-        p.x    =  1.53e+01;
-        p.y    = -2.59e+01;
-        p.z    =  1.79e-01;
-        p.vx   =  2.68e-03 * DAYS_PER_YEAR;
-        p.vy   =  1.62e-03 * DAYS_PER_YEAR;
-        p.vz   = -9.51e-05 * DAYS_PER_YEAR;
-        p.mass =  5.15e-05 * SOLAR_MASS;
-        ret p;
+        ret rec(
+            mutable x  =  1.53796971148509165e+01,
+            mutable y  = -2.59193146099879641e+01,
+            mutable z  =  1.79258772950371181e-01,
+            mutable vx =  2.68067772490389322e-03 * DAYS_PER_YEAR,
+            mutable vy =  1.62824170038242295e-03 * DAYS_PER_YEAR,
+            mutable vz = -9.51592254519715870e-05 * DAYS_PER_YEAR,
+            mass       =  5.15138902046611451e-05 * SOLAR_MASS
+            );
    }
 
    fn sun() -> Body.props {
-        let Body.props p;
-        p.mass = SOLAR_MASS;
-        ret p;
+       ret rec(
+           mutable x  =  0.0,
+           mutable y  =  0.0,
+           mutable z  =  0.0,
+           mutable vx =  0.0,
+           mutable vy =  0.0,
+           mutable vz =  0.0,
+           mass       =  SOLAR_MASS
+           );
    }
 
-   impure fn offsetMomentum(mutable Body.props props,
-                     float px, 
-                     float py, 
-                     float pz) -> Body.props {
-
+   impure fn offsetMomentum(&Body.props props,
+                            float px, 
+                            float py, 
+                            float pz) -> () {
        props.vx = -px / SOLAR_MASS;
        props.vy = -py / SOLAR_MASS;
        props.vz = -pz / SOLAR_MASS;