diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h index c88cbada0c84271fc38b0c8fefe6d3f3c305012e..987e20c92d1cdd3d2e9b7b937e7ca900b1765df3 100644 --- a/src/kernel_long_gravity.h +++ b/src/kernel_long_gravity.h @@ -30,7 +30,7 @@ /* Standard headers */ #include <math.h> -#define GADGET2_LONG_RANGE_CORRECTION +//#define GADGET2_LONG_RANGE_CORRECTION /** * @brief Computes the long-range correction term for the potential calculation @@ -51,7 +51,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval( #else const float x = 2.f * u; - const float exp_x = good_approx_expf(x); + const float exp_x = expf(x); //good_approx_expf(x); const float alpha = 1.f / (1.f + exp_x); /* We want 2 - 2 exp(x) * alpha */ @@ -84,12 +84,22 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval( *W = term1 + term2; #else - const float arg = 2.f * u; - const float exp_arg = good_approx_expf(arg); - const float term = 1.f / (1.f + exp_arg); + const float x = 2.f * u; + const float exp_x = expf(x); //good_approx_expf(x); + const float alpha = 1.f / (1.f + exp_x); - *W = arg * exp_arg * term * term - exp_arg * term + 1.f; + /* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */ + *W = 1.f - alpha; + *W = *W * x - exp_x; + *W = *W * alpha + 1.f; *W *= 2.f; + + /* const float arg = 2.f * u; */ + /* const float exp_arg = good_approx_expf(arg); */ + /* const float term = 1.f / (1.f + exp_arg); */ + + /* *W = arg * exp_arg * term * term - exp_arg * term + 1.f; */ + /* *W *= 2.f; */ #endif }