tutorials/learn-rust-the-dangerous-way
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Create ideomatic rust version

Browse Source
This commit is contained in:
Timothy Warren 2020-01-13 16:20:12 -05:00
parent fd41a5dbb0
commit 5d27ffc192
3 changed files with 134 additions and 234 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.cargo/config Normal file
View File

@ -0,0 +1,2 @@
[build]
rustflags = [ "-C", "target-cpu=core2" ]

9
Cargo.toml
View File

@ -1,9 +1,14 @@
[package]
name = "lrtdw"
name = "nbody"
version = "0.1.0"
authors = ["Timothy Warren <twarren@nabancard.com>"]
authors = ["Timothy Warren <tim@timshomepage.net>"]
edition = "2018"
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[dependencies]
[profile.release]
opt-level = 3
codegen-units = 1
debug = 2

343
src/main.rs
View File

@ -1,310 +1,203 @@
// The Computer Language Benchmarks Game
// https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
//
// Contributed by Mark C. Lewis.
// Modified slightly by Chad Whipkey.
// Converted from Java to C++ and added SSE support by Branimir Maksimovic.
// Converted from C++ to C by Alexey Medvedchikov.
// Modified by Jeremy Zerfas.
// Converted to Rust by Cliff L. Biffle
//! n-body simulation in Rust - naive version
//!
//! This program knows nothing about vector units, alignment, locality, and the
//! like. It does the math in the simplest way I could come up with, and relies
//! on the compiler to make it fast.
#![allow(non_upper_case_globals, non_camel_case_types, non_snake_case)]
use std::arch::x86_64::*;
use std::f64::consts::PI;
#[repr(C)]
struct body {
#[derive(Clone, Debug)]
struct Body {
position: [f64; 3],
velocity: [f64; 3],
mass: f64,
}
const SOLAR_MASS: f64 = 4. * PI * PI;
const DAYS_PER_YEAR: f64 = 365.24;
/// Number of bodies modeled in the simulation.
const BODIES_COUNT: usize = 5;
const STARTING_STATE: [body; BODIES_COUNT] = [
body {
const SOLAR_MASS: f64 = (4. * PI * PI);
const DAYS_PER_YEAR: f64 = 365.24;
/// Number of body-body interactions.
const INTERACTIONS: usize = BODIES_COUNT * (BODIES_COUNT - 1) / 2;
/// Initial state of the simulation.
const STARTING_STATE: [Body; BODIES_COUNT] = [
// Sun
Body {
mass: SOLAR_MASS,
position: [0.; 3],
velocity: [0.; 3],
},
body {
// Jupiter
mass: 9.54791938424326609e-04 * SOLAR_MASS,
Body {
position: [
4.84143144246472090e+00,
-1.16032004402742839e+00,
-1.03622044471123109e-01,
4.841_431_442_464_72e0,
-1.160_320_044_027_428_4e0,
-1.036_220_444_711_231_1e-1,
],
velocity: [
1.66007664274403694e-03 * DAYS_PER_YEAR,
7.69901118419740425e-03 * DAYS_PER_YEAR,
-6.90460016972063023e-05 * DAYS_PER_YEAR,
1.660_076_642_744_037e-3 * DAYS_PER_YEAR,
7.699_011_184_197_404e-3 * DAYS_PER_YEAR,
-6.904_600_169_720_63e-5 * DAYS_PER_YEAR,
],
mass: 9.547_919_384_243_266e-4 * SOLAR_MASS,
},
body {
// Saturn
mass: 2.85885980666130812e-04 * SOLAR_MASS,
Body {
position: [
8.34336671824457987e+00,
4.12479856412430479e+00,
-4.03523417114321381e-01,
8.343_366_718_244_58e0,
4.124_798_564_124_305e0,
-4.035_234_171_143_214e-1,
],
velocity: [
-2.76742510726862411e-03 * DAYS_PER_YEAR,
4.99852801234917238e-03 * DAYS_PER_YEAR,
2.30417297573763929e-05 * DAYS_PER_YEAR,
-2.767_425_107_268_624e-3 * DAYS_PER_YEAR,
4.998_528_012_349_172e-3 * DAYS_PER_YEAR,
2.304_172_975_737_639_3e-5 * DAYS_PER_YEAR,
],
mass: 2.858_859_806_661_308e-4 * SOLAR_MASS,
},
body {
// Uranus
mass: 4.36624404335156298e-05 * SOLAR_MASS,
Body {
position: [
1.28943695621391310e+01,
-1.51111514016986312e+01,
-2.23307578892655734e-01,
1.289_436_956_213_913_1e1,
-1.511_115_140_169_863_1e1,
-2.233_075_788_926_557_3e-1,
],
velocity: [
2.96460137564761618e-03 * DAYS_PER_YEAR,
2.37847173959480950e-03 * DAYS_PER_YEAR,
-2.96589568540237556e-05 * DAYS_PER_YEAR,
2.964_601_375_647_616e-3 * DAYS_PER_YEAR,
2.378_471_739_594_809_5e-3 * DAYS_PER_YEAR,
-2.965_895_685_402_375_6e-5 * DAYS_PER_YEAR,
],
mass: 4.366_244_043_351_563e-5 * SOLAR_MASS,
},
body {
// Neptune
mass: 5.15138902046611451e-05 * SOLAR_MASS,
Body {
position: [
1.53796971148509165e+01,
-2.59193146099879641e+01,
1.79258772950371181e-01,
1.537_969_711_485_091_1e1,
-2.591_931_460_998_796_4e1,
1.792_587_729_503_711_8e-1,
],
velocity: [
2.68067772490389322e-03 * DAYS_PER_YEAR,
1.62824170038242295e-03 * DAYS_PER_YEAR,
-9.51592254519715870e-05 * DAYS_PER_YEAR,
2.680_677_724_903_893_2e-3 * DAYS_PER_YEAR,
1.628_241_700_382_423e-3 * DAYS_PER_YEAR,
-9.515_922_545_197_159e-5 * DAYS_PER_YEAR,
],
mass: 5.151_389_020_466_114_5e-5 * SOLAR_MASS,
},
];
// Figure out how many total different interactions there are between each
// body and every other body. Some of the calculations for these
// interactions will be calculated two at a time by using x86 SSE
// instructions and because of that it will also be useful to have a
// ROUNDED_INTERACTIONS_COUNT that is equal to the next highest even number
// which is equal to or greater than INTERACTIONS_COUNT.
const INTERACTIONS_COUNT: usize = BODIES_COUNT * (BODIES_COUNT - 1) / 2;
const ROUNDED_INTERACTIONS_COUNT: usize = INTERACTIONS_COUNT + INTERACTIONS_COUNT % 2;
/// Adjust the Sun's velocity to offset system momentum.
fn offset_momentum(bodies: &mut [Body; BODIES_COUNT]) {
let (sun, planets) = bodies.split_first_mut().unwrap();
// It's useful to have two arrays to keep track of the position_Deltas
// and magnitudes of force between the bodies for each interaction. For the
// position_Deltas array, instead of using a one dimensional array of
// structures that each contain the X, Y, and Z components for a position
// delta, a two dimensional array is used instead which consists of three
// arrays that each contain all of the X, Y, and Z components for all of the
// position_Deltas. This allows for more efficient loading of this data into
// SSE registers. Both of these arrays are also set to contain
// ROUNDED_INTERACTIONS_COUNT elements to simplify one of the following
// loops and to also keep the second and third arrays in position_Deltas
// aligned properly.
#[derive(Copy, Clone)]
#[repr(C)]
union Interactions {
scalars: [f64; ROUNDED_INTERACTIONS_COUNT],
vectors: [__m128d; ROUNDED_INTERACTIONS_COUNT / 2],
}
impl Interactions {
/// Returns a refrence to the storage as `f64`s.
pub fn as_scalars(&mut self) -> &mut [f64; ROUNDED_INTERACTIONS_COUNT] {
// Safety: the in-memory representation of `f64` and `__m128d` is
// compatible, so accesses to the union members is safe in any
// order..
unsafe { &mut self.scalars }
}
/// Returns a reference to the storage as `__m128d`s.
pub fn as_vectors(&mut self) -> &mut [__m128d; ROUNDED_INTERACTIONS_COUNT / 2] {
// Safety: the in-memory representation of `f64` and `__m128d` is
// compatible, so accesses to the union members is safe in any
// order..
unsafe { &mut self.vectors }
}
}
// Calculate the momentum of each body and conserve momentum of the system by
// adding to the Sun's velocity the appropriate opposite velocity needed in
// order to offset that body's momentum.
fn offset_Momentum(bodies: &mut [body; BODIES_COUNT]) {
for i in 0..BODIES_COUNT {
sun.velocity = [0.; 3];
for planet in planets {
for m in 0..3 {
bodies[0].velocity[m] -= bodies[i].velocity[m] * bodies[i].mass / SOLAR_MASS;
sun.velocity[m] -= planet.velocity[m] * planet.mass / SOLAR_MASS;
}
}
}
// Output the total energy of the system.
fn output_Energy(bodies: &mut [body; BODIES_COUNT]) {
/// A convenient way of computing `x` squared
fn sqr(x: f64) -> f64 {
x * x
}
/// Print the system energy.
fn output_energy(bodies: &mut [Body; BODIES_COUNT]) {
let mut energy = 0.;
for i in 0..BODIES_COUNT {
for (i, body) in bodies.iter().enumerate() {
// Add the kinetic energy for each body.
energy += 0.5
* bodies[i].mass
* (bodies[i].velocity[0] * bodies[i].velocity[0]
+ bodies[i].velocity[1] * bodies[i].velocity[1]
+ bodies[i].velocity[2] * bodies[i].velocity[2]);
energy +=
0.5 * body.mass * (
sqr(body.velocity[0])
+ sqr(body.velocity[1])
+ sqr(body.velocity[2])
);
// Add the potential energy between this body and
// every other body
for j in i + 1..BODIES_COUNT {
let mut position_Delta = [0.; 3];
for m in 0..3 {
position_Delta[m] = bodies[i].position[m] - bodies[j].position[m];
}
energy -= bodies[i].mass * bodies[j].mass
/ f64::sqrt(
position_Delta[0] * position_Delta[0]
+ position_Delta[1] * position_Delta[1]
+ position_Delta[2] * position_Delta[2],
// Add the potential energy between this body and every other body.
for body2 in &bodies[i + 1..BODIES_COUNT] {
energy -= body.mass * body2.mass /
f64::sqrt(
sqr(body.position[0] - body2.position[0])
+ sqr(body.position[1] - body2.position[1])
+ sqr(body.position[2] - body2.position[2]),
);
}
}
// Output the total energy of the system
println!("{:.9}", energy);
}
}
// Advance all the bodies in the system by one timestep. Calculate the
// interactions between all the bodies, update each body's velocity based on
// those interactions, and update each body's position by the distance it
// travels in a timestep at it's updated velocity.
#[cfg(target_feature = "sse2")]
fn advance(
bodies: &mut [body; BODIES_COUNT],
position_Deltas: &mut [Interactions; 3],
magnitudes: &mut Interactions,
) {
// Calculate the position_Deltas between the bodies for each interaction.
/// Steps the simulation foward by one time step
fn advance(bodies: &mut [Body; BODIES_COUNT]) {
// Compute point-to-point vectors between each unique pair of bodies.
// Note: this array is large enough that computing it mutable and returning
// it from a block, as I did with magnitudes below, generates a memcpy.
// Sigh. So I'll leave it mutable.
let mut position_deltas = [[0.;3]; INTERACTIONS];
{
let mut k = 0;
for i in 0..BODIES_COUNT - 1 {
for j in i + 1..BODIES_COUNT {
for m in 0..3 {
position_Deltas[m].as_scalars()[k] =
bodies[i].position[m] - bodies[j].position[m];
for i in 0..BODIES_COUNT-1 {
for j in i+1..BODIES_COUNT {
for (m, pd) in position_deltas[k].iter_mut().enumerate() {
*pd = bodies[i].position[m] - bodies[j].position[m];
}
k += 1;
}
}
}
// Calculate the magnitudes of force between the bodies for each
// interaction. This loop processes two interactions at a time which is why
// ROUNDED_INTERACTIONS_COUNT/2 iterations are done.
for i in 0..ROUNDED_INTERACTIONS_COUNT / 2 {
// Load position_Deltas of two bodies into position_Delta.
let mut position_Delta = [unsafe { _mm_setzero_pd() }; 3];
// Compute the `i/d^3` magnitude between each pair of bodies.
let magnitudes = {
let mut magnitudes = [0.; INTERACTIONS];
for (i, mag) in magnitudes.iter_mut().enumerate() {
let distance_squared = sqr(position_deltas[i][0])
+ sqr(position_deltas[i][1])
+ sqr(position_deltas[i][2]);
for m in 0..3 {
position_Delta[m] = position_Deltas[m].as_vectors()[i];
*mag = 0.01 / (distance_squared * distance_squared.sqrt());
}
let distance_Squared: __m128d = unsafe {
_mm_add_pd(
_mm_add_pd(
_mm_mul_pd(position_Delta[0], position_Delta[0]),
_mm_mul_pd(position_Delta[1], position_Delta[1]),
),
_mm_mul_pd(position_Delta[2], position_Delta[2]),
)
magnitudes
};
// Doing square roots normally using double precision floating point
// math can be quite time consuming so SSE's much faster single
// precision reciprocal square root approximation instruction is used as
// a starting point instead. The precision isn't quite sufficient to get
// acceptable results so two iterations of the Newton–Raphson method are
// done to improve precision further.
let mut distance_Reciprocal: __m128d =
unsafe { _mm_cvtps_pd(_mm_rsqrt_ps(_mm_cvtpd_ps(distance_Squared))) };
for _ in 0..2 {
// Normally the last four multiplications in this equation would
// have to be done sequentially but by placing the last
// multiplication in parentheses, a compiler can then schedule that
// multiplication earlier.
distance_Reciprocal = unsafe {
_mm_sub_pd(
_mm_mul_pd(distance_Reciprocal, _mm_set1_pd(1.5)),
_mm_mul_pd(
_mm_mul_pd(
_mm_mul_pd(_mm_set1_pd(0.5), distance_Squared),
distance_Reciprocal,
),
_mm_mul_pd(distance_Reciprocal, distance_Reciprocal),
),
)
};
}
// Calculate the magnitudes of force between the bodies. Typically this
// calculation would probably be done by using a division by the cube of
// the distance (or similarly a multiplication by the cube of its
// reciprocal) but for better performance on modern computers it often
// will make sense to do part of the calculation using a division by the
// distance_Squared which was already calculated earlier. Additionally
// this method is probably a little more accurate due to less rounding
// as well.
magnitudes.as_vectors()[i] = unsafe {
_mm_mul_pd(
_mm_div_pd(_mm_set1_pd(0.01), distance_Squared),
distance_Reciprocal,
)
};
}
// Use the calculated magnitudes of force to update the velocities for all
// of the bodies.
// Apply every other body's gravitation to each body's velocity
{
let mut k = 0;
for i in 0..BODIES_COUNT - 1 {
for j in i + 1..BODIES_COUNT {
let i_mass_magnitude = bodies[i].mass * magnitudes.as_scalars()[k];
let j_mass_magnitude = bodies[j].mass * magnitudes.as_scalars()[k];
for m in 0..3 {
bodies[i].velocity[m] -= position_Deltas[m].as_scalars()[k] * j_mass_magnitude;
bodies[j].velocity[m] += position_Deltas[m].as_scalars()[k] * i_mass_magnitude;
for i in 0..BODIES_COUNT-1 {
for j in i+1..BODIES_COUNT {
let i_mass_mag = bodies[i].mass * magnitudes[k];
let j_mass_mag = bodies[j].mass * magnitudes[k];
for (m, pd) in position_deltas[k].iter().enumerate() {
bodies[i].velocity[m] -= *pd * j_mass_mag;
bodies[j].velocity[m] += *pd * i_mass_mag;
}
k += 1;
}
}
}
// Use the updated velocities to update the positions for all of the bodies.
for i in 0..BODIES_COUNT {
for m in 0..3 {
bodies[i].position[m] += 0.01 * bodies[i].velocity[m];
// Update each body's position.
for body in bodies {
for (m, pos) in body.position.iter_mut().enumerate() {
*pos += 0.01 * body.velocity[m];
}
}
}
fn main() {
let mut solar_Bodies = STARTING_STATE;
let mut position_Deltas: [Interactions; 3] = [Interactions {
scalars: [0.; ROUNDED_INTERACTIONS_COUNT],
}; 3];
let mut magnitudes: Interactions = Interactions {
scalars: [0.; ROUNDED_INTERACTIONS_COUNT],
};
offset_Momentum(&mut solar_Bodies);
output_Energy(&mut solar_Bodies);
let c = std::env::args().nth(1).unwrap().parse().unwrap();
let mut solar_bodies = STARTING_STATE;
offset_momentum(&mut solar_bodies);
output_energy(&mut solar_bodies);
for _ in 0..c {
advance(&mut solar_Bodies, &mut position_Deltas, &mut magnitudes);
advance(&mut solar_bodies)
}
output_Energy(&mut solar_Bodies);
output_energy(&mut solar_bodies);
}