From 36b2be3b95e2de9c6e868bd6c6fa7452bdeabf55 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Fri, 24 Apr 2020 21:24:45 +0200
Subject: [PATCH] Avoid FPE in the M2L MAC
---
src/multipole_accept.h | 10 ++++++----
1 file changed, 6 insertions(+), 4 deletions(-)
diff --git a/src/multipole_accept.h b/src/multipole_accept.h
index 899e1c6412..cf78bb3bb0 100644
--- a/src/multipole_accept.h
+++ b/src/multipole_accept.h
@@ -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);
}
E_BA_term *= 8.f;
- E_BA_term *= max(rho_A, rho_B);
- E_BA_term /= (rho_A + rho_B);
+ if (rho_A != 0. && rho_B != 0.) {
+ 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
const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
#else
@@ -163,7 +165,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
/* 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) */
#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
--
GitLab