// forward declaration needed by function Li_projection and C below
-numeric S_num(int n, int p, const numeric& x);
+const cln::cl_N S_num(int n, int p, const cln::cl_N& x);
// helper function for classical polylog Li
} else {
cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
for (int j=0; j<n-1; j++) {
- result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
+ result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
* cln::expt(cln::log(x), j) / cln::factorial(j);
}
return result;
}
}
-
// helper function for classical polylog Li
-numeric Lin_numeric(int n, const numeric& x)
+const cln::cl_N Lin_numeric(const int n, const cln::cl_N& x)
{
if (n == 1) {
// just a log
- return -cln::log(1-x.to_cl_N());
+ return -cln::log(1-x);
}
- if (x.is_zero()) {
+ if (zerop(x)) {
return 0;
}
if (x == 1) {
// [Kol] (2.22)
return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
}
- if (abs(x.real()) < 0.4 && abs(abs(x)-1) < 0.01) {
- cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
- cln::cl_N result = -cln::expt(cln::log(x_), n-1) * cln::log(1-x_) / cln::factorial(n-1);
+ if (cln::abs(realpart(x)) < 0.4 && cln::abs(cln::abs(x)-1) < 0.01) {
+ cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
for (int j=0; j<n-1; j++) {
- result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x_).to_cl_N())
- * cln::expt(cln::log(x_), j) / cln::factorial(j);
+ result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
+ * cln::expt(cln::log(x), j) / cln::factorial(j);
}
return result;
}
// 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();
+ const cln::cl_N value = x;
// second guess: the argument's format
- if (!x.real().is_rational())
+ if (!instanceof(realpart(x), cln::cl_RA_ring))
prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
- else if (!x.imag().is_rational())
+ else if (!instanceof(imagpart(x), cln::cl_RA_ring))
prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
// [Kol] (5.15)
cln::cl_N add;
for (int j=0; j<n-1; j++) {
add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
- * Lin_numeric(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
+ * Lin_numeric(n-j,1) * cln::expt(cln::log(-value),j) / cln::factorial(j);
}
result = result - add;
return result;
t0buf = t[0];
q++;
t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
+ for (int k=j-2; k>=0; k--) {
+ t[k] = t[k] + t[k+1] * cln::expt(x[k], q+j-1-k) / cln::expt(cln::cl_I(q+j-1-k), s[k]);
+ }
+ q++;
+ t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
for (int k=j-2; k>=0; k--) {
flag_accidental_zero = cln::zerop(t[k+1]);
t[k] = t[k] + t[k+1] * cln::expt(x[k], q+j-1-k) / cln::expt(cln::cl_I(q+j-1-k), s[k]);
}
- } while ( (t[0] != t0buf) || flag_accidental_zero );
+ } while ( (t[0] != t0buf) || cln::zerop(t[0]) || flag_accidental_zero );
return t[0];
}
if (x.op(x.nops()-1).is_zero()) {
need_trafo = true;
}
- if (depth == 1 && !need_trafo) {
+ if (depth == 1 && x.nops() == 2 && !need_trafo) {
return -Li(x.nops(), y / x.op(x.nops()-1)).evalf();
}
// classical polylogs
if (m_.info(info_flags::posint)) {
if (x_.info(info_flags::numeric)) {
- return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_));
+ int m__ = ex_to<numeric>(m_).to_int();
+ const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
+ const cln::cl_N result = Lin_numeric(m__, x__);
+ return numeric(result);
} else {
// try to numerically evaluate second argument
ex x_val = x_.evalf();
if (x_val.info(info_flags::numeric)) {
- return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_val));
+ int m__ = ex_to<numeric>(m_).to_int();
+ const cln::cl_N x__ = ex_to<numeric>(x_val).to_cl_N();
+ const cln::cl_N result = Lin_numeric(m__, x__);
+ return numeric(result);
}
}
}
}
}
if (m_.info(info_flags::posint) && x_.info(info_flags::numeric) && !x_.info(info_flags::crational)) {
- return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_));
+ int m__ = ex_to<numeric>(m_).to_int();
+ const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
+ const cln::cl_N result = Lin_numeric(m__, x__);
+ return numeric(result);
}
return Li(m_, x_).hold();
if (k == 0) {
if (n & 1) {
if (j & 1) {
- result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
+ result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
}
else {
- result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
+ result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
}
}
}
if (k & 1) {
if (j & 1) {
result = result + cln::factorial(n+k-1)
- * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+ * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
/ (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
}
else {
result = result - cln::factorial(n+k-1)
- * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+ * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
/ (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
}
}
else {
if (j & 1) {
- result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+ result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
/ (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
}
else {
result = result + cln::factorial(n+k-1)
- * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+ * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
/ (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
}
}
// helper function for S(n,p,x)
cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
{
+ static cln::float_format_t oldprec = cln::default_float_format;
+
if (p==1) {
return Li_projection(n+1, x, prec);
}
-
+
+ // precision has changed, we need to clear lookup table Yn
+ if ( oldprec != prec ) {
+ Yn.clear();
+ ynsize = 0;
+ ynlength = 100;
+ oldprec = prec;
+ }
+
// check if precalculated values are sufficient
if (p > ynsize+1) {
for (int i=ynsize; i<p-1; i++) {
res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
* S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
}
- result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
+ result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
}
return result;
// helper function for S(n,p,x)
-numeric S_num(int n, int p, const numeric& x)
+const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
{
if (x == 1) {
if (n == 1) {
// 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();
+ const cln::cl_N value = x;
// second guess: the argument's format
- if (!x.real().is_rational())
+ if (!instanceof(realpart(value), cln::cl_RA_ring))
prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
- else if (!x.imag().is_rational())
+ else if (!instanceof(imagpart(value), cln::cl_RA_ring))
prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
// [Kol] (5.3)
cln::cl_N res2;
for (int r=0; r<p; r++) {
res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
- * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
+ * S_num(p-r,n-s,1-value) / cln::factorial(r);
}
- result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
+ result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
}
return result;
for (int r=0; r<=s; r++) {
result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
/ cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
- * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
+ * S_num(n+s-r,p-s,cln::recip(value));
}
}
result = result * cln::expt(cln::cl_I(-1),n);
static ex S_evalf(const ex& n, const ex& p, const ex& x)
{
if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
+ const int n_ = ex_to<numeric>(n).to_int();
+ const int p_ = ex_to<numeric>(p).to_int();
if (is_a<numeric>(x)) {
- return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
+ const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
+ const cln::cl_N result = S_num(n_, p_, x_);
+ return numeric(result);
} else {
ex x_val = x.evalf();
if (is_a<numeric>(x_val)) {
- return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x_val));
+ const cln::cl_N x_val_ = ex_to<numeric>(x_val).to_cl_N();
+ const cln::cl_N result = S_num(n_, p_, x_val_);
+ return numeric(result);
}
}
}
return Li(n+1, x);
}
if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
- return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
+ int n_ = ex_to<numeric>(n).to_int();
+ int p_ = ex_to<numeric>(p).to_int();
+ const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
+ const cln::cl_N result = S_num(n_, p_, x_);
+ return numeric(result);
}
}
if (n.is_zero()) {