+
+/*static ::cl_N Li2_series(const ::cl_N & x,
+ const ::cl_float_format_t & prec)
+{
+ // Note: argument must be in the unit circle
+ // This is very inefficient unless we have fast floating point Bernoulli
+ // numbers implemented!
+ ::cl_N c1 = -::log(1-x);
+ ::cl_N c2 = c1;
+ // hard-wire the first two Bernoulli numbers
+ ::cl_N acc = c1 - ::square(c1)/4;
+ ::cl_N aug;
+ ::cl_F pisq = ::square(::cl_pi(prec)); // pi^2
+ ::cl_F piac = ::cl_float(1, prec); // accumulator: pi^(2*i)
+ unsigned i = 1;
+ c1 = ::square(c1);
+ do {
+ c2 = c1 * c2;
+ piac = piac * pisq;
+ aug = c2 * (*(bernoulli(numeric(2*i)).clnptr())) / ::factorial(2*i+1);
+ // aug = c2 * ::cl_I(i%2 ? 1 : -1) / ::cl_I(2*i+1) * ::cl_zeta(2*i, prec) / piac / (::cl_I(1)<<(2*i-1));
+ acc = acc + aug;
+ ++i;
+ } while (acc != acc+aug);
+ return acc;
+}*/
+
+/** Numeric evaluation of Dilogarithm within circle of convergence (unit
+ * circle) using a power series. */
+static ::cl_N Li2_series(const ::cl_N & x,
+ const ::cl_float_format_t & prec)
+{
+ // Note: argument must be in the unit circle
+ ::cl_N aug, acc;
+ ::cl_N num = ::complex(::cl_float(1, prec), 0);
+ ::cl_I den = 0;
+ unsigned i = 1;
+ do {
+ num = num * x;
+ den = den + i; // 1, 4, 9, 16, ...
+ i += 2;
+ aug = num / den;
+ acc = acc + aug;
+ } while (acc != acc+aug);
+ return acc;
+}
+
+/** Folds Li2's argument inside a small rectangle to enhance convergence. */
+static ::cl_N Li2_projection(const ::cl_N & x,
+ const ::cl_float_format_t & prec)
+{
+ const ::cl_R re = ::realpart(x);
+ const ::cl_R im = ::imagpart(x);
+ if (re > ::cl_F(".5"))
+ // zeta(2) - Li2(1-x) - log(x)*log(1-x)
+ return(::cl_zeta(2)
+ - Li2_series(1-x, prec)
+ - ::log(x)*::log(1-x));
+ if ((re <= 0 && ::abs(im) > ::cl_F(".75")) || (re < ::cl_F("-.5")))
+ // -log(1-x)^2 / 2 - Li2(x/(x-1))
+ return(-::square(::log(1-x))/2
+ - Li2_series(x/(x-1), prec));
+ if (re > 0 && ::abs(im) > ::cl_LF(".75"))
+ // Li2(x^2)/2 - Li2(-x)
+ return(Li2_projection(::square(x), prec)/2
+ - Li2_projection(-x, prec));
+ return Li2_series(x, prec);
+}
+
+/** Numeric evaluation of Dilogarithm. The domain is the entire complex plane,
+ * the branch cut lies along the positive real axis, starting at 1 and
+ * continuous with quadrant IV.
+ *
+ * @return arbitrary precision numerical Li2(x). */
+const numeric Li2(const numeric & x)