Skip to content
Snippets Groups Projects
Commit 36b2be3b authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Avoid FPE in the M2L MAC

parent 4a9c4642
No related branches found
No related tags found
1 merge request!1077Improved multipole acceptance criterion (MAC)
Hide whitespace changes
Inline Side-by-side
Showing
with 6 additions and 4 deletions
src/multipole_accept.h
+ 6
− 4
View file @ 36b2be3b
...@@ -63,10 +63,12 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept( ...@@ -63,10 +63,12 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
binomial(p, n) * B->m_pole.power[n] * integer_powf(rho_A, p - n); binomial(p, n) * B->m_pole.power[n] * integer_powf(rho_A, p - n);
} }
E_BA_term *= 8.f; E_BA_term *= 8.f;
E_BA_term *= max(rho_A, rho_B); if (rho_A != 0. && rho_B != 0.) {
E_BA_term /= (rho_A + rho_B); E_BA_term *= max(rho_A, rho_B);
E_BA_term /= (rho_A + rho_B);
}
/* Compute r^(p+2) */ /* Compute r^(p+2) */
#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1 #if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2)); const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
#else #else
...@@ -163,7 +165,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept( ...@@ -163,7 +165,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
const int p = SELF_GRAVITY_MULTIPOLE_ORDER; const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
/* Compute the error estimator (without the 1/M_B term that cancels out) */ /* Compute the error estimator (without the 1/M_B term that cancels out) */
float E_BA_term = 8.f * B->m_pole.power[p]; const float E_BA_term = 8.f * B->m_pole.power[p];
/* Compute r^(p+2) */ /* Compute r^(p+2) */
#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1 #if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment