*
* The functions are:
* classical polylogarithm Li(n,x)
- * multiple polylogarithm Li(lst(m_1,...,m_k),lst(x_1,...,x_k))
- * G(lst(a_1,...,a_k),y) or G(lst(a_1,...,a_k),lst(s_1,...,s_k),y)
+ * multiple polylogarithm Li(lst{m_1,...,m_k},lst{x_1,...,x_k})
+ * G(lst{a_1,...,a_k},y) or G(lst{a_1,...,a_k},lst{s_1,...,s_k},y)
* Nielsen's generalized polylogarithm S(n,p,x)
- * harmonic polylogarithm H(m,x) or H(lst(m_1,...,m_k),x)
- * multiple zeta value zeta(m) or zeta(lst(m_1,...,m_k))
- * alternating Euler sum zeta(m,s) or zeta(lst(m_1,...,m_k),lst(s_1,...,s_k))
+ * harmonic polylogarithm H(m,x) or H(lst{m_1,...,m_k},x)
+ * multiple zeta value zeta(m) or zeta(lst{m_1,...,m_k})
+ * alternating Euler sum zeta(m,s) or zeta(lst{m_1,...,m_k},lst{s_1,...,s_k})
*
* Some remarks:
*
* 0, 1 and -1 --- or in compactified --- a string with zeros in front of 1 or -1 is written as a single
* number --- notation.
*
- * - All functions can be nummerically evaluated with arguments in the whole complex plane. The parameters
+ * - All functions can be numerically evaluated with arguments in the whole complex plane. The parameters
* for Li, zeta and S must be positive integers. If you want to have an alternating Euler sum, you have
* to give the signs of the parameters as a second argument s to zeta(m,s) containing 1 and -1.
*
*/
/*
- * GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
// lookup table for factors built from Bernoulli numbers
// see fill_Xn()
-std::vector<std::vector<cln::cl_N> > Xn;
+std::vector<std::vector<cln::cl_N>> Xn;
// initial size of Xn that should suffice for 32bit machines (must be even)
const int xninitsizestep = 26;
int xninitsize = xninitsizestep;
return result;
}
+class less_object_for_cl_N
+{
+public:
+ bool operator() (const cln::cl_N & a, const cln::cl_N & b) const
+ {
+ // absolute value?
+ if (abs(a) != abs(b))
+ return (abs(a) < abs(b)) ? true : false;
+
+ // complex phase?
+ if (phase(a) != phase(b))
+ return (phase(a) < phase(b)) ? true : false;
+
+ // equal, therefore "less" is not true
+ return false;
+ }
+};
+
+
// convergence transformation, used for numerical evaluation of G function.
// the parameter x, s and y must only contain numerics
static cln::cl_N
const cln::cl_N& y, bool flag_trailing_zeros_only)
{
// sort (|x|<->position) to determine indices
- typedef std::multimap<cln::cl_R, std::size_t> sortmap_t;
+ typedef std::multimap<cln::cl_N, std::size_t, less_object_for_cl_N> sortmap_t;
sortmap_t sortmap;
std::size_t size = 0;
for (std::size_t i = 0; i < x.size(); ++i) {
if (!zerop(x[i])) {
- sortmap.insert(std::make_pair(abs(x[i]), i));
+ sortmap.insert(std::make_pair(x[i], i));
++size;
}
}
// include upper limit (scale)
- sortmap.insert(std::make_pair(abs(y), x.size()));
+ sortmap.insert(std::make_pair(y, x.size()));
// generate missing dummy-symbols
int i = 1;
Gparameter pendint;
ex result = G_transform(pendint, a, scale, gsyms, flag_trailing_zeros_only);
// replace dummy symbols with their values
- result = result.eval().expand();
+ result = result.expand();
result = result.subs(subslst).evalf();
if (!is_a<numeric>(result))
throw std::logic_error("G_do_trafo: G_transform returned non-numeric result");
static ex G2_evalf(const ex& x_, const ex& y)
{
- if (!y.info(info_flags::positive)) {
+ if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
return G(x_, y).hold();
}
- lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
+ lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
if (x.nops() == 0) {
return _ex1;
}
{
//TODO eval to MZV or H or S or Lin
- if (!y.info(info_flags::positive)) {
+ if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
return G(x_, y).hold();
}
- lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
+ lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
if (x.nops() == 0) {
return _ex1;
}
}
+// option do_not_evalf_params() removed.
unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2).
evalf_func(G2_evalf).
eval_func(G2_eval).
- do_not_evalf_params().
overloaded(2));
//TODO
// derivative_func(G2_deriv).
static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
{
- if (!y.info(info_flags::positive)) {
+ if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
return G(x_, s_, y).hold();
}
- lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
- lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst(s_);
+ lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
+ lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst{s_};
if (x.nops() != s.nops()) {
return G(x_, s_, y).hold();
}
{
//TODO eval to MZV or H or S or Lin
- if (!y.info(info_flags::positive)) {
+ if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
return G(x_, s_, y).hold();
}
- lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
- lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst(s_);
+ lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
+ lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst{s_};
if (x.nops() != s.nops()) {
return G(x_, s_, y).hold();
}
}
+// option do_not_evalf_params() removed.
+// This is safe: in the code above it only matters if s_ > 0 or s_ < 0,
+// s_ is allowed to be of floating type.
unsigned G3_SERIAL::serial = function::register_new(function_options("G", 3).
evalf_func(G3_evalf).
eval_func(G3_eval).
- do_not_evalf_params().
overloaded(2));
//TODO
// derivative_func(G3_deriv).
{
if (is_a<lst>(m) || is_a<lst>(x)) {
// multiple polylog
- epvector seq;
- seq.push_back(expair(Li(m, x), 0));
- return pseries(rel, seq);
+ epvector seq { expair(Li(m, x), 0) };
+ return pseries(rel, std::move(seq));
}
// classical polylog
// substitute the argument's series expansion
ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
// maybe that was terminating, so add a proper order term
- epvector nseq;
- nseq.push_back(expair(Order(_ex1), order));
- ser += pseries(rel, nseq);
+ epvector nseq { expair(Order(_ex1), order) };
+ ser += pseries(rel, std::move(nseq));
// reexpanding it will collapse the series again
return ser.series(rel, order);
}
if (is_a<lst>(m_)) {
m = ex_to<lst>(m_);
} else {
- m = lst(m_);
+ m = lst{m_};
}
lst x;
if (is_a<lst>(x_)) {
x = ex_to<lst>(x_);
} else {
- x = lst(x_);
+ x = lst{x_};
}
c.s << "\\mathrm{Li}_{";
lst::const_iterator itm = m.begin();
// lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
// see fill_Yn()
-std::vector<std::vector<cln::cl_N> > Yn;
+std::vector<std::vector<cln::cl_N>> Yn;
int ynsize = 0; // number of Yn[]
int ynlength = 100; // initial length of all Yn[i]
return _ex0;
}
if (x == 1) {
- lst m(n+1);
+ lst m{n+1};
for (int i=ex_to<numeric>(p).to_int()-1; i>0; i--) {
m.append(1);
}
// substitute the argument's series expansion
ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
// maybe that was terminating, so add a proper order term
- epvector nseq;
- nseq.push_back(expair(Order(_ex1), order));
- ser += pseries(rel, nseq);
+ epvector nseq { expair(Order(_ex1), order) };
+ ser += pseries(rel, std::move(nseq));
// reexpanding it will collapse the series again
return ser.series(rel, order);
}
// recursivly transforms H to corresponding multiple polylogarithms
struct map_trafo_H_convert_to_Li : public map_function
{
- ex operator()(const ex& e)
+ ex operator()(const ex& e) override
{
if (is_a<add>(e) || is_a<mul>(e)) {
return e.map(*this);
if (is_a<lst>(e.op(0))) {
parameter = ex_to<lst>(e.op(0));
} else {
- parameter = lst(e.op(0));
+ parameter = lst{e.op(0)};
}
ex arg = e.op(1);
// recursivly transforms H to corresponding zetas
struct map_trafo_H_convert_to_zeta : public map_function
{
- ex operator()(const ex& e)
+ ex operator()(const ex& e) override
{
if (is_a<add>(e) || is_a<mul>(e)) {
return e.map(*this);
if (is_a<lst>(e.op(0))) {
parameter = ex_to<lst>(e.op(0));
} else {
- parameter = lst(e.op(0));
+ parameter = lst{e.op(0)};
}
lst m;
// remove trailing zeros from H-parameters
struct map_trafo_H_reduce_trailing_zeros : public map_function
{
- ex operator()(const ex& e)
+ ex operator()(const ex& e) override
{
if (is_a<add>(e) || is_a<mul>(e)) {
return e.map(*this);
if (is_a<lst>(e.op(0))) {
parameter = ex_to<lst>(e.op(0));
} else {
- parameter = lst(e.op(0));
+ parameter = lst{e.op(0)};
}
ex arg = e.op(1);
if (parameter.op(parameter.nops()-1) == 0) {
if (h2nops > 1) {
hlong = ex_to<lst>(h2.op(0));
} else {
- hlong = h2.op(0).op(0);
+ hlong = lst{h2.op(0).op(0)};
}
}
for (std::size_t i=0; i<=hlong.nops(); i++) {
// applies trafo_H_mult recursively on expressions
struct map_trafo_H_mult : public map_function
{
- ex operator()(const ex& e)
+ ex operator()(const ex& e) override
{
if (is_a<add>(e)) {
return e.map(*this);
ex addzeta = convert_H_to_zeta(newparameter);
return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
} else {
- return e * (-H(lst(ex(0)),1/arg).hold());
+ return e * (-H(lst{ex(0)},1/arg).hold());
}
}
newparameter.prepend(1);
return e.subs(h == H(newparameter, h.op(1)).hold());
} else {
- return e * H(lst(ex(1)),1-arg).hold();
+ return e * H(lst{ex(1)},1-arg).hold();
}
}
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(ex(-1)));
- return (e * (addzeta - H(lst(ex(-1)),1/arg).hold())).expand();
+ ex addzeta = convert_H_to_zeta(lst{ex(-1)});
+ return (e * (addzeta - H(lst{ex(-1)},1/arg).hold())).expand();
}
}
newparameter.prepend(-1);
return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
} else {
- return (e * H(lst(ex(-1)),(1-arg)/(1+arg)).hold()).expand();
+ return (e * H(lst{ex(-1)},(1-arg)/(1+arg)).hold()).expand();
}
}
newparameter.prepend(1);
return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
} else {
- return (e * H(lst(ex(1)),(1-arg)/(1+arg)).hold()).expand();
+ return (e * H(lst{ex(1)},(1-arg)/(1+arg)).hold()).expand();
}
}
// do x -> 1-x transformation
struct map_trafo_H_1mx : public map_function
{
- ex operator()(const ex& e)
+ ex operator()(const ex& e) override
{
if (is_a<add>(e) || is_a<mul>(e)) {
return e.map(*this);
// leading one
map_trafo_H_1mx recursion;
map_trafo_H_mult unify;
- ex res = H(lst(ex(1)), arg).hold() * H(newparameter, arg).hold();
+ ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
std::size_t firstzero = 0;
while (parameter.op(firstzero) == 1) {
firstzero++;
// do x -> 1/x transformation
struct map_trafo_H_1overx : public map_function
{
- ex operator()(const ex& e)
+ ex operator()(const ex& e) override
{
if (is_a<add>(e) || is_a<mul>(e)) {
return e.map(*this);
}
if (allthesame) {
map_trafo_H_mult unify;
- return unify((pow(H(lst(ex(-1)),1/arg).hold() - H(lst(ex(0)),1/arg).hold(), parameter.nops())
+ return unify((pow(H(lst{ex(-1)},1/arg).hold() - H(lst{ex(0)},1/arg).hold(), parameter.nops())
/ factorial(parameter.nops())).expand());
}
} else {
}
if (allthesame) {
map_trafo_H_mult unify;
- return unify((pow(H(lst(ex(1)),1/arg).hold() + H(lst(ex(0)),1/arg).hold() + H_polesign, parameter.nops())
+ return unify((pow(H(lst{ex(1)},1/arg).hold() + H(lst{ex(0)},1/arg).hold() + H_polesign, parameter.nops())
/ factorial(parameter.nops())).expand());
}
}
// leading one
map_trafo_H_1overx recursion;
map_trafo_H_mult unify;
- ex res = H(lst(ex(1)), arg).hold() * H(newparameter, arg).hold();
+ ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
std::size_t firstzero = 0;
while (parameter.op(firstzero) == 1) {
firstzero++;
// do x -> (1-x)/(1+x) transformation
struct map_trafo_H_1mxt1px : public map_function
{
- ex operator()(const ex& e)
+ ex operator()(const ex& e) override
{
if (is_a<add>(e) || is_a<mul>(e)) {
return e.map(*this);
}
if (allthesame) {
map_trafo_H_mult unify;
- return unify((pow(-H(lst(ex(1)),(1-arg)/(1+arg)).hold() - H(lst(ex(-1)),(1-arg)/(1+arg)).hold(), parameter.nops())
+ return unify((pow(-H(lst{ex(1)},(1-arg)/(1+arg)).hold() - H(lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
/ factorial(parameter.nops())).expand());
}
} else if (parameter.op(0) == -1) {
}
if (allthesame) {
map_trafo_H_mult unify;
- return unify((pow(log(2) - H(lst(ex(-1)),(1-arg)/(1+arg)).hold(), parameter.nops())
+ return unify((pow(log(2) - H(lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
/ factorial(parameter.nops())).expand());
}
} else {
}
if (allthesame) {
map_trafo_H_mult unify;
- return unify((pow(-log(2) - H(lst(ex(0)),(1-arg)/(1+arg)).hold() + H(lst(ex(-1)),(1-arg)/(1+arg)).hold(), parameter.nops())
+ return unify((pow(-log(2) - H(lst{ex(0)},(1-arg)/(1+arg)).hold() + H(lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
/ factorial(parameter.nops())).expand());
}
}
// leading one
map_trafo_H_1mxt1px recursion;
map_trafo_H_mult unify;
- ex res = H(lst(ex(1)), arg).hold() * H(newparameter, arg).hold();
+ ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
std::size_t firstzero = 0;
while (parameter.op(firstzero) == 1) {
firstzero++;
if (is_a<lst>(m_)) {
m = ex_to<lst>(m_);
} else {
- m = lst(m_);
+ m = lst{m_};
}
if (m.nops() == 0) {
return _ex1;
static ex H_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
{
- epvector seq;
- seq.push_back(expair(H(m, x), 0));
- return pseries(rel, seq);
+ epvector seq { expair(H(m, x), 0) };
+ return pseries(rel, std::move(seq));
}
if (is_a<lst>(m_)) {
m = ex_to<lst>(m_);
} else {
- m = lst(m_);
+ m = lst{m_};
}
ex mb = *m.begin();
if (mb > _ex1) {
if (is_a<lst>(m_)) {
m = ex_to<lst>(m_);
} else {
- m = lst(m_);
+ m = lst{m_};
}
c.s << "\\mathrm{H}_{";
lst::const_iterator itm = m.begin();
if (is_a<lst>(m)) {
return filter2(filter(H(m, x).hold()));
} else {
- return filter2(filter(H(lst(m), x).hold()));
+ return filter2(filter(H(lst{m}, x).hold()));
}
}
int Sm = 0;
int Smp1 = 0;
- std::vector<std::vector<cln::cl_N> > crG(s.size() - 1, std::vector<cln::cl_N>(L2 + 1));
+ std::vector<std::vector<cln::cl_N>> crG(s.size() - 1, std::vector<cln::cl_N>(L2 + 1));
for (int m=0; m < (int)s.size() - 1; m++) {
Sm += s[m];
Smp1 = Sm + s[m+1];
// [Cra] section 4
-static void calc_f(std::vector<std::vector<cln::cl_N> >& f_kj,
+static void calc_f(std::vector<std::vector<cln::cl_N>>& f_kj,
const int maxr, const int L1)
{
cln::cl_N t0, t1, t2, t3, t4;
int i, j, k;
- std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
+ std::vector<std::vector<cln::cl_N>>::iterator it = f_kj.begin();
cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
t0 = cln::exp(-lambda);
// [Cra] (3.1)
static cln::cl_N crandall_Z(const std::vector<int>& s,
- const std::vector<std::vector<cln::cl_N> >& f_kj)
+ const std::vector<std::vector<cln::cl_N>>& f_kj)
{
const int j = s.size();
}
}
- std::vector<std::vector<cln::cl_N> > f_kj(L1);
+ std::vector<std::vector<cln::cl_N>> f_kj(L1);
calc_f(f_kj, maxr, L1);
const cln::cl_N r0factorial = cln::factorial(r[0]-1);
if (is_a<lst>(m_)) {
m = ex_to<lst>(m_);
} else {
- m = lst(m_);
+ m = lst{m_};
}
lst s;
if (is_a<lst>(s_)) {
s = ex_to<lst>(s_);
} else {
- s = lst(s_);
+ s = lst{s_};
}
c.s << "\\zeta(";
lst::const_iterator itm = m.begin();
git://www.ginac.de / ginac.git / blobdiff