Minimize unsafe in advance, make main safe

This commit is contained in:
Timothy Warren 2020-01-13 13:02:30 -05:00
parent fd9ac09655
commit a61113711f
1 changed files with 36 additions and 31 deletions

View File

@ -23,7 +23,7 @@ const SOLAR_MASS: f64 = 4. * PI * PI;
const DAYS_PER_YEAR: f64 = 365.24; const DAYS_PER_YEAR: f64 = 365.24;
const BODIES_COUNT: usize = 5; const BODIES_COUNT: usize = 5;
static mut solar_Bodies: [body; BODIES_COUNT] = [ const STARTING_STATE: [body; BODIES_COUNT] = [
body { body {
// Sun // Sun
mass: SOLAR_MASS, mass: SOLAR_MASS,
@ -181,7 +181,8 @@ fn output_Energy(bodies: &mut [body; BODIES_COUNT]) {
// interactions between all the bodies, update each body's velocity based on // interactions between all the bodies, update each body's velocity based on
// those interactions, and update each body's position by the distance it // those interactions, and update each body's position by the distance it
// travels in a timestep at it's updated velocity. // travels in a timestep at it's updated velocity.
unsafe fn advance( #[cfg(target_feature = "sse2")]
fn advance(
bodies: &mut [body; BODIES_COUNT], bodies: &mut [body; BODIES_COUNT],
position_Deltas: &mut [Interactions; 3], position_Deltas: &mut [Interactions; 3],
magnitudes: &mut Interactions, magnitudes: &mut Interactions,
@ -205,19 +206,21 @@ unsafe fn advance(
// ROUNDED_INTERACTIONS_COUNT/2 iterations are done. // ROUNDED_INTERACTIONS_COUNT/2 iterations are done.
for i in 0..ROUNDED_INTERACTIONS_COUNT / 2 { for i in 0..ROUNDED_INTERACTIONS_COUNT / 2 {
// Load position_Deltas of two bodies into position_Delta. // Load position_Deltas of two bodies into position_Delta.
let mut position_Delta = [_mm_setzero_pd(); 3]; let mut position_Delta = [unsafe { _mm_setzero_pd() }; 3];
for m in 0..3 { for m in 0..3 {
position_Delta[m] = position_Deltas[m].as_vectors()[i]; position_Delta[m] = position_Deltas[m].as_vectors()[i];
} }
let distance_Squared: __m128d = _mm_add_pd( let distance_Squared: __m128d = unsafe {
_mm_add_pd( _mm_add_pd(
_mm_mul_pd(position_Delta[0], position_Delta[0]), _mm_add_pd(
_mm_mul_pd(position_Delta[1], position_Delta[1]), _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]), ),
); _mm_mul_pd(position_Delta[2], position_Delta[2]),
)
};
// Doing square roots normally using double precision floating point // Doing square roots normally using double precision floating point
// math can be quite time consuming so SSE's much faster single // math can be quite time consuming so SSE's much faster single
@ -226,23 +229,25 @@ unsafe fn advance(
// acceptable results so two iterations of the Newton–Raphson method are // acceptable results so two iterations of the Newton–Raphson method are
// done to improve precision further. // done to improve precision further.
let mut distance_Reciprocal: __m128d = let mut distance_Reciprocal: __m128d =
_mm_cvtps_pd(_mm_rsqrt_ps(_mm_cvtpd_ps(distance_Squared))); unsafe { _mm_cvtps_pd(_mm_rsqrt_ps(_mm_cvtpd_ps(distance_Squared))) };
for _ in 0..2 { for _ in 0..2 {
// Normally the last four multiplications in this equation would // Normally the last four multiplications in this equation would
// have to be done sequentially but by placing the last // have to be done sequentially but by placing the last
// multiplication in parentheses, a compiler can then schedule that // multiplication in parentheses, a compiler can then schedule that
// multiplication earlier. // multiplication earlier.
distance_Reciprocal = _mm_sub_pd( distance_Reciprocal = unsafe {
_mm_mul_pd(distance_Reciprocal, _mm_set1_pd(1.5)), _mm_sub_pd(
_mm_mul_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), _mm_mul_pd(
distance_Reciprocal, _mm_mul_pd(_mm_set1_pd(0.5), distance_Squared),
distance_Reciprocal,
),
_mm_mul_pd(distance_Reciprocal, distance_Reciprocal),
), ),
_mm_mul_pd(distance_Reciprocal, distance_Reciprocal), )
), };
);
} }
// Calculate the magnitudes of force between the bodies. Typically this // Calculate the magnitudes of force between the bodies. Typically this
@ -253,10 +258,12 @@ unsafe fn advance(
// distance_Squared which was already calculated earlier. Additionally // distance_Squared which was already calculated earlier. Additionally
// this method is probably a little more accurate due to less rounding // this method is probably a little more accurate due to less rounding
// as well. // as well.
magnitudes.as_vectors()[i] = _mm_mul_pd( magnitudes.as_vectors()[i] = unsafe {
_mm_div_pd(_mm_set1_pd(0.01), distance_Squared), _mm_mul_pd(
distance_Reciprocal, _mm_div_pd(_mm_set1_pd(0.01), distance_Squared),
); distance_Reciprocal,
)
};
} }
// Use the calculated magnitudes of force to update the velocities for all // Use the calculated magnitudes of force to update the velocities for all
@ -285,7 +292,7 @@ unsafe fn advance(
} }
fn main() { fn main() {
// These are new: let mut solar_Bodies = STARTING_STATE;
let mut position_Deltas: [Interactions; 3] = [Interactions { let mut position_Deltas: [Interactions; 3] = [Interactions {
scalars: [0.; ROUNDED_INTERACTIONS_COUNT], scalars: [0.; ROUNDED_INTERACTIONS_COUNT],
}; 3]; }; 3];
@ -293,13 +300,11 @@ fn main() {
scalars: [0.; ROUNDED_INTERACTIONS_COUNT], scalars: [0.; ROUNDED_INTERACTIONS_COUNT],
}; };
unsafe { offset_Momentum(&mut solar_Bodies);
offset_Momentum(&mut solar_Bodies); output_Energy(&mut solar_Bodies);
output_Energy(&mut solar_Bodies); let c = std::env::args().nth(1).unwrap().parse().unwrap();
let c = std::env::args().nth(1).unwrap().parse().unwrap(); for _ in 0..c {
for _ in 0..c { advance(&mut solar_Bodies, &mut position_Deltas, &mut magnitudes);
advance(&mut solar_Bodies, &mut position_Deltas, &mut magnitudes);
}
output_Energy(&mut solar_Bodies);
} }
output_Energy(&mut solar_Bodies);
} }