+/*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)
+{
+ if (::zerop(*x.value))
+ return x;
+
+ // what is the desired float format?
+ // first guess: default format
+ ::cl_float_format_t prec = ::cl_default_float_format;
+ // second guess: the argument's format
+ if (!::instanceof(::realpart(*x.value),cl_RA_ring))
+ prec = ::cl_float_format(The(::cl_F)(::realpart(*x.value)));
+ else if (!::instanceof(::imagpart(*x.value),cl_RA_ring))
+ prec = ::cl_float_format(The(::cl_F)(::imagpart(*x.value)));
+
+ if (*x.value==1) // may cause trouble with log(1-x)
+ return ::cl_zeta(2, prec);
+
+ if (::abs(*x.value) > 1)
+ // -log(-x)^2 / 2 - zeta(2) - Li2(1/x)
+ return(-::square(::log(-*x.value))/2
+ - ::cl_zeta(2, prec)
+ - Li2_projection(::recip(*x.value), prec));
+ else
+ return Li2_projection(*x.value, prec);
+}
+
+