diff --git a/qpms/ewald.c b/qpms/ewald.c
index d8915e6..bcac6fd 100644
--- a/qpms/ewald.c
+++ b/qpms/ewald.c
@@ -139,6 +139,7 @@ int ewald32_sigma0(complex double *result, double *err,
 {
   qpms_csf_result gam;
   int retval = complex_gamma_inc_e(-0.5, -sq(k/(2*eta)), &gam);
+  gam.val = conj(gam.val); // We take the other branch, cf. [Linton, p. 642 in the middle]
   if (0 != retval)
     abort();
   *result = gam.val * c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI;
@@ -191,6 +192,9 @@ int ewald32_sigma_long_shiftedpoints (
     complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné
     for(qpms_l_t j = 0; j <= lMax/2; ++j) {
       int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j);
+      // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above
+      if(creal(z) < 0) 
+        Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above
       if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort();
     }
     // R-DEPENDENT END
diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c
index 0f71725..20879e5 100644
--- a/qpms/ewaldsf.c
+++ b/qpms/ewaldsf.c
@@ -35,7 +35,7 @@ void IgnoreUnderflowsGSLErrorHandler (const char * reason,
 
 // DLMF 8.7.3 (latter expression) for complex second argument
 // BTW if a is large negative, it might take a while to evaluate.
-int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result) {
+int cx_gamma_inc_series_e(const double a, const complex double z, qpms_csf_result * result) {
   if (a <= 0 && a == (int) a) {
     result->val = NAN + NAN*I;
     result->err = NAN;
@@ -62,7 +62,7 @@ int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result) {
       result->val = NAN + NAN*I; result->err = NAN;
       return retval;
     }
-    complex double summand = - cpow(z, k) / fullgamma_ak.val;
+    complex double summand = - cpow(z, k) / fullgamma_ak.val; // TODO test branch selection here with cimag(z) = -0.0
     ckahanadd(&sum, &sumc, summand);
     double summanderr = fabs(fullgamma_ak.err * cabs(summand / fullgamma_ak.val));
     // TODO add also the rounding error
diff --git a/tests/ewalds.c b/tests/ewalds.c
index 9d7908f..bfd3e2e 100644
--- a/tests/ewalds.c
+++ b/tests/ewalds.c
@@ -52,17 +52,34 @@ ewaldtest_triang_params paramslist[] = {
   { 6, {1.1, 0.23}, 2.3, 0.97, 2.3,  100, 100, 1., TRIANGULAR_VERTICAL},
   { 6, {1.1, 0.23}, 2.3, 0.97, 2.9,  100, 100, 1., TRIANGULAR_VERTICAL},
 #endif
-  { 2, {1.1, 0.23}, 2.3, 0.97, 0.5,  20, 20, 1., TRIANGULAR_VERTICAL},
-  { 2, {-1.1, -0.23}, 2.3, 0.97, 0.5,  20, 20, 1., TRIANGULAR_VERTICAL},
-  { 2, {0, 1.1}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
-  { 2, {0, -1.1}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
-  { 2, {-1.1, 0}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
-  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
-  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 10, 10, 1., TRIANGULAR_VERTICAL},
-  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 5, 5, 1., TRIANGULAR_VERTICAL},
-  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL},
-  { 2, {0.3, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL},
-  { 2, {2.7, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0.23}, 2.3, 0.97, 0.5,  20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {-1.1, -0.23}, 2.3, 0.97, 0.5,  20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {0, 1.1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {0, -1.1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {-1.1, 0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 5.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 5.5, 10, 80, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 5.5, 5, 40, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 5.5, 2, 16, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 2.5, 10, 80, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 2.5, 5, 40, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 2.5, 2, 16, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 10, 80, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 5, 40, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 2, 16, 1., TRIANGULAR_VERTICAL},
+  { 2, {0.3, 0}, 2.3, 0.97, 0.5, 2, 16, 1., TRIANGULAR_VERTICAL},
+  { 2, {2.7, 1}, 2.3, 0.97, 0.5, 5, 40, 1., TRIANGULAR_VERTICAL},
+  { 2, {2.7, 1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {2.7, 1}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {2.7, 1}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {2.7, 1}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 1}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 1}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 1}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
+
 
 // end:
 //  { 0,  {0, 0}, 0, 0, 0, 0, 0, 0, 0}
@@ -91,6 +108,10 @@ void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) {
 }
 
 
+static inline double san(double x) {
+  return fabs(x) < 1e-13 ? 0 : x;
+}
+
 ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p);
 
 int main() {
@@ -110,15 +131,15 @@ int main() {
         qpms_y_t y_conj = qpms_mn2y_sc(-m,n);
         // y n m sigma_total (err), regsigmas_416 regsigmas_415_recon
         printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n| predict %.16g%+.16gj \n| actual  %.16g%+.16gj\n",
-            y, n, m, creal(r->sigmas_total[y]), cimag(r->sigmas_total[y]),
+            y, n, m, creal(san(r->sigmas_total[y])), san(cimag(r->sigmas_total[y])),
             r->err_sigmas_total[y],
-            creal(r->sigmas_long[y]), cimag(r->sigmas_long[y]),
+            san(creal(r->sigmas_long[y])), san(cimag(r->sigmas_long[y])),
             r->err_sigmas_long[y],
-            creal(r->sigmas_short[y]), cimag(r->sigmas_short[y]),
+            san(creal(r->sigmas_short[y])), san(cimag(r->sigmas_short[y])),
             r->err_sigmas_short[y],
-            creal(r->regsigmas_416[y]), cimag(r->regsigmas_416[y]),
-            creal(r->sigmas_total[y])+creal(r->sigmas_total[y_conj]),
-            cimag(r->sigmas_total[y])-cimag(r->sigmas_total[y_conj])
+            san(creal(r->regsigmas_416[y])), san(cimag(r->regsigmas_416[y])),
+            san(creal(r->sigmas_total[y]) + creal(r->sigmas_total[y_conj])),
+            san(cimag(r->sigmas_total[y]) - cimag(r->sigmas_total[y_conj]))
         );
       }
     }