fix integer division bug

andrewfowlie · Oct 31, 2019 · 8d6505b · 8d6505b
1 parent c0b3e4a
commit 8d6505b
Showing 1 changed file with 3 additions and 1 deletion.
diff --git a/src/zeta.cpp b/src/zeta.cpp
@@ -12,6 +12,7 @@
 #include <bernoulli.h>
 
 #include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
 #include <complex>
 #include <stdexcept>
 
@@ -109,7 +110,8 @@ cdouble polylog(double s, cdouble a, int N) {
     cdouble term = a;
     cdouble sum = term;
     for (int i = 2; i <= N; i += 1) {
-      term *= a * pow((i - 1) / i, 2.5);
+      const double x = (i - 1) / static_cast<double>(i);
+      term *= a * gsl_pow_2(x) * sqrt(x);
       sum += term;
     }
     return sum;