diff --git a/qpms/ewald.c b/qpms/ewald.c
index 26bda5d..7720a33 100644
--- a/qpms/ewald.c
+++ b/qpms/ewald.c
@@ -189,6 +189,17 @@ int ewald3_sigma0(complex double *result, double *err,
   return 0;
 }
 
+// Wrapper over cpow to correctly handle the k->0 behaviour
+static inline complex double cpow_0lim_zi(complex double x, long z)
+{
+  if (x == 0 && z == 0)
+    return 1;
+  else if (x == 0 && z > 0) // lol branch cut in cpow
+    return 0;
+  else
+    return cpow(x, z);
+}
+
 int ewald3_21_xy_sigma_long (
     complex double *target, // must be c->nelem_sc long
     double *err,
@@ -293,7 +304,7 @@ int ewald3_21_xy_sigma_long (
         double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
         assert((n-abs(m))/2 == c->s1_jMaxes[y]);
         for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
-          complex double summand = cpow(rbeta_pq/k, n-2*j) 
+          complex double summand = cpow_0lim_zi(rbeta_pq/k, n-2*j)
             * e_imalpha_pq  * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop
             * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation
             * c->s1_constfacs[y][j];
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 6e6797b..3d1ef0a 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -48,3 +48,8 @@ add_test( NAME scalar_ewald32_square_cplxk_negative_Im2 COMMAND test_scalar_ewal
 #	lMax b1.x b1.y b2.x b2.y wavenum.real wavenum.imag k.x k.y particle_shift.x particle_shift.y csphase rtol atol maxR maxK eta1 ...
 	3    1    0    0    1    2.3          -0.1          2.7   1              0.0           0.0        -1 1e-8 1e-10  20  160  1.5   1.6 2.5 2.6
 	)
+
+add_test( NAME scalar_ewald32_square_zerok COMMAND test_scalar_ewald32 
+#	lMax b1.x b1.y b2.x b2.y wavenum.real wavenum.imag k.x k.y particle_shift.x particle_shift.y csphase rtol atol maxR maxK eta1 ...
+	3    1    0    0    1    2.3          0            0   0                0.5           0.1325      -1 1e-8 1e-10  40  160  1.5   1.6 2.5 2.6
+	)
diff --git a/tests/test_scalar_ewald32.c b/tests/test_scalar_ewald32.c
index 7673876..679f663 100644
--- a/tests/test_scalar_ewald32.c
+++ b/tests/test_scalar_ewald32.c
@@ -17,6 +17,8 @@
 #include <string.h>
 #include <float.h>
 #include <gsl/gsl_sf_legendre.h>
+#include <fenv.h>
+
 typedef struct ewaldtest2d_params {
   qpms_l_t lMax;
   point2d b1, b2;
@@ -66,6 +68,7 @@ int isclose_cmplx(complex double a, complex double b, double rtol, double atol)
 ewaldtest2d_results *ewaldtest2d(const ewaldtest2d_params p);
 
 int main(int argc, char **argv) {
+  feenableexcept(FE_INVALID | FE_OVERFLOW);
   bool verbose = !!getenv("QPMS_VERBOSE_TESTS");
   gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler);
   QPMS_ENSURE(argc >= 18, "At least 16 arguments expected, I see only %d.", argc-1);