+const numeric atanh(const numeric &x)
+{
+ return cln::atanh(x.to_cl_N());
+}
+
+
+/*static cln::cl_N Li2_series(const ::cl_N &x,
+ const ::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!
+ cln::cl_N c1 = -cln::log(1-x);
+ cln::cl_N c2 = c1;
+ // hard-wire the first two Bernoulli numbers
+ cln::cl_N acc = c1 - cln::square(c1)/4;
+ cln::cl_N aug;
+ cln::cl_F pisq = cln::square(cln::cl_pi(prec)); // pi^2
+ cln::cl_F piac = cln::cl_float(1, prec); // accumulator: pi^(2*i)
+ unsigned i = 1;
+ c1 = cln::square(c1);
+ do {
+ c2 = c1 * c2;
+ piac = piac * pisq;
+ aug = c2 * (*(bernoulli(numeric(2*i)).clnptr())) / cln::factorial(2*i+1);
+ // aug = c2 * cln::cl_I(i%2 ? 1 : -1) / cln::cl_I(2*i+1) * cln::cl_zeta(2*i, prec) / piac / (cln::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 cln::cl_N Li2_series(const cln::cl_N &x,
+ const cln::float_format_t &prec)
+{
+ // Note: argument must be in the unit circle
+ cln::cl_N aug, acc;
+ cln::cl_N num = cln::complex(cln::cl_float(1, prec), 0);
+ cln::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 cln::cl_N Li2_projection(const cln::cl_N &x,
+ const cln::float_format_t &prec)
+{
+ const cln::cl_R re = cln::realpart(x);
+ const cln::cl_R im = cln::imagpart(x);
+ if (re > cln::cl_F(".5"))
+ // zeta(2) - Li2(1-x) - log(x)*log(1-x)
+ return(cln::zeta(2)
+ - Li2_series(1-x, prec)
+ - cln::log(x)*cln::log(1-x));
+ if ((re <= 0 && cln::abs(im) > cln::cl_F(".75")) || (re < cln::cl_F("-.5")))
+ // -log(1-x)^2 / 2 - Li2(x/(x-1))
+ return(- cln::square(cln::log(1-x))/2
+ - Li2_series(x/(x-1), prec));
+ if (re > 0 && cln::abs(im) > cln::cl_LF(".75"))
+ // Li2(x^2)/2 - Li2(-x)
+ return(Li2_projection(cln::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 (x.is_zero())
+ return _num0;
+
+ // what is the desired float format?
+ // first guess: default format
+ cln::float_format_t prec = cln::default_float_format;
+ const cln::cl_N value = x.to_cl_N();
+ // second guess: the argument's format
+ if (!x.real().is_rational())
+ prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
+ else if (!x.imag().is_rational())
+ prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
+
+ if (cln::the<cln::cl_N>(value)==1) // may cause trouble with log(1-x)
+ return cln::zeta(2, prec);
+
+ if (cln::abs(value) > 1)
+ // -log(-x)^2 / 2 - zeta(2) - Li2(1/x)
+ return(- cln::square(cln::log(-value))/2
+ - cln::zeta(2, prec)
+ - Li2_projection(cln::recip(value), prec));
+ else
+ return Li2_projection(x.to_cl_N(), prec);