+ ex operator()(const ex& e)
+ {
+ if (is_a<add>(e) || is_a<mul>(e)) {
+ return e.map(*this);
+ }
+ if (is_a<function>(e)) {
+ std::string name = ex_to<function>(e).get_name();
+ if (name == "H") {
+ lst parameter;
+ if (is_a<lst>(e.op(0))) {
+ parameter = ex_to<lst>(e.op(0));
+ } else {
+ parameter = lst(e.op(0));
+ }
+
+ lst m;
+ lst s;
+ ex pf;
+ if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
+ return pf * zeta(m, s);
+ } else {
+ return zeta(m);
+ }
+ }
+ }
+ return e;
+ }
+};
+
+
+// remove trailing zeros from H-parameters
+struct map_trafo_H_reduce_trailing_zeros : public map_function
+{
+ ex operator()(const ex& e)
+ {
+ if (is_a<add>(e) || is_a<mul>(e)) {
+ return e.map(*this);
+ }
+ if (is_a<function>(e)) {
+ std::string name = ex_to<function>(e).get_name();
+ if (name == "H") {
+ lst parameter;
+ if (is_a<lst>(e.op(0))) {
+ parameter = ex_to<lst>(e.op(0));
+ } else {
+ parameter = lst(e.op(0));
+ }
+ ex arg = e.op(1);
+ if (parameter.op(parameter.nops()-1) == 0) {
+
+ //
+ if (parameter.nops() == 1) {
+ return log(arg);
+ }
+
+ //
+ lst::const_iterator it = parameter.begin();
+ while ((it != parameter.end()) && (*it == 0)) {
+ it++;
+ }
+ if (it == parameter.end()) {
+ return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
+ }
+
+ //
+ parameter.remove_last();
+ int lastentry = parameter.nops();
+ while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
+ lastentry--;
+ }
+
+ //
+ ex result = log(arg) * H(parameter,arg).hold();
+ ex acc = 0;
+ for (ex i=0; i<lastentry; i++) {
+ if (parameter[i] > 0) {
+ parameter[i]++;
+ result -= (acc + parameter[i]-1) * H(parameter, arg).hold();
+ parameter[i]--;
+ acc = 0;
+ } else if (parameter[i] < 0) {
+ parameter[i]--;
+ result -= (acc + abs(parameter[i]+1)) * H(parameter, arg).hold();
+ parameter[i]++;
+ acc = 0;
+ } else {
+ acc++;
+ }
+ }
+
+ if (lastentry < parameter.nops()) {
+ result = result / (parameter.nops()-lastentry+1);
+ return result.map(*this);
+ } else {
+ return result;
+ }
+ }
+ }
+ }
+ return e;
+ }
+};
+
+
+// returns an expression with zeta functions corresponding to the parameter list for H
+ex convert_H_to_zeta(const lst& m)
+{
+ symbol xtemp("xtemp");
+ map_trafo_H_reduce_trailing_zeros filter;
+ map_trafo_H_convert_to_zeta filter2;
+ return filter2(filter(H(m, xtemp).hold())).subs(xtemp == 1);
+}
+
+
+// convert signs form Li to H representation
+lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf)
+{
+ lst res;
+ lst::const_iterator itm = m.begin();
+ lst::const_iterator itx = ++x.begin();
+ ex signum = _ex1;
+ pf = _ex1;
+ res.append(*itm);
+ itm++;
+ while (itx != x.end()) {
+ signum *= *itx;
+ pf *= signum;
+ res.append((*itm) * signum);
+ itm++;
+ itx++;
+ }
+ return res;
+}
+
+
+// multiplies an one-dimensional H with another H
+// [ReV] (18)
+ex trafo_H_mult(const ex& h1, const ex& h2)
+{
+ ex res;
+ ex hshort;
+ lst hlong;
+ ex h1nops = h1.op(0).nops();
+ ex h2nops = h2.op(0).nops();
+ if (h1nops > 1) {
+ hshort = h2.op(0).op(0);
+ hlong = ex_to<lst>(h1.op(0));
+ } else {
+ hshort = h1.op(0).op(0);
+ if (h2nops > 1) {
+ hlong = ex_to<lst>(h2.op(0));
+ } else {
+ hlong = h2.op(0).op(0);
+ }
+ }
+ for (int i=0; i<=hlong.nops(); i++) {
+ lst newparameter;
+ int j=0;
+ for (; j<i; j++) {
+ newparameter.append(hlong[j]);
+ }
+ newparameter.append(hshort);
+ for (; j<hlong.nops(); j++) {
+ newparameter.append(hlong[j]);
+ }
+ res += H(newparameter, h1.op(1)).hold();
+ }
+ return res;
+}
+
+
+// applies trafo_H_mult recursively on expressions
+struct map_trafo_H_mult : public map_function
+{
+ ex operator()(const ex& e)
+ {
+ if (is_a<add>(e)) {
+ return e.map(*this);
+ }
+
+ if (is_a<mul>(e)) {
+
+ ex result = 1;
+ ex firstH;
+ lst Hlst;
+ for (int pos=0; pos<e.nops(); pos++) {
+ if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
+ std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
+ if (name == "H") {
+ for (ex i=0; i<e.op(pos).op(1); i++) {
+ Hlst.append(e.op(pos).op(0));
+ }
+ continue;
+ }
+ } else if (is_a<function>(e.op(pos))) {
+ std::string name = ex_to<function>(e.op(pos)).get_name();
+ if (name == "H") {
+ if (e.op(pos).op(0).nops() > 1) {
+ firstH = e.op(pos);
+ } else {
+ Hlst.append(e.op(pos));
+ }
+ continue;
+ }
+ }
+ result *= e.op(pos);
+ }
+ if (firstH == 0) {
+ if (Hlst.nops() > 0) {
+ firstH = Hlst[Hlst.nops()-1];
+ Hlst.remove_last();
+ } else {
+ return e;
+ }
+ }
+
+ if (Hlst.nops() > 0) {
+ ex buffer = trafo_H_mult(firstH, Hlst.op(0));
+ result *= buffer;
+ for (int i=1; i<Hlst.nops(); i++) {
+ result *= Hlst.op(i);
+ }
+ result = result.expand();
+ map_trafo_H_mult recursion;
+ return recursion(result);
+ } else {
+ return e;
+ }
+
+ }
+ return e;
+ }
+};
+
+
+// do integration [ReV] (55)
+// put parameter 0 in front of existing parameters
+ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
+{
+ ex h;
+ std::string name;
+ if (is_a<function>(e)) {
+ name = ex_to<function>(e).get_name();
+ }
+ if (name == "H") {
+ h = e;
+ } else {
+ for (int i=0; i<e.nops(); i++) {
+ if (is_a<function>(e.op(i))) {
+ std::string name = ex_to<function>(e.op(i)).get_name();
+ if (name == "H") {
+ h = e.op(i);
+ }
+ }
+ }
+ }
+ if (h != 0) {
+ lst newparameter = ex_to<lst>(h.op(0));
+ newparameter.prepend(0);
+ ex addzeta = convert_H_to_zeta(newparameter);
+ return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
+ } else {
+ return e * (-H(lst(0),1/arg).hold());
+ }
+}
+
+
+// do integration [ReV] (55)
+// put parameter -1 in front of existing parameters
+ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
+{
+ ex h;
+ std::string name;
+ if (is_a<function>(e)) {
+ name = ex_to<function>(e).get_name();
+ }
+ if (name == "H") {
+ h = e;
+ } else {
+ for (int i=0; i<e.nops(); i++) {
+ if (is_a<function>(e.op(i))) {
+ std::string name = ex_to<function>(e.op(i)).get_name();
+ if (name == "H") {
+ h = e.op(i);
+ }
+ }
+ }
+ }
+ if (h != 0) {
+ lst newparameter = ex_to<lst>(h.op(0));
+ newparameter.prepend(-1);
+ ex addzeta = convert_H_to_zeta(newparameter);
+ return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
+ } else {
+ ex addzeta = convert_H_to_zeta(lst(-1));
+ return (e * (addzeta - H(lst(-1),1/arg).hold())).expand();
+ }
+}
+
+
+// do integration [ReV] (55)
+// put parameter -1 in front of existing parameters
+ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
+{
+ ex h;
+ std::string name;
+ if (is_a<function>(e)) {
+ name = ex_to<function>(e).get_name();
+ }
+ if (name == "H") {
+ h = e;
+ } else {
+ for (int i=0; i<e.nops(); i++) {
+ if (is_a<function>(e.op(i))) {
+ std::string name = ex_to<function>(e.op(i)).get_name();
+ if (name == "H") {
+ h = e.op(i);
+ }
+ }
+ }
+ }
+ if (h != 0) {
+ lst newparameter = ex_to<lst>(h.op(0));
+ newparameter.prepend(-1);
+ return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
+ } else {
+ return (e * H(lst(-1),(1-arg)/(1+arg)).hold()).expand();
+ }
+}
+
+
+// do integration [ReV] (55)
+// put parameter 1 in front of existing parameters
+ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
+{
+ ex h;
+ std::string name;
+ if (is_a<function>(e)) {
+ name = ex_to<function>(e).get_name();
+ }
+ if (name == "H") {
+ h = e;
+ } else {
+ for (int i=0; i<e.nops(); i++) {
+ if (is_a<function>(e.op(i))) {
+ std::string name = ex_to<function>(e.op(i)).get_name();
+ if (name == "H") {
+ h = e.op(i);
+ }
+ }
+ }
+ }
+ if (h != 0) {
+ lst newparameter = ex_to<lst>(h.op(0));
+ newparameter.prepend(1);
+ return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
+ } else {
+ return (e * H(lst(1),(1-arg)/(1+arg)).hold()).expand();
+ }
+}
+
+
+// do x -> 1/x transformation
+struct map_trafo_H_1overx : public map_function
+{
+ ex operator()(const ex& e)
+ {
+ if (is_a<add>(e) || is_a<mul>(e)) {
+ return e.map(*this);
+ }
+
+ if (is_a<function>(e)) {
+ std::string name = ex_to<function>(e).get_name();
+ if (name == "H") {
+
+ lst parameter = ex_to<lst>(e.op(0));
+ ex arg = e.op(1);
+
+ // special cases if all parameters are either 0, 1 or -1
+ bool allthesame = true;
+ if (parameter.op(0) == 0) {
+ for (int i=1; i<parameter.nops(); i++) {
+ if (parameter.op(i) != 0) {
+ allthesame = false;
+ break;
+ }
+ }
+ if (allthesame) {
+ return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
+ }
+ } else if (parameter.op(0) == -1) {
+ for (int i=1; i<parameter.nops(); i++) {
+ if (parameter.op(i) != -1) {
+ allthesame = false;
+ break;
+ }
+ }
+ if (allthesame) {
+ map_trafo_H_mult unify;
+ return unify((pow(H(lst(-1),1/arg).hold() - H(lst(0),1/arg).hold(), parameter.nops())
+ / factorial(parameter.nops())).expand());
+ }
+ } else {
+ for (int i=1; i<parameter.nops(); i++) {
+ if (parameter.op(i) != 1) {
+ allthesame = false;
+ break;
+ }
+ }
+ if (allthesame) {
+ map_trafo_H_mult unify;
+ return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() - I*Pi, parameter.nops())
+ / factorial(parameter.nops())).expand());
+ }
+ }
+
+ lst newparameter = parameter;
+ newparameter.remove_first();
+
+ if (parameter.op(0) == 0) {
+
+ // leading zero
+ ex res = convert_H_to_zeta(parameter);
+ map_trafo_H_1overx recursion;
+ ex buffer = recursion(H(newparameter, arg).hold());
+ if (is_a<add>(buffer)) {
+ for (int i=0; i<buffer.nops(); i++) {
+ res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
+ }
+ } else {
+ res += trafo_H_1tx_prepend_zero(buffer, arg);
+ }
+ return res;
+
+ } else if (parameter.op(0) == -1) {
+
+ // leading negative one
+ ex res = convert_H_to_zeta(parameter);
+ map_trafo_H_1overx recursion;
+ ex buffer = recursion(H(newparameter, arg).hold());
+ if (is_a<add>(buffer)) {
+ for (int i=0; i<buffer.nops(); i++) {
+ res += trafo_H_1tx_prepend_zero(buffer.op(i), arg) - trafo_H_1tx_prepend_minusone(buffer.op(i), arg);
+ }
+ } else {
+ res += trafo_H_1tx_prepend_zero(buffer, arg) - trafo_H_1tx_prepend_minusone(buffer, arg);
+ }
+ return res;
+
+ } else {
+
+ // leading one
+ map_trafo_H_1overx recursion;
+ map_trafo_H_mult unify;
+ ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
+ int firstzero = 0;
+ while (parameter.op(firstzero) == 1) {
+ firstzero++;
+ }
+ for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+ lst newparameter;
+ int j=0;
+ for (; j<=i; j++) {
+ newparameter.append(parameter[j+1]);
+ }
+ newparameter.append(1);
+ for (; j<parameter.nops()-1; j++) {
+ newparameter.append(parameter[j+1]);
+ }
+ res -= H(newparameter, arg).hold();
+ }
+ res = recursion(res).expand() / firstzero;
+ return unify(res);
+
+ }
+
+ }
+ }
+ return e;
+ }
+};
+
+
+// do x -> (1-x)/(1+x) transformation
+struct map_trafo_H_1mxt1px : public map_function
+{
+ ex operator()(const ex& e)
+ {
+ if (is_a<add>(e) || is_a<mul>(e)) {
+ return e.map(*this);
+ }
+
+ if (is_a<function>(e)) {
+ std::string name = ex_to<function>(e).get_name();
+ if (name == "H") {
+
+ lst parameter = ex_to<lst>(e.op(0));
+ ex arg = e.op(1);
+
+ // special cases if all parameters are either 0, 1 or -1
+ bool allthesame = true;
+ if (parameter.op(0) == 0) {
+ for (int i=1; i<parameter.nops(); i++) {
+ if (parameter.op(i) != 0) {
+ allthesame = false;
+ break;
+ }
+ }
+ if (allthesame) {
+ map_trafo_H_mult unify;
+ return unify((pow(-H(lst(1),(1-arg)/(1+arg)).hold() - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+ / factorial(parameter.nops())).expand());
+ }
+ } else if (parameter.op(0) == -1) {
+ for (int i=1; i<parameter.nops(); i++) {
+ if (parameter.op(i) != -1) {
+ allthesame = false;
+ break;
+ }
+ }
+ if (allthesame) {
+ map_trafo_H_mult unify;
+ return unify((pow(log(2) - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+ / factorial(parameter.nops())).expand());
+ }
+ } else {
+ for (int i=1; i<parameter.nops(); i++) {
+ if (parameter.op(i) != 1) {
+ allthesame = false;
+ break;
+ }
+ }
+ if (allthesame) {
+ map_trafo_H_mult unify;
+ return unify((pow(-log(2) - H(lst(0),(1-arg)/(1+arg)).hold() + H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+ / factorial(parameter.nops())).expand());
+ }
+ }
+
+ lst newparameter = parameter;
+ newparameter.remove_first();
+
+ if (parameter.op(0) == 0) {
+
+ // leading zero
+ ex res = convert_H_to_zeta(parameter);
+ map_trafo_H_1mxt1px recursion;
+ ex buffer = recursion(H(newparameter, arg).hold());
+ if (is_a<add>(buffer)) {
+ for (int i=0; i<buffer.nops(); i++) {
+ res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
+ }
+ } else {
+ res -= trafo_H_1mxt1px_prepend_one(buffer, arg) + trafo_H_1mxt1px_prepend_minusone(buffer, arg);
+ }
+ return res;
+
+ } else if (parameter.op(0) == -1) {
+
+ // leading negative one
+ ex res = convert_H_to_zeta(parameter);
+ map_trafo_H_1mxt1px recursion;
+ ex buffer = recursion(H(newparameter, arg).hold());
+ if (is_a<add>(buffer)) {
+ for (int i=0; i<buffer.nops(); i++) {
+ res -= trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
+ }
+ } else {
+ res -= trafo_H_1mxt1px_prepend_minusone(buffer, arg);
+ }
+ return res;
+
+ } else {
+
+ // leading one
+ map_trafo_H_1mxt1px recursion;
+ map_trafo_H_mult unify;
+ ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
+ int firstzero = 0;
+ while (parameter.op(firstzero) == 1) {
+ firstzero++;
+ }
+ for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+ lst newparameter;
+ int j=0;
+ for (; j<=i; j++) {
+ newparameter.append(parameter[j+1]);
+ }
+ newparameter.append(1);
+ for (; j<parameter.nops()-1; j++) {
+ newparameter.append(parameter[j+1]);
+ }
+ res -= H(newparameter, arg).hold();
+ }
+ res = recursion(res).expand() / firstzero;
+ return unify(res);
+
+ }
+
+ }
+ }
+ return e;
+ }
+};
+
+
+// do the actual summation.
+cln::cl_N H_do_sum(const std::vector<int>& m, const cln::cl_N& x)
+{
+ const int j = m.size();
+
+ std::vector<cln::cl_N> t(j);
+
+ cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
+ cln::cl_N factor = cln::expt(x, j) * one;
+ cln::cl_N t0buf;
+ int q = 0;
+ do {
+ t0buf = t[0];
+ q++;
+ t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),m[j-1]);
+ for (int k=j-2; k>=1; k--) {
+ t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), m[k]);
+ }
+ t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), m[0]);
+ factor = factor * x;
+ } while (t[0] != t0buf);
+
+ return t[0];
+}
+
+
+} // end of anonymous namespace
+
+
+//////////////////////////////////////////////////////////////////////
+//
+// Harmonic polylogarithm H(m,x)
+//
+// GiNaC function
+//
+//////////////////////////////////////////////////////////////////////
+