diff --git a/src/main.rs b/src/main.rs index aecc9f3..fe834c0 100644 --- a/src/main.rs +++ b/src/main.rs @@ -88,6 +88,51 @@ static mut solar_Bodies: [body; BODIES_COUNT] = [ }, ]; +// 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; + +// 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. @@ -136,68 +181,19 @@ fn output_Energy(bodies: &mut [body; BODIES_COUNT]) { // 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. -unsafe fn advance(bodies: &mut [body; BODIES_COUNT]) { - // 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; - - // 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 - } - } - } - - static mut position_Deltas: [Interactions; 3] = - [Interactions { scalars: [0.; ROUNDED_INTERACTIONS_COUNT] }; 3]; - static mut magnitudes: Interactions = - Interactions { scalars: [0.; ROUNDED_INTERACTIONS_COUNT] }; - +unsafe 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. { 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]; + position_Deltas[m].as_scalars()[k] = + bodies[i].position[m] - bodies[j].position[m]; } k += 1; } @@ -289,12 +285,20 @@ unsafe fn advance(bodies: &mut [body; BODIES_COUNT]) { } fn main() { + // These are new: + let mut position_Deltas: [Interactions; 3] = [Interactions { + scalars: [0.; ROUNDED_INTERACTIONS_COUNT], + }; 3]; + let mut magnitudes: Interactions = Interactions { + scalars: [0.; ROUNDED_INTERACTIONS_COUNT], + }; + unsafe { offset_Momentum(&mut solar_Bodies); output_Energy(&mut solar_Bodies); let c = std::env::args().nth(1).unwrap().parse().unwrap(); for _ in 0..c { - advance(&mut solar_Bodies); + advance(&mut solar_Bodies, &mut position_Deltas, &mut magnitudes); } output_Energy(&mut solar_Bodies); }