diff --git a/src/c/CMakeLists.txt b/src/c/CMakeLists.txt
index 54b97ef..58c7f98 100644
--- a/src/c/CMakeLists.txt
+++ b/src/c/CMakeLists.txt
@@ -9,6 +9,7 @@ add_library(polylogarithm_c
   Li4.c
   Li5.c
   Li6.c
+  log.c
 )
 
 target_include_directories(polylogarithm_c PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
diff --git a/src/c/Li2.c b/src/c/Li2.c
index d7a8652..1e60dd8 100644
--- a/src/c/Li2.c
+++ b/src/c/Li2.c
@@ -7,7 +7,7 @@
 #include <complex.h>
 #include <float.h>
 #include <math.h>
-#include "fast_clog.h"
+#include "log.h"
 
 
 static long double hornerl(long double x, const long double* c, int len)
@@ -331,23 +331,23 @@ float _Complex cli2f(float _Complex z)
    /* transformation to |z|<1, Re(z)<=0.5f */
    if (rz <= 0.5f) {
       if (nz > 1.0f) {
-         const float _Complex lz = fast_clogf(-z);
-         u = -fast_clogf(1.0f - 1.0f / z);
+         const float _Complex lz = clogf(-z);
+         u = -clogf(1.0f - 1.0f / z);
          rest = -0.5f*lz*lz - PI*PI/6;
          sgn = -1;
       } else { /* nz <= 1 */
-         u = -fast_clogf(1.0f - z);
+         u = -clogf(1.0f - z);
          rest = 0;
          sgn = 1;
       }
    } else { /* rz > 0.5f */
       if (nz <= 2*rz) {
-         u = -fast_clogf(z);
-         rest = u*fast_clogf(1.0f - z) + PI*PI/6;
+         u = -clogf(z);
+         rest = u*clogf(1.0f - z) + PI*PI/6;
          sgn = -1;
       } else { /* nz > 2*rz */
-         const float _Complex lz = fast_clogf(-z);
-         u = -fast_clogf(1.0f - 1.0f / z);
+         const float _Complex lz = clogf(-z);
+         u = -clogf(1.0f - 1.0f / z);
          rest = -0.5f*lz*lz - PI*PI/6;
          sgn = -1;
       }
@@ -412,23 +412,23 @@ double _Complex cli2(double _Complex z)
    /* transformation to |z|<1, Re(z)<=0.5 */
    if (rz <= 0.5) {
       if (nz > 1.0) {
-         const double _Complex lz = fast_clog(-z);
-         u = -fast_clog(1.0 - 1.0 / z);
+         const double _Complex lz = clog(-z);
+         u = -clog1p(-1.0/z);
          rest = -0.5*lz*lz - PI*PI/6;
          sgn = -1;
       } else { /* nz <= 1 */
-         u = -fast_clog(1.0 - z);
+         u = -clog1p(-z);
          rest = 0;
          sgn = 1;
       }
    } else { /* rz > 0.5 */
       if (nz <= 2*rz) {
-         u = -fast_clog(z);
-         rest = u*fast_clog(1.0 - z) + PI*PI/6;
+         u = -clog(z);
+         rest = u*clog1p(-z) + PI*PI/6;
          sgn = -1;
       } else { /* nz > 2*rz */
-         const double _Complex lz = fast_clog(-z);
-         u = -fast_clog(1.0 - 1.0 / z);
+         const double _Complex lz = clog(-z);
+         u = -clog1p(-1.0/z);
          rest = -0.5*lz*lz - PI*PI/6;
          sgn = -1;
       }
@@ -517,23 +517,23 @@ long double _Complex cli2l(long double _Complex z)
    /* transformation to |z|<1, Re(z)<=0.5 */
    if (rz <= 0.5L) {
       if (nz > 1.0L) {
-         const long double _Complex lz = fast_clogl(-z);
-         u = -fast_clogl(1.0L - 1.0L/z);
+         const long double _Complex lz = clogl(-z);
+         u = -clog1pl(-1.0L/z);
          rest = -0.5L*lz*lz - PI*PI/6;
          sgn = -1;
       } else { /* nz <= 1 */
-         u = -fast_clogl(1.0L - z);
+         u = -clog1pl(-z);
          rest = 0;
          sgn = 1;
       }
    } else { /* rz > 0.5L */
       if (nz <= 2*rz) {
-         u = -fast_clogl(z);
-         rest = u * fast_clogl(1.0L - z) + PI*PI/6;
+         u = -clogl(z);
+         rest = u * clog1pl(-z) + PI*PI/6;
          sgn = -1;
       } else { /* nz > 2*rz */
-         const long double _Complex lz = fast_clogl(-z);
-         u = -fast_clogl(1.0L - 1.0L/z);
+         const long double _Complex lz = clogl(-z);
+         u = -clog1pl(-1.0L/z);
          rest = -0.5L*lz*lz - PI*PI/6;
          sgn = -1;
       }
diff --git a/src/c/Li3.c b/src/c/Li3.c
index ff57f9f..ef1f532 100644
--- a/src/c/Li3.c
+++ b/src/c/Li3.c
@@ -7,7 +7,7 @@
 #include <complex.h>
 #include <float.h>
 #include <math.h>
-#include "fast_clog.h"
+#include "log.h"
 
 /// Li_3(x) for x in [-1,0]
 static double li3_neg(double x)
@@ -145,7 +145,7 @@ double _Complex cli3(double _Complex z)
       const double _Complex u4 = u2*u2;
       const double _Complex u8 = u4*u4;
       const double _Complex c0 = zeta3 + u*(zeta2 - u2/12.0);
-      const double _Complex c1 = 0.25 * (3.0 - 2.0*fast_pos_clog(-u));
+      const double _Complex c1 = 0.25 * (3.0 - 2.0*pos_clog(-u));
 
       const double cs[7] = {
          -3.4722222222222222e-03, 1.1574074074074074e-05,
@@ -165,11 +165,11 @@ double _Complex cli3(double _Complex z)
    double _Complex u = 0.0, rest = 0.0;
 
    if (nz <= 1.0) {
-      u = -fast_pos_clog(1.0 - z);
+      u = -pos_clog(1.0 - z);
    } else { // nz > 1
       const double arg = pz > 0.0 ? pz - PI : pz + PI;
       const double _Complex lmz = lnz + arg*I; // clog(-z)
-      u = -fast_pos_clog(1.0 - 1.0/z);
+      u = -pos_clog(1.0 - 1.0/z);
       rest = -lmz*(lmz*lmz/6.0 + zeta2);
    }
 
@@ -275,7 +275,7 @@ long double _Complex cli3l(long double _Complex z)
       const long double _Complex u  = lnz + pz*I; // clog(z)
       const long double _Complex u2 = u*u;
       const long double _Complex c0 = zeta3 + u*(zeta2 - u2/12.0L);
-      const long double _Complex c1 = 0.25L * (3.0L - 2.0L*fast_pos_clogl(-u));
+      const long double _Complex c1 = 0.25L * (3.0L - 2.0L*pos_clogl(-u));
 
       const long double cs[] = {
         -3.47222222222222222222222222222222222e-03L,
@@ -317,11 +317,11 @@ long double _Complex cli3l(long double _Complex z)
    long double _Complex u = 0.0L, rest = 0.0L;
 
    if (nz <= 1.0L) {
-      u = -fast_pos_clogl(1.0L - z);
+      u = -pos_clogl(1.0L - z);
    } else { // nz > 1
       const long double arg = pz > 0.0L ? pz - PI : pz + PI;
       const long double _Complex lmz = lnz + arg*I; // clog(-z)
-      u = -fast_pos_clogl(1.0L - 1.0L/z);
+      u = -pos_clogl(1.0L - 1.0L/z);
       rest = -lmz*(lmz*lmz/6.0L + zeta2);
    }
 
diff --git a/src/c/Li4.c b/src/c/Li4.c
index ce42e8d..0c7cf7a 100644
--- a/src/c/Li4.c
+++ b/src/c/Li4.c
@@ -7,7 +7,7 @@
 #include <complex.h>
 #include <float.h>
 #include <math.h>
-#include "fast_clog.h"
+#include "log.h"
 
 /// Li_4(x) for x in [-1,0]
 static double li4_neg(double x)
@@ -200,7 +200,7 @@ double _Complex cli4(double _Complex z)
       const double _Complex u8 = u4*u4;
       const double c1 = 1.2020569031595943; // zeta(3)
       const double c2 = 0.82246703342411322;
-      const double _Complex c3 = (11.0/6.0 - fast_pos_clog(-u))/6.0;
+      const double _Complex c3 = (11.0/6.0 - pos_clog(-u))/6.0;
       const double c4 = -1.0/48.0;
 
       const double cs[7] = {
@@ -224,12 +224,12 @@ double _Complex cli4(double _Complex z)
    double sgn = 1;
 
    if (nz <= 1.0) {
-      u = -fast_pos_clog(1.0 - z);
+      u = -pos_clog(1.0 - z);
    } else { // nz > 1
       const double arg = pz > 0.0 ? pz - PI : pz + PI;
       const double _Complex lmz = lnz + arg*I; // clog(-z)
       const double _Complex lmz2 = lmz*lmz;
-      u = -fast_pos_clog(1.0 - 1.0/z);
+      u = -pos_clog(1.0 - 1.0/z);
       r = 1.0/360.0*(-7*PI4 + lmz2*(-30.0*PI2 - 15.0*lmz2));
       sgn = -1;
    }
@@ -336,7 +336,7 @@ long double _Complex cli4l(long double _Complex z)
       const long double _Complex u2 = u*u;
       const long double c1 = 1.20205690315959428539973816151144999L; // zeta(3)
       const long double c2 = 0.822467033424113218236207583323012595L;
-      const long double _Complex c3 = (11.0L/6.0L - fast_pos_clogl(-u))/6.0L;
+      const long double _Complex c3 = (11.0L/6.0L - pos_clogl(-u))/6.0L;
       const long double c4 = -1.0L/48.0L;
 
       const long double cs[] = {
@@ -377,12 +377,12 @@ long double _Complex cli4l(long double _Complex z)
    long double sgn = 1;
 
    if (nz <= 1.0L) {
-      u = -fast_pos_clogl(1.0L - z);
+      u = -pos_clogl(1.0L - z);
    } else { // nz > 1
       const long double arg = pz > 0.0 ? pz - PI : pz + PI;
       const long double _Complex lmz = lnz + arg*I; // clog(-z)
       const long double _Complex lmz2 = lmz*lmz;
-      u = -fast_pos_clogl(1.0L - 1.0L/z);
+      u = -pos_clogl(1.0L - 1.0L/z);
       r = 1.0L/360.0L*(-7*PI4 + lmz2*(-30.0L*PI2 - 15.0L*lmz2));
       sgn = -1;
    }
diff --git a/src/c/Li5.c b/src/c/Li5.c
index 1a09c05..2249583 100644
--- a/src/c/Li5.c
+++ b/src/c/Li5.c
@@ -7,7 +7,7 @@
 #include <complex.h>
 #include <float.h>
 #include <math.h>
-#include "fast_clog.h"
+#include "log.h"
 
 
 /**
@@ -61,7 +61,7 @@ double _Complex cli5(double _Complex z)
       const double c1 = 1.0823232337111382; // zeta(4)
       const double c2 = 0.60102845157979714; // zeta(3)/2
       const double c3 = 0.27415567780803774;
-      const double _Complex c4 = (25.0/12.0 - fast_pos_clog(-u))/24.0;
+      const double _Complex c4 = (25.0/12.0 - pos_clog(-u))/24.0;
       const double c5 = -1.0/240.0;
 
       const double cs[6] = {
@@ -84,12 +84,12 @@ double _Complex cli5(double _Complex z)
    double _Complex u = 0.0, rest = 0.0;
 
    if (nz <= 1.0) {
-      u = -fast_pos_clog(1.0 - z);
+      u = -pos_clog(1.0 - z);
    } else { // nz > 1
       const double arg = pz > 0.0 ? pz - PI : pz + PI;
       const double _Complex lmz = lnz + arg*I; // clog(-z)
       const double _Complex lmz2 = lmz*lmz;
-      u = -fast_pos_clog(1.0 - 1.0/z);
+      u = -pos_clog(1.0 - 1.0/z);
       rest = -1.0/360.0*lmz*(7*PI4 + lmz2*(10.0*PI2 + 3.0*lmz2));
    }
 
@@ -196,7 +196,7 @@ long double _Complex cli5l(long double _Complex z)
       const long double c1 = 1.08232323371113819151600369654116790L; // zeta(4)
       const long double c2 = 0.601028451579797142699869080755724995L; // zeta(3)/2
       const long double c3 = 0.274155677808037739412069194441004198L;
-      const long double _Complex c4 = (25.0L/12.0L - fast_pos_clogl(-u))/24.0L;
+      const long double _Complex c4 = (25.0L/12.0L - pos_clogl(-u))/24.0L;
       const long double c5 = -1.0L/240.0L;
 
       const long double cs[] = {
@@ -236,12 +236,12 @@ long double _Complex cli5l(long double _Complex z)
    long double _Complex u = 0.0L, rest = 0.0L;
 
    if (nz <= 1.0L) {
-      u = -fast_pos_clogl(1.0L - z);
+      u = -pos_clogl(1.0L - z);
    } else { // nz > 1
       const long double arg = pz > 0.0 ? pz - PI : pz + PI;
       const long double _Complex lmz = lnz + arg*I; // clog(-z)
       const long double _Complex lmz2 = lmz*lmz;
-      u = -fast_pos_clogl(1.0L - 1.0L/z);
+      u = -pos_clogl(1.0L - 1.0L/z);
       rest = -1.0L/360.0L*lmz*(7*PI4 + lmz2*(10.0L*PI2 + 3.0L*lmz2));
    }
 
diff --git a/src/c/Li6.c b/src/c/Li6.c
index 3c4b571..555e952 100644
--- a/src/c/Li6.c
+++ b/src/c/Li6.c
@@ -7,7 +7,7 @@
 #include <complex.h>
 #include <float.h>
 #include <math.h>
-#include "fast_clog.h"
+#include "log.h"
 
 
 /**
@@ -62,7 +62,7 @@ double _Complex cli6(double _Complex z)
       const double c2 = 0.54116161685556910;
       const double c3 = 0.20034281719326571;
       const double c4 = 0.068538919452009435;
-      const double _Complex c5 = (137.0/60.0 - fast_pos_clog(-u))/120.0;
+      const double _Complex c5 = (137.0/60.0 - pos_clog(-u))/120.0;
       const double c6 = -1.0/1440.0;
 
       const double cs[5] = {
@@ -86,12 +86,12 @@ double _Complex cli6(double _Complex z)
    double sgn = 1;
 
    if (nz <= 1.0) {
-      u = -fast_pos_clog(1.0 - z);
+      u = -pos_clog(1.0 - z);
    } else { // nz > 1
       const double arg = pz > 0.0 ? pz - PI : pz + PI;
       const double _Complex lmz = lnz + arg*I; // clog(-z)
       const double _Complex lmz2 = lmz*lmz;
-      u = -fast_pos_clog(1.0 - 1.0/z);
+      u = -pos_clog(1.0 - 1.0/z);
       r = -31.0*PI6/15120.0 + lmz2*(-7.0/720.0*PI4 + lmz2*(-1.0/144.0*PI2 - 1.0/720.0*lmz2));
       sgn = -1;
    }
@@ -202,7 +202,7 @@ long double _Complex cli6l(long double _Complex z)
       const long double c2 = 0.541161616855569095758001848270583951L;
       const long double c3 = 0.200342817193265714233289693585241665L;
       const long double c4 = 0.0685389194520094348530172986102510496L;
-      const long double _Complex c5 = (137.0L/60.0L - fast_pos_clogl(-u))/120.0L;
+      const long double _Complex c5 = (137.0L/60.0L - pos_clogl(-u))/120.0L;
       const long double c6 = -1.0L/1440.0L;
 
       const long double cs[] = {
@@ -247,12 +247,12 @@ long double _Complex cli6l(long double _Complex z)
    long double sgn = 1;
 
    if (nz <= 1.0L) {
-      u = -fast_pos_clogl(1.0L - z);
+      u = -pos_clogl(1.0L - z);
    } else { // nz > 1
       const long double arg = pz > 0.0 ? pz - PI : pz + PI;
       const long double _Complex lmz = lnz + arg*I; // clog(-z)
       const long double _Complex lmz2 = lmz*lmz;
-      u = -fast_pos_clogl(1.0L - 1.0L/z);
+      u = -pos_clogl(1.0L - 1.0L/z);
       r = -31.0L*PI6/15120.0L
           + lmz2*(-7.0L/720.0L*PI4 +
                   lmz2*(-1.0L/144.0L*PI2 - 1.0L/720.0L*lmz2));
diff --git a/src/c/fast_clog.h b/src/c/fast_clog.h
deleted file mode 100644
index 412045e..0000000
--- a/src/c/fast_clog.h
+++ /dev/null
@@ -1,73 +0,0 @@
-/* ====================================================================
- * This file is part of Polylogarithm.
- *
- * Polylogarithm is licenced under the MIT License.
- * ==================================================================== */
-
-#pragma once
-
-#include <complex.h>
-#include <math.h>
-
-
-static inline float _Complex fast_clogf(float _Complex z)
-{
-   const float rz = crealf(z);
-   const float iz = cimagf(z);
-
-   if (iz == 0.0f && rz > 0.0f) {
-      return logf(rz);
-   } else if (iz == 0.0f) {
-      return logf(-rz) + I*3.14159265f;
-   }
-
-   return logf(hypotf(rz, iz)) + I*atan2f(iz, rz);
-}
-
-
-static inline double _Complex fast_clog(double _Complex z)
-{
-   const double rz = creal(z);
-   const double iz = cimag(z);
-
-   return log(hypot(rz, iz)) + I*atan2(iz, rz);
-}
-
-
-static inline double _Complex fast_pos_clog(double _Complex z)
-{
-   const double rz = creal(z);
-   const double iz = cimag(z);
-
-   if (iz == 0.0 && rz > 0.0) {
-      return log(rz);
-   } else if (iz == 0.0) {
-      return log(-rz) + I*3.1415926535897932;
-   }
-
-   return log(hypot(rz, iz)) + I*atan2(iz, rz);
-}
-
-
-static inline long double _Complex fast_clogl(long double _Complex z)
-{
-   const long double rz = creall(z);
-   const long double iz = cimagl(z);
-
-   return logl(hypotl(rz, iz)) + I*atan2l(iz, rz);
-}
-
-
-static inline long double _Complex fast_pos_clogl(long double _Complex z)
-{
-   const long double rz = creall(z);
-   const long double iz = cimagl(z);
-
-   if (iz == 0.0L && rz > 0.0L) {
-      return logl(rz);
-   } else if (iz == 0.0L) {
-      return logl(-rz) + I*3.14159265358979323846264338327950288L;
-   }
-
-   return logl(hypotl(rz, iz)) + I*atan2l(iz, rz);
-}
diff --git a/src/c/log.c b/src/c/log.c
new file mode 100644
index 0000000..a31aba2
--- /dev/null
+++ b/src/c/log.c
@@ -0,0 +1,70 @@
+/* ====================================================================
+ * This file is part of Polylogarithm.
+ *
+ * Polylogarithm is licenced under the MIT License.
+ * ==================================================================== */
+
+#include <complex.h>
+#include <math.h>
+
+
+double _Complex clog1p(double _Complex z)
+{
+   const double _Complex u = 1.0 + z;
+   const double rz = creal(u);
+   const double iz = cimag(u);
+
+   if (rz == 1.0 && iz == 0.0) {
+      return z;
+   } else if (rz <= 0.0) {
+      return clog(u);
+   }
+
+   return clog(u)*(z/(u - 1.0));
+}
+
+
+long double _Complex clog1pl(long double _Complex z)
+{
+   const long double _Complex u = 1.0L + z;
+   const long double rz = creall(u);
+   const long double iz = cimagl(u);
+
+   if (rz == 1.0L && iz == 0.0L) {
+      return z;
+   } else if (rz <= 0.0L) {
+      return clogl(u);
+   }
+
+   return clogl(u)*(z/(u - 1.0L));
+}
+
+
+double _Complex pos_clog(double _Complex z)
+{
+   const double rz = creal(z);
+   const double iz = cimag(z);
+
+   if (iz == 0.0 && rz > 0.0) {
+      return log(rz);
+   } else if (iz == 0.0) {
+      return log(-rz) + I*3.1415926535897932;
+   }
+
+   return clog(z);
+}
+
+
+long double _Complex pos_clogl(long double _Complex z)
+{
+   const long double rz = creall(z);
+   const long double iz = cimagl(z);
+
+   if (iz == 0.0L && rz > 0.0L) {
+      return logl(rz);
+   } else if (iz == 0.0L) {
+      return logl(-rz) + I*3.14159265358979323846264338327950288L;
+   }
+
+   return clogl(z);
+}
diff --git a/src/c/log.h b/src/c/log.h
new file mode 100644
index 0000000..011e2c3
--- /dev/null
+++ b/src/c/log.h
@@ -0,0 +1,19 @@
+/* ====================================================================
+ * This file is part of Polylogarithm.
+ *
+ * Polylogarithm is licenced under the MIT License.
+ * ==================================================================== */
+
+#pragma once
+
+/** returns clog(1 + z) with double precision */
+double _Complex clog1p(double _Complex z);
+
+/** returns clog(1 + z) with long double precision */
+long double _Complex clog1pl(long double _Complex z);
+
+/** returns clog(z) with double precision, treating -0.0 correctly */
+double _Complex pos_clog(double _Complex z);
+
+/** returns clog(z) with long double precision, treating -0.0L correctly */
+long double _Complex pos_clogl(long double _Complex z);
diff --git a/src/cpp/Li.cpp b/src/cpp/Li.cpp
index f3d023f..773994e 100644
--- a/src/cpp/Li.cpp
+++ b/src/cpp/Li.cpp
@@ -39,7 +39,7 @@ namespace {
       if (std::imag(z) == 0.0 && std::real(z) > 0.0) {
          return { std::log(std::real(z)), 0.0 };
       } else if (std::imag(z) == 0.0) {
-         return { std::log(-std::real(z)), 4*std::atan(1.0) };
+         return { std::log(-std::real(z)), PI };
       }
       return { std::log(std::hypot(std::real(z), std::imag(z))), std::arg(z) };
    }
diff --git a/src/cpp/Li2.cpp b/src/cpp/Li2.cpp
index e8ecc61..57bf7d1 100644
--- a/src/cpp/Li2.cpp
+++ b/src/cpp/Li2.cpp
@@ -5,45 +5,15 @@
 // ====================================================================
 
 #include "Li2.hpp"
-#include "complex.hpp"
+#include "horner.hpp"
+#include "log.hpp"
 #include <cfloat>
 #include <cmath>
+#include <complex>
 #include <limits>
 
 namespace polylogarithm {
 
-namespace {
-
-   template <typename T, int N>
-   T horner(T x, const T (&c)[N]) noexcept
-   {
-      T p = c[N - 1];
-      for (int i = N - 2; i >= 0; --i) {
-         p = p*x + c[i];
-      }
-      return p;
-   }
-
-   template <int Nstart, typename T, int N>
-   Complex<T> horner(const Complex<T>& z, const T (&coeffs)[N]) noexcept
-   {
-      static_assert(0 <= Nstart && Nstart < N && N >= 2, "invalid array bounds");
-
-      const T r = z.re + z.re;
-      const T s = z.re * z.re + z.im * z.im;
-      T a = coeffs[N - 1], b = coeffs[N - 2];
-
-      for (int i = N - 3; i >= Nstart; --i) {
-         const T t = a;
-         a = b + r * a;
-         b = coeffs[i] - s * t;
-      }
-
-      return Complex<T>(z.re*a + b, z.im*a);
-   }
-
-} // anonymous namespace
-
 /**
  * @brief Real dilogarithm \f$\operatorname{Li}_2(x)\f$
  * @param x real argument
@@ -312,16 +282,15 @@ long double Li2(long double x) noexcept
 
 /**
  * @brief Complex dilogarithm \f$\operatorname{Li}_2(z)\f$
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_2(z)\f$
  * @note Implementation translated from SPheno to C++
  * @author Werner Porod
  * @note translated to C++ by Alexander Voigt
  */
-std::complex<float> Li2(const std::complex<float>& z_) noexcept
+std::complex<float> Li2(const std::complex<float>& z) noexcept
 {
    const float PI = 3.14159265f;
-   const Complex<float> z = { std::real(z_), std::imag(z_) };
 
    // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
    const float bf[] = {
@@ -331,66 +300,69 @@ std::complex<float> Li2(const std::complex<float>& z_) noexcept
       + 1.0f/211680
    };
 
+   const float rz = std::real(z);
+   const float iz = std::imag(z);
+
    // special cases
-   if (z.im == 0) {
-      if (z.re <= 1) {
-         return { Li2(z.re), z.im };
+   if (iz == 0) {
+      if (rz <= 1) {
+         return { Li2(rz), iz };
       }
-      // z.re > 1
-      return { Li2(z.re), -PI*std::log(z.re) };
+      // rz > 1
+      return { Li2(rz), -PI*std::log(rz) };
    }
 
-   const float nz = norm_sqr(z);
+   const float nz = std::norm(z);
 
    if (nz < std::numeric_limits<float>::epsilon()) {
       return z*(1.0f + 0.25f*z);
    }
 
-   Complex<float> u(0.0f, 0.0f), rest(0.0f, 0.0f);
+   std::complex<float> u(0.0f, 0.0f), rest(0.0f, 0.0f);
    float sgn = 1;
 
    // transformation to |z|<1, Re(z)<=0.5
-   if (z.re <= 0.5f) {
+   if (rz <= 0.5f) {
       if (nz > 1) {
-         const Complex<float> lz = log(-z);
-         u = -log(1.0f - 1.0f / z);
+         const auto lz = std::log(-z);
+         u = -log1p(-1.0f/z);
          rest = -0.5f*lz*lz - PI*PI/6;
          sgn = -1;
       } else { // nz <= 1
-         u = -log(1.0f - z);
+         u = -log1p(-z);
          rest = 0;
          sgn = 1;
       }
-   } else { // z.re > 0.5
-      if (nz <= 2*z.re) {
-         u = -log(z);
-         rest = u*log(1.0f - z) + PI*PI/6;
+   } else { // rz > 0.5
+      if (nz <= 2*rz) {
+         u = -std::log(z);
+         rest = u*log1p(-z) + PI*PI/6;
          sgn = -1;
-      } else { // nz > 2*z.re
-         const Complex<float> lz = log(-z);
-         u = -log(1.0f - 1.0f / z);
+      } else { // nz > 2*rz
+         const auto lz = std::log(-z);
+         u = -log1p(-1.0f/z);
          rest = -0.5f*lz*lz - PI*PI/6;
          sgn = -1;
       }
    }
 
-   const Complex<float> u2(u*u);
+   const auto u2(u*u);
 
    return sgn*(u + u2*(bf[0] + u*horner<1>(u2, bf))) + rest;
 }
 
+
 /**
  * @brief Complex dilogarithm \f$\operatorname{Li}_2(z)\f$
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_2(z)\f$
  * @note Implementation translated from SPheno to C++
  * @author Werner Porod
  * @note translated to C++ by Alexander Voigt
  */
-std::complex<double> Li2(const std::complex<double>& z_) noexcept
+std::complex<double> Li2(const std::complex<double>& z) noexcept
 {
    const double PI = 3.1415926535897932;
-   const Complex<double> z = { std::real(z_), std::imag(z_) };
 
    // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
    // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 9}]
@@ -407,66 +379,68 @@ std::complex<double> Li2(const std::complex<double>& z_) noexcept
       + 4.5189800296199182e-16
    };
 
+   const double rz = std::real(z);
+   const double iz = std::imag(z);
+
    // special cases
-   if (z.im == 0) {
-      if (z.re <= 1) {
-         return { Li2(z.re), z.im };
+   if (iz == 0) {
+      if (rz <= 1) {
+         return { Li2(rz), iz };
       }
-      // z.re > 1
-      return { Li2(z.re), -PI*std::log(z.re) };
+      // rz > 1
+      return { Li2(rz), -PI*std::log(rz) };
    }
 
-   const double nz = norm_sqr(z);
+   const double nz = std::norm(z);
 
    if (nz < std::numeric_limits<double>::epsilon()) {
       return z*(1.0 + 0.25*z);
    }
 
-   Complex<double> u(0.0, 0.0), rest(0.0, 0.0);
+   std::complex<double> u(0.0, 0.0), rest(0.0, 0.0);
    double sgn = 1;
 
    // transformation to |z|<1, Re(z)<=0.5
-   if (z.re <= 0.5) {
+   if (rz <= 0.5) {
       if (nz > 1) {
-         const Complex<double> lz = log(-z);
-         u = -log(1.0 - 1.0 / z);
+         const auto lz = std::log(-z);
+         u = -log1p(-1.0/z);
          rest = -0.5*lz*lz - PI*PI/6;
          sgn = -1;
       } else { // nz <= 1
-         u = -log(1.0 - z);
+         u = -log1p(-z);
          rest = 0;
          sgn = 1;
       }
-   } else { // z.re > 0.5
-      if (nz <= 2*z.re) {
-         u = -log(z);
-         rest = u*log(1.0 - z) + PI*PI/6;
+   } else { // rz > 0.5
+      if (nz <= 2*rz) {
+         u = -std::log(z);
+         rest = u*log1p(-z) + PI*PI/6;
          sgn = -1;
-      } else { // nz > 2*z.re
-         const Complex<double> lz = log(-z);
-         u = -log(1.0 - 1.0 / z);
+      } else { // nz > 2*rz
+         const auto lz = std::log(-z);
+         u = -log1p(-1.0/z);
          rest = -0.5*lz*lz - PI*PI/6;
          sgn = -1;
       }
    }
 
-   const Complex<double> u2(u*u);
+   const auto u2(u*u);
 
    return sgn*(u + u2*(bf[0] + u*horner<1>(u2, bf))) + rest;
 }
 
 /**
  * @brief Complex dilogarithm \f$\operatorname{Li}_2(z)\f$ with long double precision
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_2(z)\f$
  * @note Implementation translated from SPheno to C++
  * @author Werner Porod
  * @note translated to C++ and extended to long double precision by Alexander Voigt
  */
-std::complex<long double> Li2(const std::complex<long double>& z_) noexcept
+std::complex<long double> Li2(const std::complex<long double>& z) noexcept
 {
    const long double PI = 3.14159265358979323846264338327950288L;
-   const Complex<long double> z = { std::real(z_), std::imag(z_) };
 
    // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
    // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 22}]
@@ -498,50 +472,53 @@ std::complex<long double> Li2(const std::complex<long double>& z_) noexcept
 #endif
    };
 
+   const long double rz = std::real(z);
+   const long double iz = std::imag(z);
+
    // special cases
-   if (z.im == 0) {
-      if (z.re <= 1) {
-         return { Li2(z.re), z.im };
+   if (iz == 0) {
+      if (rz <= 1) {
+         return { Li2(rz), iz };
       }
-      // z.re > 1
-      return { Li2(z.re), -PI*std::log(z.re) };
+      // rz > 1
+      return { Li2(rz), -PI*std::log(rz) };
    }
 
-   const long double nz = norm_sqr(z);
+   const long double nz = std::norm(z);
 
    if (nz < std::numeric_limits<long double>::epsilon()) {
       return z*(1.0L + 0.25L*z);
    }
 
-   Complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
+   std::complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
    long double sgn = 1;
 
    // transformation to |z|<1, Re(z)<=0.5
-   if (z.re <= 0.5L) {
+   if (rz <= 0.5L) {
       if (nz > 1) {
-         const Complex<long double> lz = log(-z);
-         u = -log(1.0L - 1.0L/z);
+         const auto lz = std::log(-z);
+         u = -log1p(-1.0L/z);
          rest = -0.5L*lz*lz - PI*PI/6;
          sgn = -1;
       } else { // nz <= 1
-         u = -log(1.0L - z);
+         u = -log1p(-z);
          rest = 0;
          sgn = 1;
       }
-   } else { // z.re > 0.5L
-      if (nz <= 2*z.re) {
-         u = -log(z);
-         rest = u*log(1.0L - z) + PI*PI/6;
+   } else { // rz > 0.5L
+      if (nz <= 2*rz) {
+         u = -std::log(z);
+         rest = u*log1p(-z) + PI*PI/6;
          sgn = -1;
-      } else { // nz > 2*z.re
-         const Complex<long double> lz = log(-z);
-         u = -log(1.0L - 1.0L/z);
+      } else { // nz > 2*rz
+         const auto lz = std::log(-z);
+         u = -log1p(-1.0L/z);
          rest = -0.5L*lz*lz - PI*PI/6;
          sgn = -1;
       }
    }
 
-   const Complex<long double> u2(u*u);
+   const auto u2(u*u);
 
    return sgn*(u + u2*(bf[0] + u*horner<1>(u2, bf))) + rest;
 }
diff --git a/src/cpp/Li3.cpp b/src/cpp/Li3.cpp
index 92ff4ff..d7fe2c5 100644
--- a/src/cpp/Li3.cpp
+++ b/src/cpp/Li3.cpp
@@ -5,32 +5,16 @@
 // ====================================================================
 
 #include "Li3.hpp"
-#include "complex.hpp"
+#include "horner.hpp"
+#include "log.hpp"
 #include <cfloat>
 #include <cmath>
+#include <complex>
 
 namespace polylogarithm {
 
 namespace {
 
-   template <typename T, int N>
-   Complex<T> horner(const Complex<T>& z, const T (&coeffs)[N]) noexcept
-   {
-      static_assert(N >= 2, "more than two coefficients required");
-
-      const T r = z.re + z.re;
-      const T s = z.re * z.re + z.im * z.im;
-      T a = coeffs[N - 1], b = coeffs[N - 2];
-
-      for (int i = N - 3; i >= 0; --i) {
-         const T t = a;
-         a = b + r * a;
-         b = coeffs[i] - s * t;
-      }
-
-      return Complex<T>(z.re*a + b, z.im*a);
-   }
-
    /// Li_3(x) for x in [-1,0]
    double li3_neg(double x) noexcept
    {
@@ -126,11 +110,11 @@ double Li3(double x) noexcept
 
 /**
  * @brief Complex trilogarithm \f$\operatorname{Li}_3(z)\f$
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_3(z)\f$
  * @author Alexander Voigt
  */
-std::complex<double> Li3(const std::complex<double>& z_) noexcept
+std::complex<double> Li3(const std::complex<double>& z) noexcept
 {
    const double PI    = 3.1415926535897932;
    const double zeta2 = 1.6449340668482264;
@@ -147,28 +131,29 @@ std::complex<double> Li3(const std::complex<double>& z_) noexcept
       3.1043578879654623e-14,  5.2617586299125061e-15
    };
 
-   const Complex<double> z = { std::real(z_), std::imag(z_) };
+   const double rz = std::real(z);
+   const double iz = std::imag(z);
 
-   if (z.im == 0) {
-      if (z.re <= 1) {
-         return { Li3(z.re), z.im };
+   if (iz == 0) {
+      if (rz <= 1) {
+         return { Li3(rz), iz };
       } else {
-         const double l = std::log(z.re);
-         return { Li3(z.re), -0.5*PI*l*l };
+         const double l = std::log(rz);
+         return { Li3(rz), -0.5*PI*l*l };
       }
    }
 
-   const double nz  = norm(z);
-   const double pz  = arg(z);
+   const double nz  = std::abs(z);
+   const double pz  = std::arg(z);
    const double lnz = std::log(nz);
 
    if (lnz*lnz + pz*pz < 1) { // |log(z)| < 1
-      const Complex<double> u(lnz, pz); // log(z)
-      const Complex<double> u2 = u*u;
-      const Complex<double> u4 = u2*u2;
-      const Complex<double> u8 = u4*u4;
-      const Complex<double> c0 = zeta3 + u*(zeta2 - u2/12.0);
-      const Complex<double> c1 = 0.25 * (3.0 - 2.0*log(-u));
+      const std::complex<double> u(lnz, pz); // log(z)
+      const auto u2 = u*u;
+      const auto u4 = u2*u2;
+      const auto u8 = u4*u4;
+      const auto c0 = zeta3 + u*(zeta2 - u2/12.0);
+      const auto c1 = 0.25 * (3.0 - 2.0*pos_log(-u));
 
       const double cs[7] = {
          -3.4722222222222222e-03, 1.1574074074074074e-05,
@@ -185,20 +170,20 @@ std::complex<double> Li3(const std::complex<double>& z_) noexcept
          u8*u8*cs[6];
    }
 
-   Complex<double> u(0.0, 0.0), rest(0.0, 0.0);
+   std::complex<double> u(0.0, 0.0), rest(0.0, 0.0);
 
    if (nz <= 1) {
-      u = -log(1.0 - z);
+      u = -log1p(-z);
    } else { // nz > 1
       const double arg = pz > 0.0 ? pz - PI : pz + PI;
-      const Complex<double> lmz(lnz, arg); // log(-z)
-      u = -log(1.0 - 1.0/z);
+      const std::complex<double> lmz(lnz, arg); // log(-z)
+      u = -log1p(-1.0/z);
       rest = -lmz*(lmz*lmz/6.0 + zeta2);
    }
 
-   const Complex<double> u2 = u*u;
-   const Complex<double> u4 = u2*u2;
-   const Complex<double> u8 = u4*u4;
+   const auto u2 = u*u;
+   const auto u4 = u2*u2;
+   const auto u8 = u4*u4;
 
    return
       rest +
@@ -212,11 +197,11 @@ std::complex<double> Li3(const std::complex<double>& z_) noexcept
 
 /**
  * @brief Complex trilogarithm \f$\operatorname{Li}_3(z)\f$ with long double precision
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_3(z)\f$
  * @author Alexander Voigt
  */
-std::complex<long double> Li3(const std::complex<long double>& z_) noexcept
+std::complex<long double> Li3(const std::complex<long double>& z) noexcept
 {
    const long double PI    = 3.14159265358979323846264338327950288L;
    const long double zeta2 = 1.64493406684822643647241516664602519L;
@@ -271,32 +256,33 @@ std::complex<long double> Li3(const std::complex<long double>& z_) noexcept
 #endif
    };
 
-   const Complex<long double> z = { std::real(z_), std::imag(z_) };
+   const long double rz = std::real(z);
+   const long double iz = std::imag(z);
 
-   if (z.im == 0) {
-      if (z.re == 0) {
-         return { z.re, z.im };
+   if (iz == 0) {
+      if (rz == 0) {
+         return { rz, iz };
       }
-      if (z.re == 1) {
-         return { zeta3, z.im };
+      if (rz == 1) {
+         return { zeta3, iz };
       }
-      if (z.re == -1) {
-         return { -0.75L*zeta3, z.im };
+      if (rz == -1) {
+         return { -0.75L*zeta3, iz };
       }
-      if (z.re == 0.5L) {
-         return { 0.537213193608040200940623225594965827L, z.im };
+      if (rz == 0.5L) {
+         return { 0.537213193608040200940623225594965827L, iz };
       }
    }
 
-   const long double nz  = norm(z);
-   const long double pz  = arg(z);
+   const long double nz  = std::abs(z);
+   const long double pz  = std::arg(z);
    const long double lnz = std::log(nz);
 
    if (lnz*lnz + pz*pz < 1) { // |log(z)| < 1
-      const Complex<long double> u(lnz, pz); // log(z)
-      const Complex<long double> u2 = u*u;
-      const Complex<long double> c0 = zeta3 + u*(zeta2 - u2/12.0L);
-      const Complex<long double> c1 = 0.25L * (3.0L - 2.0L*log(-u));
+      const std::complex<long double> u(lnz, pz); // log(z)
+      const auto u2 = u*u;
+      const auto c0 = zeta3 + u*(zeta2 - u2/12.0L);
+      const auto c1 = 0.25L * (3.0L - 2.0L*pos_log(-u));
 
       const long double cs[] = {
         -3.47222222222222222222222222222222222e-03L,
@@ -323,21 +309,21 @@ std::complex<long double> Li3(const std::complex<long double>& z_) noexcept
 #endif
       };
 
-      return c0 + u2*(c1 + u2*horner(u2, cs));
+      return c0 + u2*(c1 + u2*horner<0>(u2, cs));
    }
 
-   Complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
+   std::complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
 
    if (nz <= 1) {
-      u = -log(1.0L - z);
+      u = -log1p(-z);
    } else { // nz > 1
       const long double arg = pz > 0.0 ? pz - PI : pz + PI;
-      const Complex<long double> lmz(lnz, arg); // log(-z)
-      u = -log(1.0L - 1.0L/z);
+      const std::complex<long double> lmz(lnz, arg); // log(-z)
+      u = -log1p(-1.0L/z);
       rest = -lmz*(lmz*lmz/6.0L + zeta2);
    }
 
-   return rest + u*horner(u, bf);
+   return rest + u*horner<0>(u, bf);
 }
 
 } // namespace polylogarithm
diff --git a/src/cpp/Li4.cpp b/src/cpp/Li4.cpp
index 3ccb33e..df16968 100644
--- a/src/cpp/Li4.cpp
+++ b/src/cpp/Li4.cpp
@@ -5,32 +5,16 @@
 // ====================================================================
 
 #include "Li4.hpp"
-#include "complex.hpp"
+#include "horner.hpp"
+#include "log.hpp"
 #include <cfloat>
 #include <cmath>
+#include <complex>
 
 namespace polylogarithm {
 
 namespace {
 
-   template <typename T, int N>
-   Complex<T> horner(const Complex<T>& z, const T (&coeffs)[N]) noexcept
-   {
-      static_assert(N >= 2, "more than two coefficients required");
-
-      const T r = z.re + z.re;
-      const T s = z.re * z.re + z.im * z.im;
-      T a = coeffs[N - 1], b = coeffs[N - 2];
-
-      for (int i = N - 3; i >= 0; --i) {
-         const T t = a;
-         a = b + r * a;
-         b = coeffs[i] - s * t;
-      }
-
-      return Complex<T>(z.re*a + b, z.im*a);
-   }
-
    /// Li_4(x) for x in [-1,0]
    double li4_neg(double x) noexcept
    {
@@ -178,11 +162,11 @@ double Li4(double x) noexcept
 
 /**
  * @brief Complex polylogarithm \f$\operatorname{Li}_4(z)\f$
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_4(z)\f$
  * @author Alexander Voigt
  */
-std::complex<double> Li4(const std::complex<double>& z_) noexcept
+std::complex<double> Li4(const std::complex<double>& z) noexcept
 {
    const double PI    = 3.1415926535897932;
    const double PI2   = PI*PI;
@@ -200,57 +184,57 @@ std::complex<double> Li4(const std::complex<double>& z_) noexcept
       2.3647571168618257e-14, -7.9231351220311617e-15
    };
 
-   const Complex<double> z = { std::real(z_), std::imag(z_) };
+   const double rz = std::real(z);
+   const double iz = std::imag(z);
 
-   if (z.im == 0) {
-      if (z.re <= 1) {
-         return { Li4(z.re), z.im };
+   if (iz == 0) {
+      if (rz <= 1) {
+         return { Li4(rz), iz };
       } else {
-         const double l = std::log(z.re);
-         return { Li4(z.re), -1.0/6*PI*l*l*l };
+         const double l = std::log(rz);
+         return { Li4(rz), -1.0/6*PI*l*l*l };
       }
    }
 
-   const double nz  = norm(z);
-   const double pz  = arg(z);
+   const double nz  = std::abs(z);
+   const double pz  = std::arg(z);
    const double lnz = std::log(nz);
 
    if (lnz*lnz + pz*pz < 1) { // |log(z)| < 1
-      const Complex<double> u(lnz, pz); // log(z)
-      const Complex<double> u2 = u*u;
+      const std::complex<double> u(lnz, pz); // log(z)
       const double c1 = 1.2020569031595943; // zeta(3)
       const double c2 = 0.82246703342411322;
-      const Complex<double> c3 = (11.0/6.0 - log(-u))/6.0;
+      const auto c3 = (11.0/6.0 - pos_log(-u))/6.0;
       const double c4 = -1.0/48.0;
-
       const double cs[7] = {
          -6.9444444444444444e-04, 1.6534391534391534e-06,
          -1.0935444136502338e-08, 1.0438378493934049e-10,
          -1.2165942300622435e-12, 1.6130006528350101e-14,
          -2.3428810452879340e-16
       };
+      const auto u2 = u*u;
 
       return zeta4 + u2*(c2 + u2*c4) +
-         u*(c1 + u2*(c3 + u2*horner(u2, cs)));
+         u*(c1 + u2*(c3 + u2*horner<0>(u2, cs)));
    }
 
-   Complex<double> u(0.0, 0.0), rest(0.0, 0.0);
+   std::complex<double> u(0.0, 0.0), rest(0.0, 0.0);
    double sgn = 1;
 
    if (nz <= 1) {
-      u = -log(1.0 - z);
+      u = -log1p(-z);
    } else { // nz > 1
       const double arg = pz > 0.0 ? pz - PI : pz + PI;
-      const Complex<double> lmz(lnz, arg); // log(-z)
-      const Complex<double> lmz2 = lmz*lmz;
-      u = -log(1.0 - 1.0/z);
+      const std::complex<double> lmz(lnz, arg); // log(-z)
+      const auto lmz2 = lmz*lmz;
+      u = -log1p(-1.0/z);
       rest = 1.0/360.0*(-7*PI4 + lmz2*(-30.0*PI2 - 15.0*lmz2));
       sgn = -1;
    }
 
-   const Complex<double> u2 = u*u;
-   const Complex<double> u4 = u2*u2;
-   const Complex<double> u8 = u4*u4;
+   const auto u2 = u*u;
+   const auto u4 = u2*u2;
+   const auto u8 = u4*u4;
 
    return
       rest + sgn * (
@@ -265,11 +249,11 @@ std::complex<double> Li4(const std::complex<double>& z_) noexcept
 
 /**
  * @brief Complex polylogarithm \f$\operatorname{Li}_4(z)\f$ with long double precision
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_4(z)\f$
  * @author Alexander Voigt
  */
-std::complex<long double> Li4(const std::complex<long double>& z_) noexcept
+std::complex<long double> Li4(const std::complex<long double>& z) noexcept
 {
    const long double PI    = 3.14159265358979323846264338327950288L;
    const long double PI2   = PI*PI;
@@ -325,32 +309,31 @@ std::complex<long double> Li4(const std::complex<long double>& z_) noexcept
 #endif
    };
 
-   const Complex<long double> z = { std::real(z_), std::imag(z_) };
+   const long double rz = std::real(z);
+   const long double iz = std::imag(z);
 
-   if (z.im == 0) {
-      if (z.re == 0) {
-         return { z.re, z.im };
+   if (iz == 0) {
+      if (rz == 0) {
+         return { rz, iz };
       }
-      if (z.re == 1) {
-         return { zeta4, z.im };
+      if (rz == 1) {
+         return { zeta4, iz };
       }
-      if (z.re == -1) {
-         return { -7.0L*PI4/720.0L, z.im };
+      if (rz == -1) {
+         return { -7.0L*PI4/720.0L, iz };
       }
    }
 
-   const long double nz  = norm(z);
-   const long double pz  = arg(z);
+   const long double nz  = std::abs(z);
+   const long double pz  = std::arg(z);
    const long double lnz = std::log(nz);
 
    if (lnz*lnz + pz*pz < 1) { // |log(z)| < 1
-      const Complex<long double> u(lnz, pz); // log(z)
-      const Complex<long double> u2 = u*u;
+      const std::complex<long double> u(lnz, pz); // log(z)
       const long double c1 = 1.20205690315959428539973816151144999L; // zeta(3)
       const long double c2 = 0.822467033424113218236207583323012595L;
-      const Complex<long double> c3 = (11.0L/6.0L - log(-u))/6.0L;
+      const auto c3 = (11.0L/6.0L - pos_log(-u))/6.0L;
       const long double c4 = -1.0L/48.0L;
-
       const long double cs[] = {
         -6.94444444444444444444444444444444444e-04L,
          1.65343915343915343915343915343915344e-06L,
@@ -374,26 +357,27 @@ std::complex<long double> Li4(const std::complex<long double>& z_) noexcept
         -3.84151865355610606630482668147264051e-37L
 #endif
       };
+      const auto u2 = u*u;
 
       return zeta4 + u2*(c2 + u2*c4) +
-         u*(c1 + u2*(c3 + u2*horner(u2, cs)));
+         u*(c1 + u2*(c3 + u2*horner<0>(u2, cs)));
    }
 
-   Complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
+   std::complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
    long double sgn = 1;
 
    if (nz <= 1) {
-      u = -log(1.0L - z);
+      u = -log1p(-z);
    } else { // nz > 1
       const long double arg = pz > 0.0 ? pz - PI : pz + PI;
-      const Complex<long double> lmz(lnz, arg); // log(-z)
-      const Complex<long double> lmz2 = lmz*lmz;
-      u = -log(1.0L - 1.0L/z);
+      const std::complex<long double> lmz(lnz, arg); // log(-z)
+      const auto lmz2 = lmz*lmz;
+      u = -log1p(-1.0L/z);
       rest = 1.0L/360.0L*(-7*PI4 + lmz2*(-30.0L*PI2 - 15.0L*lmz2));
       sgn = -1;
    }
 
-   return rest + sgn*u*horner(u, bf);
+   return rest + sgn*u*horner<0>(u, bf);
 }
 
 } // namespace polylogarithm
diff --git a/src/cpp/Li5.cpp b/src/cpp/Li5.cpp
index 975eb9c..9706eb5 100644
--- a/src/cpp/Li5.cpp
+++ b/src/cpp/Li5.cpp
@@ -5,41 +5,21 @@
 // ====================================================================
 
 #include "Li5.hpp"
-#include "complex.hpp"
+#include "horner.hpp"
+#include "log.hpp"
 #include <cfloat>
 #include <cmath>
+#include <complex>
 
 namespace polylogarithm {
 
-namespace {
-
-   template <typename T, int N>
-   Complex<T> horner(const Complex<T>& z, const T (&coeffs)[N]) noexcept
-   {
-      static_assert(N >= 2, "more than two coefficients required");
-
-      const T r = z.re + z.re;
-      const T s = z.re * z.re + z.im * z.im;
-      T a = coeffs[N - 1], b = coeffs[N - 2];
-
-      for (int i = N - 3; i >= 0; --i) {
-         const T t = a;
-         a = b + r * a;
-         b = coeffs[i] - s * t;
-      }
-
-      return Complex<T>(z.re*a + b, z.im*a);
-   }
-
-} // anonymous namespace
-
 /**
  * @brief Complex polylogarithm \f$\operatorname{Li}_5(z)\f$
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_5(z)\f$
  * @author Alexander Voigt
  */
-std::complex<double> Li5(const std::complex<double>& z_) noexcept
+std::complex<double> Li5(const std::complex<double>& z) noexcept
 {
    const double PI    = 3.1415926535897932;
    const double PI2   = PI*PI;
@@ -58,61 +38,61 @@ std::complex<double> Li5(const std::complex<double>& z_) noexcept
       1.0903545401333394e-15
    };
 
-   const Complex<double> z = { std::real(z_), std::imag(z_) };
+   const double rz = std::real(z);
+   const double iz = std::imag(z);
 
-   if (z.im == 0) {
-      if (z.re == 0) {
-         return { z.re, z.im };
+   if (iz == 0) {
+      if (rz == 0) {
+         return { rz, iz };
       }
-      if (z.re == 1) {
-         return { zeta5, z.im };
+      if (rz == 1) {
+         return { zeta5, iz };
       }
-      if (z.re == -1) {
-         return { -15.0*zeta5/16.0, z.im };
+      if (rz == -1) {
+         return { -15.0*zeta5/16.0, iz };
       }
    }
 
-   const double nz  = norm(z);
-   const double pz  = arg(z);
+   const double nz  = std::abs(z);
+   const double pz  = std::arg(z);
    const double lnz = std::log(nz);
 
    if (lnz*lnz + pz*pz < 1) { // |log(z)| < 1
-      const Complex<double> u(lnz, pz); // log(z)
-      const Complex<double> u2 = u*u;
+      const std::complex<double> u(lnz, pz); // log(z)
       const double c0 = zeta5;
       const double c1 = 1.0823232337111382; // zeta(4)
       const double c2 = 0.60102845157979714; // zeta(3)/2
       const double c3 = 0.27415567780803774;
-      const Complex<double> c4 = (25.0/12.0 - log(-u))/24.0;
+      const auto c4 = (25.0/12.0 - pos_log(-u))/24.0;
       const double c5 = -1.0/240.0;
-
       const double cs[6] = {
          -1.1574074074074074e-04, 2.0667989417989418e-07,
          -1.0935444136502338e-09, 8.6986487449450412e-12,
          -8.6899587861588824e-14, 1.0081254080218813e-15
       };
+      const auto u2 = u*u;
 
       return c0 + u * c1 +
          u2 * (c2 + u * c3 +
          u2 * (c4 + u * c5 +
-         u2 * horner(u2, cs)));
+         u2 * horner<0>(u2, cs)));
    }
 
-   Complex<double> u(0.0, 0.0), rest(0.0, 0.0);
+   std::complex<double> u(0.0, 0.0), rest(0.0, 0.0);
 
    if (nz <= 1) {
-      u = -log(1.0 - z);
+      u = -log1p(-z);
    } else { // nz > 1
       const double arg = pz > 0.0 ? pz - PI : pz + PI;
-      const Complex<double> lmz(lnz, arg); // log(-z)
-      const Complex<double> lmz2 = lmz*lmz;
-      u = -log(1.0 - 1.0/z);
+      const std::complex<double> lmz(lnz, arg); // log(-z)
+      const auto lmz2 = lmz*lmz;
+      u = -log1p(-1.0/z);
       rest = -1.0/360.0*lmz*(7*PI4 + lmz2*(10.0*PI2 + 3.0*lmz2));
    }
 
-   const Complex<double> u2 = u*u;
-   const Complex<double> u4 = u2*u2;
-   const Complex<double> u8 = u4*u4;
+   const auto u2 = u*u;
+   const auto u4 = u2*u2;
+   const auto u8 = u4*u4;
 
    return
       rest +
@@ -126,11 +106,11 @@ std::complex<double> Li5(const std::complex<double>& z_) noexcept
 
 /**
  * @brief Complex polylogarithm \f$\operatorname{Li}_5(z)\f$ with long double precision
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_5(z)\f$
  * @author Alexander Voigt
  */
-std::complex<long double> Li5(const std::complex<long double>& z_) noexcept
+std::complex<long double> Li5(const std::complex<long double>& z) noexcept
 {
    const long double PI    = 3.14159265358979323846264338327950288L;
    const long double PI2   = PI*PI;
@@ -186,34 +166,33 @@ std::complex<long double> Li5(const std::complex<long double>& z_) noexcept
 #endif
    };
 
-   const Complex<long double> z = { std::real(z_), std::imag(z_) };
+   const long double rz = std::real(z);
+   const long double iz = std::imag(z);
 
-   if (z.im == 0) {
-      if (z.re == 0) {
-         return { z.re, z.im };
+   if (iz == 0) {
+      if (rz == 0) {
+         return { rz, iz };
       }
-      if (z.re == 1) {
-         return { zeta5, z.im };
+      if (rz == 1) {
+         return { zeta5, iz };
       }
-      if (z.re == -1) {
-         return { -15.0L*zeta5/16.0L, z.im };
+      if (rz == -1) {
+         return { -15.0L*zeta5/16.0L, iz };
       }
    }
 
-   const long double nz  = norm(z);
-   const long double pz  = arg(z);
+   const long double nz  = std::abs(z);
+   const long double pz  = std::arg(z);
    const long double lnz = std::log(nz);
 
    if (lnz*lnz + pz*pz < 1) { // |log(z)| < 1
-      const Complex<long double> u(lnz, pz); // log(z)
-      const Complex<long double> u2 = u*u;
+      const std::complex<long double> u(lnz, pz); // log(z)
       const long double c0 = zeta5;
       const long double c1 = 1.08232323371113819151600369654116790L; // zeta(4)
       const long double c2 = 0.601028451579797142699869080755724995L; // zeta(3)/2
       const long double c3 = 0.274155677808037739412069194441004198L;
-      const Complex<long double> c4 = (25.0L/12.0L - log(-u))/24.0L;
+      const auto c4 = (25.0L/12.0L - pos_log(-u))/24.0L;
       const long double c5 = -1.0L/240.0L;
-
       const long double cs[] = {
         -1.15740740740740740740740740740740741e-04L,
          2.06679894179894179894179894179894180e-07L,
@@ -236,26 +215,27 @@ std::complex<long double> Li5(const std::complex<long double>& z_) noexcept
          4.66812326074215685597260023169508118e-37L
 #endif
       };
+      const auto u2 = u*u;
 
       return c0 + u * c1 +
          u2 * (c2 + u * c3 +
          u2 * (c4 + u * c5 +
-         u2 * horner(u2, cs)));
+         u2 * horner<0>(u2, cs)));
    }
 
-   Complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
+   std::complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
 
    if (nz <= 1) {
-      u = -log(1.0L - z);
+      u = -log1p(-z);
    } else { // nz > 1
       const long double arg = pz > 0.0 ? pz - PI : pz + PI;
-      const Complex<long double> lmz(lnz, arg); // log(-z)
-      const Complex<long double> lmz2 = lmz*lmz;
-      u = -log(1.0L - 1.0L/z);
+      const std::complex<long double> lmz(lnz, arg); // log(-z)
+      const auto lmz2 = lmz*lmz;
+      u = -log1p(-1.0L/z);
       rest = -1.0L/360.0L*lmz*(7*PI4 + lmz2*(10.0L*PI2 + 3.0L*lmz2));
    }
 
-   return rest + u*horner(u, bf);
+   return rest + u*horner<0>(u, bf);
 }
 
 } // namespace polylogarithm
diff --git a/src/cpp/Li6.cpp b/src/cpp/Li6.cpp
index 7098805..aa53d6c 100644
--- a/src/cpp/Li6.cpp
+++ b/src/cpp/Li6.cpp
@@ -5,41 +5,21 @@
 // ====================================================================
 
 #include "Li6.hpp"
-#include "complex.hpp"
+#include "horner.hpp"
+#include "log.hpp"
 #include <cfloat>
 #include <cmath>
+#include <complex>
 
 namespace polylogarithm {
 
-namespace {
-
-   template <typename T, int N>
-   Complex<T> horner(const Complex<T>& z, const T (&coeffs)[N]) noexcept
-   {
-      static_assert(N >= 2, "more than two coefficients required");
-
-      const T r = z.re + z.re;
-      const T s = z.re * z.re + z.im * z.im;
-      T a = coeffs[N - 1], b = coeffs[N - 2];
-
-      for (int i = N - 3; i >= 0; --i) {
-         const T t = a;
-         a = b + r * a;
-         b = coeffs[i] - s * t;
-      }
-
-      return Complex<T>(z.re*a + b, z.im*a);
-   }
-
-} // anonymous namespace
-
 /**
  * @brief Complex polylogarithm \f$\operatorname{Li}_6(z)\f$
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_6(z)\f$
  * @author Alexander Voigt
  */
-std::complex<double> Li6(const std::complex<double>& z_) noexcept
+std::complex<double> Li6(const std::complex<double>& z) noexcept
 {
    const double PI    = 3.1415926535897932;
    const double PI2   = PI*PI;
@@ -58,65 +38,65 @@ std::complex<double> Li6(const std::complex<double>& z_) noexcept
      -6.8408811719011698e-15,  4.8691178462005581e-15
    };
 
-   const Complex<double> z = { std::real(z_), std::imag(z_) };
+   const double rz = std::real(z);
+   const double iz = std::imag(z);
 
-   if (z.im == 0) {
-      if (z.re == 0) {
-         return { z.re, z.im };
+   if (iz == 0) {
+      if (rz == 0) {
+         return { rz, iz };
       }
-      if (z.re == 1) {
-         return { zeta6, z.im };
+      if (rz == 1) {
+         return { zeta6, iz };
       }
-      if (z.re == -1) {
-         return { -31.0*zeta6/32.0, z.im };
+      if (rz == -1) {
+         return { -31.0*zeta6/32.0, iz };
       }
    }
 
-   const double nz  = norm(z);
-   const double pz  = arg(z);
+   const double nz  = std::abs(z);
+   const double pz  = std::arg(z);
    const double lnz = std::log(nz);
 
    if (lnz*lnz + pz*pz < 1) { // |log(z)| < 1
-      const Complex<double> u(lnz, pz); // log(z)
-      const Complex<double> u2 = u*u;
+      const std::complex<double> u(lnz, pz); // log(z)
       const double c0 = zeta6;
       const double c1 = 1.0369277551433699; // zeta(5)
       const double c2 = 0.54116161685556910;
       const double c3 = 0.20034281719326571;
       const double c4 = 0.068538919452009435;
-      const Complex<double> c5 = (137.0/60.0 - log(-u))/120.0;
+      const auto c5 = (137.0/60.0 - pos_log(-u))/120.0;
       const double c6 = -1.0/1440.0;
-
       const double cs[5] = {
          -1.6534391534391534e-05, 2.2964432686654909e-08,
          -9.9413128513657614e-11, 6.6912682653423394e-13,
          -5.7933058574392549e-15
       };
+      const auto u2 = u*u;
 
       return c0 + u * c1 +
          u2 * (c2 + u * c3 +
          u2 * (c4 + u * c5 +
          u2 * (c6 +
-         u * horner(u2, cs))));
+         u * horner<0>(u2, cs))));
    }
 
-   Complex<double> u(0.0, 0.0), rest(0.0, 0.0);
+   std::complex<double> u(0.0, 0.0), rest(0.0, 0.0);
    double sgn = 1;
 
    if (nz <= 1) {
-      u = -log(1.0 - z);
+      u = -log1p(-z);
    } else { // nz > 1
       const double arg = pz > 0.0 ? pz - PI : pz + PI;
-      const Complex<double> lmz(lnz, arg); // log(-z)
-      const Complex<double> lmz2 = lmz*lmz;
-      u = -log(1.0 - 1.0/z);
+      const std::complex<double> lmz(lnz, arg); // log(-z)
+      const auto lmz2 = lmz*lmz;
+      u = -log1p(-1.0/z);
       rest = -31.0*PI6/15120.0 + lmz2*(-7.0/720.0*PI4 + lmz2*(-1.0/144.0*PI2 - 1.0/720.0*lmz2));
       sgn = -1;
    }
 
-   const Complex<double> u2 = u*u;
-   const Complex<double> u4 = u2*u2;
-   const Complex<double> u8 = u4*u4;
+   const auto u2 = u*u;
+   const auto u4 = u2*u2;
+   const auto u8 = u4*u4;
 
    return
       rest + sgn * (
@@ -131,11 +111,11 @@ std::complex<double> Li6(const std::complex<double>& z_) noexcept
 
 /**
  * @brief Complex polylogarithm \f$\operatorname{Li}_6(z)\f$ with long double precision
- * @param z_ complex argument
+ * @param z complex argument
  * @return \f$\operatorname{Li}_6(z)\f$
  * @author Alexander Voigt
  */
-std::complex<long double> Li6(const std::complex<long double>& z_) noexcept
+std::complex<long double> Li6(const std::complex<long double>& z) noexcept
 {
    const long double PI    = 3.14159265358979323846264338327950288L;
    const long double PI2   = PI*PI;
@@ -192,35 +172,34 @@ std::complex<long double> Li6(const std::complex<long double>& z_) noexcept
 #endif
    };
 
-   const Complex<long double> z = { std::real(z_), std::imag(z_) };
+   const long double rz = std::real(z);
+   const long double iz = std::imag(z);
 
-   if (z.im == 0) {
-      if (z.re == 0) {
-         return { z.re, z.im };
+   if (iz == 0) {
+      if (rz == 0) {
+         return { rz, iz };
       }
-      if (z.re == 1) {
-         return { zeta6, z.im };
+      if (rz == 1) {
+         return { zeta6, iz };
       }
-      if (z.re == -1) {
-         return { -31.0L*zeta6/32.0L, z.im };
+      if (rz == -1) {
+         return { -31.0L*zeta6/32.0L, iz };
       }
    }
 
-   const long double nz  = norm(z);
-   const long double pz  = arg(z);
+   const long double nz  = std::abs(z);
+   const long double pz  = std::arg(z);
    const long double lnz = std::log(nz);
 
    if (lnz*lnz + pz*pz < 1) { // |log(z)| < 1
-      const Complex<long double> u(lnz, pz); // log(z)
-      const Complex<long double> u2 = u*u;
+      const std::complex<long double> u(lnz, pz); // log(z)
       const long double c0 = zeta6;
       const long double c1 = 1.03692775514336992633136548645703417L; // zeta(5)
       const long double c2 = 0.541161616855569095758001848270583951L;
       const long double c3 = 0.200342817193265714233289693585241665L;
       const long double c4 = 0.0685389194520094348530172986102510496L;
-      const Complex<long double> c5 = (137.0L/60.0L - log(-u))/120.0L;
+      const auto c5 = (137.0L/60.0L - pos_log(-u))/120.0L;
       const long double c6 = -1.0L/1440.0L;
-
       const long double cs[] = {
         -1.65343915343915343915343915343915344e-05L,
          2.29644326866549088771310993533215755e-08L,
@@ -242,31 +221,32 @@ std::complex<long double> Li6(const std::complex<long double>& z_) noexcept
         -6.19462586636097236736900573587284992e-37L
 #endif
       };
+      const auto u2 = u*u;
 
       return c0 + u * c1 +
          u2 * (c2 + u * c3 +
          u2 * (c4 + u * c5 +
          u2 * (c6 +
-         u * horner(u2, cs))));
+         u * horner<0>(u2, cs))));
    }
 
-   Complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
+   std::complex<long double> u(0.0L, 0.0L), rest(0.0L, 0.0L);
    long double sgn = 1;
 
    if (nz <= 1) {
-      u = -log(1.0L - z);
+      u = -log1p(-z);
    } else { // nz > 1
       const long double arg = pz > 0.0 ? pz - PI : pz + PI;
-      const Complex<long double> lmz(lnz, arg); // log(-z)
-      const Complex<long double> lmz2 = lmz*lmz;
-      u = -log(1.0L - 1.0L/z);
+      const std::complex<long double> lmz(lnz, arg); // log(-z)
+      const auto lmz2 = lmz*lmz;
+      u = -log1p(-1.0L/z);
       rest = -31.0L*PI6/15120.0L
              + lmz2*(-7.0L/720.0L*PI4 +
                      lmz2*(-1.0L/144.0L*PI2 - 1.0L/720.0L*lmz2));
       sgn = -1;
    }
 
-   return rest + sgn*u*horner(u, bf);
+   return rest + sgn*u*horner<0>(u, bf);
 }
 
 } // namespace polylogarithm
diff --git a/src/cpp/complex.hpp b/src/cpp/complex.hpp
deleted file mode 100644
index 567658f..0000000
--- a/src/cpp/complex.hpp
+++ /dev/null
@@ -1,129 +0,0 @@
-// ====================================================================
-// This file is part of Polylogarithm.
-//
-// Polylogarithm is licenced under the MIT License.
-// ====================================================================
-
-#pragma once
-
-#include <cmath>
-#include <complex>
-
-namespace polylogarithm {
-
-template <typename T>
-struct Complex {
-   constexpr Complex(T re_ = T{}, T im_ = T{}) : re(re_), im(im_) {}
-   operator std::complex<T>() const noexcept { return std::complex<T>(re, im); }
-   T re{};
-   T im{};
-};
-
-template <typename T>
-constexpr T arg(const Complex<T>& z) noexcept
-{
-   return std::atan2(z.im, z.re);
-}
-
-template <typename T>
-constexpr Complex<T> conj(const Complex<T>& z) noexcept
-{
-   return { z.re, -z.im };
-}
-
-template <typename T>
-Complex<T> log(const Complex<T>& z) noexcept
-{
-   if (z.im == T(0) && z.re > T(0)) {
-      return { std::log(z.re), T(0) };
-   } else if (z.im == T(0)) {
-      return { std::log(-z.re), 4*std::atan(T(1)) };
-   }
-   return { std::log(norm(z)), arg(z) };
-}
-
-template <typename T>
-constexpr T norm(const Complex<T>& z) noexcept
-{
-   return std::hypot(z.re, z.im);
-}
-
-template <typename T>
-constexpr T norm_sqr(const Complex<T>& z) noexcept
-{
-   return z.re*z.re + z.im*z.im;
-}
-
-template <typename T>
-constexpr Complex<T> operator+(const Complex<T>& a, const Complex<T>& b) noexcept
-{
-   return { a.re + b.re, a.im + b.im };
-}
-
-template <typename T>
-constexpr Complex<T> operator+(const Complex<T>& z, T x) noexcept
-{
-   return { z.re + x, z.im };
-}
-
-template <typename T>
-constexpr Complex<T> operator+(T x, const Complex<T>& z) noexcept
-{
-   return { x + z.re, z.im };
-}
-
-template <typename T>
-constexpr Complex<T> operator-(const Complex<T>& a, const Complex<T>& b) noexcept
-{
-   return { a.re - b.re, a.im - b.im };
-}
-
-template <typename T>
-constexpr Complex<T> operator-(T x, const Complex<T>& z) noexcept
-{
-   return { x - z.re, -z.im };
-}
-
-template <typename T>
-constexpr Complex<T> operator-(const Complex<T>& z, T x) noexcept
-{
-   return { z.re - x, z.im };
-}
-
-template <typename T>
-constexpr Complex<T> operator-(const Complex<T>& z) noexcept
-{
-   return { -z.re, -z.im };
-}
-
-template <typename T>
-constexpr Complex<T> operator*(const Complex<T>& a, const Complex<T>& b) noexcept
-{
-   return { a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re };
-}
-
-template <typename T>
-constexpr Complex<T> operator*(T x, const Complex<T>& z) noexcept
-{
-   return { x*z.re, x*z.im };
-}
-
-template <typename T>
-constexpr Complex<T> operator*(const Complex<T>& z, T x) noexcept
-{
-   return x*z;
-}
-
-template <typename T>
-constexpr Complex<T> operator/(T x, const Complex<T>& z) noexcept
-{
-   return x*conj(z)/norm_sqr(z);
-}
-
-template <typename T>
-constexpr Complex<T> operator/(const Complex<T>& z, T x) noexcept
-{
-   return { z.re/x, z.im/x };
-}
-
-} // namespace polylogarithm
diff --git a/src/cpp/horner.hpp b/src/cpp/horner.hpp
new file mode 100644
index 0000000..e7be272
--- /dev/null
+++ b/src/cpp/horner.hpp
@@ -0,0 +1,44 @@
+// ====================================================================
+// This file is part of Polylogarithm.
+//
+// Polylogarithm is licenced under the MIT License.
+// ====================================================================
+
+#pragma once
+
+#include <complex>
+
+namespace polylogarithm {
+
+template <typename T, int N>
+T horner(T x, const T (&c)[N]) noexcept
+{
+   T p = c[N - 1];
+   for (int i = N - 2; i >= 0; --i) {
+      p = p*x + c[i];
+   }
+   return p;
+}
+
+
+template <int Nstart, typename T, int N>
+std::complex<T> horner(const std::complex<T>& z, const T (&coeffs)[N]) noexcept
+{
+   static_assert(0 <= Nstart && Nstart < N && N >= 2, "invalid array bounds");
+
+   const T rz = std::real(z);
+   const T iz = std::imag(z);
+   const T r = rz + rz;
+   const T s = std::norm(z);
+   T a = coeffs[N - 1], b = coeffs[N - 2];
+
+   for (int i = N - 3; i >= Nstart; --i) {
+      const T t = a;
+      a = b + r*a;
+      b = coeffs[i] - s*t;
+   }
+
+   return { rz*a + b, iz*a };
+}
+
+} // namespace polylogarithm
diff --git a/src/cpp/log.hpp b/src/cpp/log.hpp
new file mode 100644
index 0000000..d3e7167
--- /dev/null
+++ b/src/cpp/log.hpp
@@ -0,0 +1,43 @@
+// ====================================================================
+// This file is part of Polylogarithm.
+//
+// Polylogarithm is licenced under the MIT License.
+// ====================================================================
+
+#pragma once
+
+#include <cmath>
+#include <complex>
+
+namespace polylogarithm {
+
+
+template <typename T>
+std::complex<T> log1p(const std::complex<T>& z) noexcept
+{
+   const std::complex<T> u = T(1) + z;
+
+   if (std::real(u) == T(1) && std::imag(u) == T(0)) {
+      return z;
+   } else if (std::real(u) <= T(0)) {
+      return std::log(u);
+   }
+
+   return std::log(u)*(z/(u - T(1)));
+}
+
+
+template <typename T>
+std::complex<T> pos_log(const std::complex<T>& z) noexcept
+{
+   if (std::imag(z) == T(0) && std::real(z) > T(0)) {
+      return { std::log(std::real(z)), T(0) };
+   } else if (std::imag(z) == T(0)) {
+      return { std::log(-std::real(z)), 4*std::atan(T(1)) };
+   }
+
+   return std::log(z);
+}
+
+
+} // namespace polylogarithm
diff --git a/src/fortran/CMakeLists.txt b/src/fortran/CMakeLists.txt
index 463a967..fe71f42 100644
--- a/src/fortran/CMakeLists.txt
+++ b/src/fortran/CMakeLists.txt
@@ -1,6 +1,5 @@
 if(CMAKE_Fortran_COMPILER)
   add_library(polylogarithm_fortran
-    fast_clog.f90
     Cl2.f90
     Cl3.f90
     Cl4.f90
@@ -11,6 +10,7 @@ if(CMAKE_Fortran_COMPILER)
     Li4.f90
     Li5.f90
     Li6.f90
+    log.f90
     )
   target_include_directories(polylogarithm_fortran PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
   set(FORTRAN_HEADERS
diff --git a/src/fortran/Li2.f90 b/src/fortran/Li2.f90
index b83a5fc..8eaa847 100644
--- a/src/fortran/Li2.f90
+++ b/src/fortran/Li2.f90
@@ -95,7 +95,7 @@ end function dli2
 
 double complex function cdli2(z)
   implicit none
-  double complex :: z, rest, u, u2, u4, sum, fast_cdlog
+  double complex :: z, rest, u, u2, u4, sum, cdlog1p
   double precision :: rz, iz, nz, sgn, dli2
   double precision, parameter :: PI = 3.14159265358979324D0
   double precision, parameter :: bf(10) = (/ &
@@ -130,22 +130,22 @@ double complex function cdli2(z)
   ! transformation to |z| < 1, Re(z) <= 0.5
   if (rz .le. 0.5D0) then
      if (nz .gt. 1) then
-        u = -fast_cdlog(1 - 1/z)
-        rest = -0.5D0*fast_cdlog(-z)**2 - PI**2/6
+        u = -cdlog1p(-1/z)
+        rest = -0.5D0*log(-z)**2 - PI**2/6
         sgn = -1
      else ! nz <= 1
-        u = -fast_cdlog(1 - z)
+        u = -cdlog1p(-z)
         rest = 0
         sgn = 1
      endif
   else ! rz > 0.5D0
      if (nz .le. 2*rz) then
-        u = -fast_cdlog(z)
-        rest = u*fast_cdlog(1 - z) + PI**2/6
+        u = -log(z)
+        rest = u*cdlog1p(-z) + PI**2/6
         sgn = -1
      else ! nz > 2*rz
-        u = -fast_cdlog(1 - 1/z)
-        rest = -0.5D0*fast_cdlog(-z)**2 - PI**2/6
+        u = -cdlog1p(-1/z)
+        rest = -0.5D0*log(-z)**2 - PI**2/6
         sgn = -1
      endif
   endif
diff --git a/src/fortran/Li3.f90 b/src/fortran/Li3.f90
index a616e28..6ce02a8 100644
--- a/src/fortran/Li3.f90
+++ b/src/fortran/Li3.f90
@@ -109,7 +109,7 @@ end function dli3
 !*********************************************************************
 double complex function cdli3(z)
   implicit none
-  double complex :: z, u, u2, u4, u8, c0, c1, lmz, rest, fast_pos_cdlog
+  double complex :: z, u, u2, u4, u8, c0, c1, lmz, rest, pos_cdlog
   double precision :: rz, iz, nz, pz, lnz, arg, dli3
   double precision, parameter :: PI    = 3.1415926535897932D0
   double precision, parameter :: zeta2 = 1.6449340668482264D0
@@ -154,7 +154,7 @@ double complex function cdli3(z)
       u4 = u2**2
       u8 = u4**2
       c0 = zeta3 + u*(zeta2 - u2/12)
-      c1 = 0.25D0*(3 - 2*fast_pos_cdlog(-u))
+      c1 = 0.25D0*(3 - 2*pos_cdlog(-u))
       cdli3 =                                            &
          c0 +                                            &
          c1*u2 +                                         &
@@ -165,7 +165,7 @@ double complex function cdli3(z)
    endif
 
    if (nz .le. 1) then
-      u = -fast_pos_cdlog(1 - z)
+      u = -pos_cdlog(1 - z)
       rest = 0
    else ! nz > 1
       if (pz .gt. 0) then
@@ -174,7 +174,7 @@ double complex function cdli3(z)
          arg = pz + PI
       endif
       lmz = dcmplx(lnz, arg) ! log(-z)
-      u = -fast_pos_cdlog(1 - 1/z)
+      u = -pos_cdlog(1 - 1/z)
       rest = -lmz*(lmz**2/6 + zeta2)
    endif
 
diff --git a/src/fortran/Li4.f90 b/src/fortran/Li4.f90
index d2fe6d9..3e1f257 100644
--- a/src/fortran/Li4.f90
+++ b/src/fortran/Li4.f90
@@ -164,7 +164,7 @@ end function dli4
 !*********************************************************************
 double complex function cdli4(z)
   implicit none
-  double complex :: z, u, u2, u4, u8, c3, lmz, r, fast_pos_cdlog
+  double complex :: z, u, u2, u4, u8, c3, lmz, r, pos_cdlog
   double precision :: rz, iz, nz, pz, lnz, arg, sgn, dli4
   double precision, parameter :: PI    = 3.1415926535897932D0
   double precision, parameter :: PI2   = 9.8696044010893586D0
@@ -212,7 +212,7 @@ double complex function cdli4(z)
       u2 = u**2
       u4 = u2**2
       u8 = u4**2
-      c3 = (11D0/6 - fast_pos_cdlog(-u))/6
+      c3 = (11D0/6 - pos_cdlog(-u))/6
       cdli4 = zeta4 + u2 * (c2 + u2 * c4) +                   &
           u * (                                               &
               c1 +                                            &
@@ -225,7 +225,7 @@ double complex function cdli4(z)
    endif
 
    if (nz .le. 1) then
-      u = -fast_pos_cdlog(1 - z)
+      u = -pos_cdlog(1 - z)
       r = 0
       sgn = 1
    else ! nz > 1
@@ -235,7 +235,7 @@ double complex function cdli4(z)
          arg = pz + PI
       endif
       lmz = dcmplx(lnz, arg) ! log(-z)
-      u = -fast_pos_cdlog(1 - 1/z)
+      u = -pos_cdlog(1 - 1/z)
       r = (-7*PI4 + lmz**2*(-30*PI2 - 15*lmz**2))/360
       sgn = -1
    endif
diff --git a/src/fortran/Li5.f90 b/src/fortran/Li5.f90
index 754b8ba..46900f1 100644
--- a/src/fortran/Li5.f90
+++ b/src/fortran/Li5.f90
@@ -13,7 +13,7 @@
 !*********************************************************************
 double complex function cdli5(z)
   implicit none
-  double complex :: z, u, u2, u4, u8, c4, lmz, rest, fast_pos_cdlog
+  double complex :: z, u, u2, u4, u8, c4, lmz, rest, pos_cdlog
   double precision :: rz, iz, nz, pz, lnz, arg
   double precision, parameter :: PI    = 3.1415926535897932D0
   double precision, parameter :: PI2   = 9.8696044010893586D0
@@ -64,7 +64,7 @@ double complex function cdli5(z)
    if (lnz**2 + pz**2 .lt. 1) then ! |log(z)| < 1
       u = dcmplx(lnz, pz) ! log(z)
       u2 = u**2
-      c4 = (25D0/12 - fast_pos_cdlog(-u))/24
+      c4 = (25D0/12 - pos_cdlog(-u))/24
       cdli5 =                &
          zeta5 + u * c1 +    &
          u2 * (c2 + u * c3 + &
@@ -79,7 +79,7 @@ double complex function cdli5(z)
    endif
 
    if (nz .le. 1) then
-      u = -fast_pos_cdlog(1 - z)
+      u = -pos_cdlog(1 - z)
       rest = 0
    else ! nz > 1
       if (pz .gt. 0) then
@@ -88,7 +88,7 @@ double complex function cdli5(z)
          arg = pz + PI
       endif
       lmz = dcmplx(lnz, arg) ! log(-z)
-      u = -fast_pos_cdlog(1 - 1/z)
+      u = -pos_cdlog(1 - 1/z)
       rest = -lmz*(7*PI4 + lmz**2*(10*PI2 + 3*lmz**2))/360
    endif
 
diff --git a/src/fortran/Li6.f90 b/src/fortran/Li6.f90
index d065530..4c05c69 100644
--- a/src/fortran/Li6.f90
+++ b/src/fortran/Li6.f90
@@ -13,7 +13,7 @@
 !*********************************************************************
 double complex function cdli6(z)
   implicit none
-  double complex :: z, u, u2, u4, u8, c5, lmz, r, fast_pos_cdlog
+  double complex :: z, u, u2, u4, u8, c5, lmz, r, pos_cdlog
   double precision :: rz, iz, nz, pz, lnz, arg, sgn
   double precision, parameter :: PI    = 3.1415926535897932D0
   double precision, parameter :: PI2   = 9.8696044010893586D0
@@ -65,7 +65,7 @@ double complex function cdli6(z)
    if (lnz**2 + pz**2 .lt. 1) then ! |log(z)| < 1
       u = dcmplx(lnz, pz) ! log(z)
       u2 = u**2
-      c5 = (137D0/60 - fast_pos_cdlog(-u))/120
+      c5 = (137D0/60 - pos_cdlog(-u))/120
       cdli6 = zeta6 + u * c1 + &
          u2 * (c2 + u * c3 +   &
          u2 * (c4 + u * c5 +   &
@@ -79,7 +79,7 @@ double complex function cdli6(z)
    endif
 
    if (nz .le. 1) then
-      u = -fast_pos_cdlog(1 - z)
+      u = -pos_cdlog(1 - z)
       r = 0
       sgn = 1
    else ! nz > 1
@@ -89,7 +89,7 @@ double complex function cdli6(z)
          arg = pz + PI
       endif
       lmz = dcmplx(lnz, arg) ! log(-z)
-      u = -fast_pos_cdlog(1 - 1/z)
+      u = -pos_cdlog(1 - 1/z)
       r = -31*PI6/15120 + lmz**2*(-7*PI4/720 + lmz**2*(-PI2/144 - lmz**2/720))
       sgn = -1
    endif
diff --git a/src/fortran/fast_clog.f90 b/src/fortran/log.f90
similarity index 62%
rename from src/fortran/fast_clog.f90
rename to src/fortran/log.f90
index 9cc9149..325e804 100644
--- a/src/fortran/fast_clog.f90
+++ b/src/fortran/log.f90
@@ -6,30 +6,38 @@
 
 
 !*********************************************************************
-!> @brief Fast implementation of complex logarithm
+!> @brief Implementation of log(1 + z) for complex z
 !> @param z complex argument
-!> @return log(z)
+!> @return log(1 + z)
 !*********************************************************************
-double complex function fast_cdlog(z)
+double complex function cdlog1p(z)
   implicit none
-  double complex :: z
+  double complex :: z, u
   double precision :: re, im
 
-  re = real(z)
-  im = aimag(z)
-  fast_cdlog = dcmplx(log(hypot(re, im)), datan2(im, re))
+  u = 1 + z
+  re = real(u)
+  im = aimag(u)
+
+  if (re .eq. 1 .and. im .eq. 0) then
+     cdlog1p = z
+  elseif (re .le. 0) then
+     cdlog1p = log(u)
+  else
+     cdlog1p = log(u)*(z/(u - 1));
+  endif
 
-end function fast_cdlog
+end function cdlog1p
 
 
 !*********************************************************************
-!> @brief Fast implementation of complex logarithm
+!> @brief Implementation of complex logarithm
 !> @param z complex argument
 !> @return log(z)
 !> @note Points on the branch cut are treated differently from log(z):
 !> Points with Im(z) == -0D0 are mapped to Im(z) == 0D0
 !*********************************************************************
-double complex function fast_pos_cdlog(z)
+double complex function pos_cdlog(z)
   implicit none
   double complex :: z
   double precision :: re, im, arg
@@ -38,11 +46,11 @@ double complex function fast_pos_cdlog(z)
   im = aimag(z)
 
   if (im .eq. 0 .and. re .gt. 0) then
-     fast_pos_cdlog = dcmplx(log(re), 0.0D0)
+     pos_cdlog = dcmplx(log(re), 0.0D0)
   elseif (im .eq. 0) then
-     fast_pos_cdlog = dcmplx(log(-re), 3.14159265358979324D0)
+     pos_cdlog = dcmplx(log(-re), 3.14159265358979324D0)
   else
-     fast_pos_cdlog = dcmplx(log(hypot(re, im)), datan2(im, re))
+     pos_cdlog = log(z)
   endif
 
-end function fast_pos_cdlog
+end function pos_cdlog
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index bccac34..2260a94 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -48,6 +48,7 @@ add_polylogarithm_test(test_Li3)
 add_polylogarithm_test(test_Li4)
 add_polylogarithm_test(test_Li5)
 add_polylogarithm_test(test_Li6)
+add_polylogarithm_test(test_log)
 add_polylogarithm_test(test_Sl)
 add_polylogarithm_test(test_zeta)
 add_polylogarithm_test(test_version)
diff --git a/test/test.hpp b/test/test.hpp
index e388e1b..f9b5d60 100644
--- a/test/test.hpp
+++ b/test/test.hpp
@@ -33,9 +33,17 @@
 #define CHECK_SMALL(a,eps) CHECK(std::abs(a) < (eps))
 
 
+#define CHECK_CLOSE_REL(a,b,eps) do {                                   \
+      const bool pred = is_close_rel((a), (b), (eps));                  \
+      INFO("Comparing numbers " << std::setprecision(17) << (a) << " =?= " << (b) << " with relative precision " << (eps)); \
+      CHECK(pred);                                                      \
+   } while (0)
+
+
 inline bool has_inf() noexcept
 {
-   return std::isinf(std::numeric_limits<double>::max() + 1.0);
+   const auto fn = [] (bool return_inf) { return return_inf ? std::numeric_limits<double>::infinity() : 0.0; };
+   return std::isinf(fn(true));
 }
 
 
@@ -46,6 +54,22 @@ inline bool has_signed_zero() noexcept
 }
 
 
+inline bool is_ieee754_compliant() noexcept
+{
+   return has_inf() && has_signed_zero();
+}
+
+
+inline bool is_close_rel(double x, double y, double eps) noexcept
+{
+   const double ma = std::max(std::abs(x), std::abs(y));
+   if (ma == 0.0) {
+      return true;
+   }
+   return std::abs(x - y)/ma < std::abs(eps);
+}
+
+
 template <class T> T sqr(T x) {
    return x*x;
 }
diff --git a/test/test_Li.cpp b/test/test_Li.cpp
index f4fca99..7784680 100644
--- a/test/test_Li.cpp
+++ b/test/test_Li.cpp
@@ -107,6 +107,10 @@ struct Data {
 
 TEST_CASE("test_overflow")
 {
+   if (!has_inf()) {
+      return;
+   }
+
    using polylogarithm::Li;
 
    const double eps64 = std::pow(10.0, -std::numeric_limits<double>::digits10);
diff --git a/test/test_Li2.cpp b/test/test_Li2.cpp
index be712f7..ece1567 100644
--- a/test/test_Li2.cpp
+++ b/test/test_Li2.cpp
@@ -306,6 +306,22 @@ TEST_CASE("test_special_values")
       CHECK_CLOSE_COMPLEX(sherpa_Li2(z0), li0, 2*eps);
       CHECK_CLOSE_COMPLEX(spheno_Li2(z0), li0, eps);
    }
+
+   if (is_ieee754_compliant()) {
+      // special point with small real part
+      const std::complex<double> z(4.831285545908206e-6, 0.004396919500211628);
+      const std::complex<double> expected(-1.94166578202937687444628936853e-9, 0.00439692067657240512726530759719387623);
+#ifndef __APPLE__
+      CHECK_CLOSE_REL(std::real(Li2(z))             , std::real(expected), 1e-12);
+#endif
+      CHECK_CLOSE_REL(std::imag(Li2(z))             , std::imag(expected), 1e-14);
+      CHECK_CLOSE_REL(std::real(poly_Li2(z))        , std::real(expected), 1e-12);
+      CHECK_CLOSE_REL(std::imag(poly_Li2(z))        , std::imag(expected), 1e-14);
+#ifdef ENABLE_FORTRAN
+      CHECK_CLOSE_REL(std::real(poly_Li2_fortran(z)), std::real(expected), 1e-12);
+      CHECK_CLOSE_REL(std::imag(poly_Li2_fortran(z)), std::imag(expected), 1e-14);
+#endif
+   }
 }
 
 // tests signbit for 0.0 and -0.0 arguments
diff --git a/test/test_Li3.cpp b/test/test_Li3.cpp
index 79633d2..90a827c 100644
--- a/test/test_Li3.cpp
+++ b/test/test_Li3.cpp
@@ -209,6 +209,10 @@ struct Data {
 
 TEST_CASE("test_overflow")
 {
+   if (!has_inf()) {
+      return;
+   }
+
    using polylogarithm::Li3;
 
    const double eps64 = std::pow(10.0, -std::numeric_limits<double>::digits10);
diff --git a/test/test_Li4.cpp b/test/test_Li4.cpp
index c550955..c8478f5 100644
--- a/test/test_Li4.cpp
+++ b/test/test_Li4.cpp
@@ -144,6 +144,10 @@ struct Data {
 
 TEST_CASE("test_overflow")
 {
+   if (!has_inf()) {
+      return;
+   }
+
    using polylogarithm::Li4;
 
    const double eps64 = std::pow(10.0, -std::numeric_limits<double>::digits10);
diff --git a/test/test_Li5.cpp b/test/test_Li5.cpp
index c9a3c61..616aaa0 100644
--- a/test/test_Li5.cpp
+++ b/test/test_Li5.cpp
@@ -123,6 +123,10 @@ struct Data {
 
 TEST_CASE("test_overflow")
 {
+   if (!has_inf()) {
+      return;
+   }
+
    using polylogarithm::Li5;
 
    const double eps64 = std::pow(10.0, -std::numeric_limits<double>::digits10);
diff --git a/test/test_Li6.cpp b/test/test_Li6.cpp
index 6956d4a..7208c4f 100644
--- a/test/test_Li6.cpp
+++ b/test/test_Li6.cpp
@@ -123,6 +123,10 @@ struct Data {
 
 TEST_CASE("test_overflow")
 {
+   if (!has_inf()) {
+      return;
+   }
+
    using polylogarithm::Li6;
 
    const double eps64 = std::pow(10.0, -std::numeric_limits<double>::digits10);
diff --git a/test/test_log.cpp b/test/test_log.cpp
new file mode 100644
index 0000000..dff09bd
--- /dev/null
+++ b/test/test_log.cpp
@@ -0,0 +1,40 @@
+#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN 1
+
+#include "doctest.h"
+#include "log.hpp"
+#include "test.hpp"
+#include <complex>
+#include <limits>
+
+// tests log(z) and log1p(1+z) for small complex z
+TEST_CASE("test_log_and_log1p")
+{
+   if (!is_ieee754_compliant()) {
+      return;
+   }
+
+   using namespace polylogarithm;
+
+   const double eps = std::numeric_limits<double>::epsilon();
+   // z = 9007155738389423 / 2^53 - (5069303045819176 / 2^60) I
+   const std::complex<double> z(0.999995168714454, -0.004396919500211628);
+   const std::complex<double> cpp_log(std::log(z));
+   const std::complex<double> pl_log = pos_log(z);
+   const std::complex<double> pl_log1p = log1p(z);
+
+   const std::complex<double> expected_log(4.8351532915892516848328240700759518475389916e-6, 0.00439691240793594643979611108925073970426108463713022280825);
+   const std::complex<double> expected_log1p(0.693147181532726411790714997046464927858457, -0.002198461518912911475356949677977809228317);
+
+   INFO("z             = " << std::setprecision(17) << z);
+   INFO("std::log(z)   = " << cpp_log);
+   INFO("pos_log(z)    = " << std::complex<double>(pl_log));
+   INFO("exp. log(z)   = " << expected_log);
+   INFO("log1p(z)      = " << std::complex<double>(pl_log1p));
+   INFO("exp. log1p(z) = " << expected_log1p);
+
+#ifndef __APPLE__
+   CHECK_CLOSE_REL(std::real(cpp_log), std::real(expected_log), eps);
+   CHECK_CLOSE_REL(std::real(pl_log), std::real(expected_log), eps);
+#endif
+   CHECK_CLOSE_REL(std::real(pl_log1p), std::real(expected_log1p), eps);
+}