The Computer Language
23.03 Benchmarks Game

n-body Rust #4 program

source code

//
// The Computer Language Benchmarks Game
// https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
//
// contributed by Alex Drozhak
// modified by Tung Duong
//

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;
const N_PAIRS: usize = N_BODIES*(N_BODIES-1)/2;

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

#[derive(Clone, Copy)]
struct Pair(pub usize, pub usize);

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],
    pairs : [Pair; N_PAIRS],
    offsets : [Vec3; N_PAIRS],
    distances_square :[f64;N_PAIRS],
    distances :[f64;N_PAIRS],
    mags :[f64;N_PAIRS]
}

impl NBSystem {
    fn new() -> Self {
        let mut system = NBSystem {
            planets: [Planet::sun(),
                      Planet::jupiter(),
                      Planet::saturn(),
                      Planet::uranus(),
                      Planet::neptune()],
            pairs : [Pair(0,0); N_PAIRS],
            offsets : [Vec3(0.0,0.0,0.0);N_PAIRS],
            distances_square : [0.0;N_PAIRS],
            distances : [0.0;N_PAIRS],
            mags : [0.0;N_PAIRS],
        };
        let mut index = 0;
        for i in 0..N_BODIES {
            for j in (i + 1)..N_BODIES {
				system.pairs[index].0 = i;
				system.pairs[index].1 = j;
				index += 1;
			}
		}
        system
    }

    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 (index, &Pair(i,j)) in self.pairs.iter().enumerate() {
			self.offsets[index] = planets[i].pos - planets[j].pos;
		}
		for index in 0..N_PAIRS {
			self.distances_square[index] = self.offsets[index].squared_norm();
			self.distances[index] = self.distances_square[index].sqrt();
		}
		for index in 0..N_PAIRS {
			self.mags[index] = dt / (self.distances_square[index] * self.distances[index]);
		}
        for (index, &Pair(i,j)) in self.pairs.iter().enumerate() {
		    planets[i].vel -= self.offsets[index] * planets[j].mass * self.mags[index];
		    planets[j].vel += self.offsets[index] * planets[i].mass * self.mags[index];
		}
        for planet in planets.iter_mut() {
            planet.pos += planet.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
1.67.0
(fc594f156
2023-01-24)
LLVM version: 15.0.6


Thu, 26 Jan 2023 21:01:56 GMT

MAKE:
/opt/src/rust-1.67.0/bin/rustc -C opt-level=3 -C target-cpu=ivybridge -C codegen-units=1  nbody.rs -o nbody.rust-4.rust_run

8.66s to complete and log all make actions

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

PROGRAM OUTPUT:
-0.169075164
-0.169059907