The Computer Language
Benchmarks Game

n-body Rust #3 program

source code

//
// The Computer Language Benchmarks Game
// http://benchmarksgame.alioth.debian.org/
//
// Contributed by Alex Drozhak
//

use std::ops::{Add, Sub, Mul, AddAssign, SubAssign};
use std::f64::consts::PI;

const SOLAR_MASS: f64 = 4.0 * PI * PI;
const YEAR: f64 = 365.24;
const N_BODIES: usize = 5;

#[derive(Clone, Copy)]
struct Vec3(pub f64, pub f64, pub f64);

impl Vec3 {
    fn new() -> Self {
        Vec3(0.0, 0.0, 0.0)
    }

    fn norm(&self) -> f64 {
        self.squared_norm().sqrt()
    }

    fn squared_norm(&self) -> f64 {
        self.0 * self.0 + self.1 * self.1 + self.2 * self.2
    }
}

impl Add for Vec3 {
    type Output = Self;
    fn add(self, rhs: Self) -> Self {
        Vec3(self.0 + rhs.0, self.1 + rhs.1, self.2 + rhs.2)
    }
}

impl Sub for Vec3 {
    type Output = Self;
    fn sub(self, rhs: Self) -> Self {
        Vec3(self.0 - rhs.0, self.1 - rhs.1, self.2 - rhs.2)
    }
}

impl AddAssign for Vec3 {
    fn add_assign(&mut self, rhs: Self) {
        *self = Vec3(self.0 + rhs.0, self.1 + rhs.1, self.2 + rhs.2);
    }
}

impl SubAssign for Vec3 {
    fn sub_assign(&mut self, rhs: Self) {
        *self = Vec3(self.0 - rhs.0, self.1 - rhs.1, self.2 - rhs.2);
    }
}

impl Mul<f64> for Vec3 {
    type Output = Self;
    fn mul(self, rhs: f64) -> Self {
        Vec3(self.0 * rhs, self.1 * rhs, self.2 * rhs)
    }
}

#[derive(Clone, Copy)]
struct Planet {
    pos: Vec3,
    vel: Vec3,
    mass: f64,
}

impl Planet {
    fn sun() -> Self {
        Planet {
            pos: Vec3(0.0, 0.0, 0.0),
            vel: Vec3(0.0, 0.0, 0.0),
            mass: SOLAR_MASS,
        }
    }

    fn jupiter() -> Self {
        Planet {
            pos: Vec3(4.84143144246472090e+00,
                      -1.16032004402742839e+00,
                      -1.03622044471123109e-01),
            vel: Vec3(1.66007664274403694e-03 * YEAR,
                      7.69901118419740425e-03 * YEAR,
                      -6.90460016972063023e-05 * YEAR),
            mass: 9.54791938424326609e-04 * SOLAR_MASS,
        }
    }

    fn saturn() -> Self {
        Planet {
            pos: Vec3(8.34336671824457987e+00,
                      4.12479856412430479e+00,
                      -4.03523417114321381e-01),
            vel: Vec3(-2.76742510726862411e-03 * YEAR,
                      4.99852801234917238e-03 * YEAR,
                      2.30417297573763929e-05 * YEAR),
            mass: 2.85885980666130812e-04 * SOLAR_MASS,
        }
    }

    fn uranus() -> Self {
        Planet {
            pos: Vec3(1.28943695621391310e+01,
                      -1.51111514016986312e+01,
                      -2.23307578892655734e-01),
            vel: Vec3(2.96460137564761618e-03 * YEAR,
                      2.37847173959480950e-03 * YEAR,
                      -2.96589568540237556e-05 * YEAR),
            mass: 4.36624404335156298e-05 * SOLAR_MASS,
        }
    }

    fn neptune() -> Self {
        Planet {
            pos: Vec3(1.53796971148509165e+01,
                      -2.59193146099879641e+01,
                      1.79258772950371181e-01),
            vel: Vec3(2.68067772490389322e-03 * YEAR,
                      1.62824170038242295e-03 * YEAR,
                      -9.51592254519715870e-05 * YEAR),
            mass: 5.15138902046611451e-05 * SOLAR_MASS,
        }
    }
}

struct NBSystem {
    planets: [Planet; N_BODIES],
}

impl NBSystem {
    fn new() -> Self {
        NBSystem {
            planets: [Planet::sun(),
                      Planet::jupiter(),
                      Planet::saturn(),
                      Planet::uranus(),
                      Planet::neptune()],
        }
    }

    fn offset_momentum(&mut self) {
        let p = self.planets.iter_mut().fold(Vec3::new(), |v, b| v + b.vel * b.mass);
        self.planets[0].vel = p * (-1.0 / self.planets[0].mass);
    }

    fn energy(&self) -> f64 {
        let mut e = 0.0;
        let mut bodies = self.planets.iter();
        while let Some(bi) = bodies.next() {
            e +=
                bi.vel.squared_norm() * bi.mass / 2.0 -
                bi.mass *
                bodies.clone().map(|bj| bj.mass / (bi.pos - bj.pos).norm()).fold(0.0, |a, b| a + b);
        }
        e
    }

    fn advance(&mut self, dt: f64) {
        let planets = &mut self.planets;
        for i in 0..N_BODIES {
            for j in (i + 1)..N_BODIES {
                let dp = planets[i].pos - planets[j].pos;

                let distance = dp.norm();
                let mag = dt / (distance * distance * distance);

                planets[i].vel -= dp * planets[j].mass * mag;
                planets[j].vel += dp * planets[i].mass * mag;
            }
            planets[i].pos += planets[i].vel * dt;
        }
    }
}

fn main() {
    let n = std::env::args_os()
        .nth(1)
        .and_then(|s| s.into_string().ok())
        .and_then(|n| n.parse().ok())
        .unwrap_or(1000);
    let mut system = NBSystem::new();
    system.offset_momentum();
    println!("{:.9}", system.energy());
    for _ in 0..n {
        system.advance(0.01);
    }
    println!("{:.9}", system.energy());
}
    

notes, command-line, and program output

NOTES:
64-bit Ubuntu quad core
rustc 1.21.0 (3b72af97e 2017-10-09)


Sun, 15 Oct 2017 18:15:32 GMT

MAKE:
/usr/local/src/rust-1.21.0/bin/rustc -C opt-level=3 -C target-cpu=core2 -C lto -C llvm-args='-unroll-threshold=500' nbody.rs -o nbody.rust-3.rust_run

5.51s to complete and log all make actions

COMMAND LINE:
./nbody.rust-3.rust_run 50000000

PROGRAM OUTPUT:
-0.169075164
-0.169059907