- GiNaC-devel - ginac.de

GiNaC and Asymptote
by Vladimir Kisil 25 Apr '06

25 Apr '06

Dear All, With so many exciting activity going around GiNaC already I wish to propose even more. Why do not let GiNaC classes to visualise themselves? There is an Asymptote (http://asymptote.sourceforge.net/) "is a powerful descriptive vector graphics language that provides a natural coordinate-based framework for technical drawing. Labels and equations are typeset with LaTeX, for high-quality PostScript output." Since Asymptote has "C++-like programming syntax"---it looks like a sister package for GiNAC. This can be used to draw anything from Feynman diagram to 2D and 3D graphs of functions. I used Asymptote as standalone application to visualise cycles in GiNaC library (http://arxiv.org/abs/cs.MS/0512073) I think it would be worth to integrate GiNaC with Asymptote. As a start I think to add a method .asy() to the GiNaC::basic which simply put a label at specific position with LaTeX representation of the object. For functions it would be rewritten by a method which try to draw its graph (if this is possible). If the general perception is that this is worthy I may try to produce an initial patch. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv(a)maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/

3 5

Real and imaginary parts.
by Chris Dams 24 Apr '06

24 Apr '06

Dear GiNaCers, I have added the methods .real_part() and .imag_part() to the library for taking real and imaginary parts of expressions. Functions know how to take real and imaginary parts of themselves through two new function options. I did not add the required methods for this to all standard functions but the trigoniometric functions have been done. These new features are documented in the tutorial. Also, I noticed that the capitalization of words in section titles in the tutorial is inconsistent. I made this consistent by making all (except the first) words that are not names start with a lower case letter. Best, Chris

2 3

Bug in ::match.
by Chris Dams 20 Apr '06

20 Apr '06

Dear all, I found out that sometimes match fails even when the pattern matches. This depends on the canonical ordering. Imagine that we have an expression that looks like x^5*y^(-1). and a pattern that looks like ($0)^(-1)*x^($2). Then the first thing that is attempted is to match x^5 with ($0)^(-1). This match will fail, however, the side effect of this match will be that the repl_lst now contains $0==x. This repl_lst will prevent the match of y^(-1) to ($0)^(-1) to succeed. A patch is in CVS. I also patched 1.3. Best, Chris

1 0

function classes
by Jens Vollinga 19 Apr '06

19 Apr '06

Hi, a long time ago I posted on this list an idea about function classes. Finally I some working code to demonstrate this idea (patch for HEAD attached). Apart for some minor issues maybe this code might well go into CVS, I think. But before spoiling the ginac code base, I'd like to ask for opinions or suggestions from the other developers. Things I don't like too much about it are: - function classes have to follow a naming scheme like sin_function, cos_function, etc. Reasons to have such a naming scheme: 1. cmath puts its functions mercilessly into the global namespace. 2. function names and function method names must not be the same. This is important for conjugate for example. 3. To avoid problems with the built-in global ginac functions with numeric arguments. - in order to make input_parser work, there has to be a look-up between function names and some kind of factory. This kind of reproduces the old function system within the new function class system (the patch (ab)uses class_info mechanism to save these information). I am currently not sure whether it is worth introducing this new function system because of the problems mentioned above, the code incompatibility it introduces and maybe some other disadvantages not foreseen by me. Regards, Jens Index: check/exam_indexed.cpp =================================================================== RCS file: /home/cvs/GiNaC/check/exam_indexed.cpp,v retrieving revision 1.26 diff -r1.26 exam_indexed.cpp 151,154c151,189 < DECLARE_FUNCTION_2P(symm_fcn) < REGISTER_FUNCTION(symm_fcn, set_symmetry(sy_symm(0, 1))); < DECLARE_FUNCTION_2P(anti_fcn) < REGISTER_FUNCTION(anti_fcn, set_symmetry(sy_anti(0, 1))); --- > class symm_fcn : public function > { > GINAC_DECLARE_FUNCTION_2P(symm_fcn) > public: > virtual ex eval(int level = 0) const > { > // Canonicalize argument order according to the symmetry properties > exvector v = seq; > int sig = canonicalize(v.begin(), sy_symm(0, 1)); > if (sig != INT_MAX) { > // Something has changed while sorting arguments, more evaluations later > if (sig == 0) > return 0; > return ex(sig) * thiscontainer(v); > } > return this->hold(); > } > }; > GINAC_IMPLEMENT_FUNCTION(symm_fcn) > > class anti_fcn : public function > { > GINAC_DECLARE_FUNCTION_2P(anti_fcn) > public: > virtual ex eval(int level = 0) const > { > // Canonicalize argument order according to the symmetry properties > exvector v = seq; > int sig = canonicalize(v.begin(), sy_anti(0, 1)); > if (sig != INT_MAX) { > // Something has changed while sorting arguments, more evaluations later > if (sig == 0) > return 0; > return ex(sig) * thiscontainer(v); > } > return this->hold(); > } > }; > GINAC_IMPLEMENT_FUNCTION(anti_fcn) Index: check/exam_pseries.cpp =================================================================== RCS file: /home/cvs/GiNaC/check/exam_pseries.cpp,v retrieving revision 1.27 diff -r1.27 exam_pseries.cpp 272,273c272,273 < ex e = Li2(pow(x,2)); < ex d = Li2(4) + (-log(3) + I*Pi*csgn(I-I*pow(x,2))) * (x-2) --- > ex e = Li(2, pow(x,2)); > ex d = Li(2, 4) + (-log(3) + I*Pi*csgn(I-I*pow(x,2))) * (x-2) Index: ginac/Makefile.am =================================================================== RCS file: /home/cvs/GiNaC/ginac/Makefile.am,v retrieving revision 1.36 diff -r1.36 Makefile.am 6,7c6,7 < fail.cpp fderivative.cpp function.cpp idx.cpp indexed.cpp inifcns.cpp \ < inifcns_trans.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp \ --- > fail.cpp function.cpp idx.cpp indexed.cpp inifcns.cpp inifcns_exp.cpp \ > inifcns_polylog.cpp inifcns_trig.cpp \ 9c9 < operators.cpp power.cpp registrar.cpp relational.cpp remember.cpp \ --- > operators.cpp power.cpp registrar.cpp relational.cpp \ 12c12 < input_lexer.h remember.h tostring.h utils.h --- > input_lexer.h tostring.h utils.h 17,18c17,19 < exprseq.h fail.h fderivative.h flags.h function.h hash_map.h idx.h indexed.h \ < inifcns.h integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \ --- > exprseq.h fail.h flags.h function.h hash_map.h idx.h indexed.h inifcns.h \ > inifcns_exp.h inifcns_polylog.h inifcns_trig.h \ > integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \ 23,30c24 < EXTRA_DIST = function.pl input_parser.h version.h.in < < # Files which are generated by perl scripts < $(srcdir)/function.h $(srcdir)/function.cpp: $(srcdir)/function.pl < cd $(srcdir) && perl -w function.pl < < # Force build of headers before compilation < $(srcdir)/add.cpp: $(srcdir)/function.h --- > EXTRA_DIST = input_parser.h version.h.in Index: ginac/fderivative.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/fderivative.cpp,v retrieving revision 1.16 diff -r1.16 fderivative.cpp 1,220d0 < /** @file fderivative.cpp < * < * Implementation of abstract derivatives of functions. */ < < /* < * GiNaC Copyright (C) 1999-2006 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #include <iostream> < < #include "fderivative.h" < #include "operators.h" < #include "archive.h" < #include "utils.h" < < namespace GiNaC { < < GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(fderivative, function, < print_func<print_context>(&fderivative::do_print). < print_func<print_tree>(&fderivative::do_print_tree)) < < ////////// < // default constructor < ////////// < < fderivative::fderivative() < { < tinfo_key = &fderivative::tinfo_static; < } < < ////////// < // other constructors < ////////// < < fderivative::fderivative(unsigned ser, unsigned param, const exvector & args) : function(ser, args) < { < parameter_set.insert(param); < tinfo_key = &fderivative::tinfo_static; < } < < fderivative::fderivative(unsigned ser, const paramset & params, const exvector & args) : function(ser, args), parameter_set(params) < { < tinfo_key = &fderivative::tinfo_static; < } < < fderivative::fderivative(unsigned ser, const paramset & params, std::auto_ptr<exvector> vp) : function(ser, vp), parameter_set(params) < { < tinfo_key = &fderivative::tinfo_static; < } < < ////////// < // archiving < ////////// < < fderivative::fderivative(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) < { < unsigned i = 0; < while (true) { < unsigned u; < if (n.find_unsigned("param", u, i)) < parameter_set.insert(u); < else < break; < ++i; < } < } < < void fderivative::archive(archive_node &n) const < { < inherited::archive(n); < paramset::const_iterator i = parameter_set.begin(), end = parameter_set.end(); < while (i != end) { < n.add_unsigned("param", *i); < ++i; < } < } < < DEFAULT_UNARCHIVE(fderivative) < < ////////// < // functions overriding virtual functions from base classes < ////////// < < void fderivative::print(const print_context & c, unsigned level) const < { < // class function overrides print(), but we don't want that < basic::print(c, level); < } < < void fderivative::do_print(const print_context & c, unsigned level) const < { < c.s << "D["; < paramset::const_iterator i = parameter_set.begin(), end = parameter_set.end(); < --end; < while (i != end) < c.s << *i++ << ","; < c.s << *i << "](" << registered_functions()[serial].name << ")"; < printseq(c, '(', ',', ')', exprseq::precedence(), function::precedence()); < } < < void fderivative::do_print_tree(const print_tree & c, unsigned level) const < { < c.s << std::string(level, ' ') << class_name() << " " < << registered_functions()[serial].name << " @" << this < << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec < << ", nops=" << nops() < << ", params="; < paramset::const_iterator i = parameter_set.begin(), end = parameter_set.end(); < --end; < while (i != end) < c.s << *i++ << ","; < c.s << *i << std::endl; < for (size_t i=0; i<seq.size(); ++i) < seq[i].print(c, level + c.delta_indent); < c.s << std::string(level + c.delta_indent, ' ') << "=====" << std::endl; < } < < ex fderivative::eval(int level) const < { < if (level > 1) { < // first evaluate children, then we will end up here again < return fderivative(serial, parameter_set, evalchildren(level)); < } < < // No parameters specified? Then return the function itself < if (parameter_set.empty()) < return function(serial, seq); < < // If the function in question actually has a derivative, return it < if (registered_functions()[serial].has_derivative() && parameter_set.size() == 1) < return pderivative(*(parameter_set.begin())); < < return this->hold(); < } < < /** Numeric evaluation falls back to evaluation of arguments. < * @see basic::evalf */ < ex fderivative::evalf(int level) const < { < return basic::evalf(level); < } < < /** The series expansion of derivatives falls back to Taylor expansion. < * @see basic::series */ < ex fderivative::series(const relational & r, int order, unsigned options) const < { < return basic::series(r, order, options); < } < < ex fderivative::thiscontainer(const exvector & v) const < { < return fderivative(serial, parameter_set, v); < } < < ex fderivative::thiscontainer(std::auto_ptr<exvector> vp) const < { < return fderivative(serial, parameter_set, vp); < } < < /** Implementation of ex::diff() for derivatives. It applies the chain rule. < * @see ex::diff */ < ex fderivative::derivative(const symbol & s) const < { < ex result; < for (size_t i=0; i<seq.size(); i++) { < ex arg_diff = seq[i].diff(s); < if (!arg_diff.is_zero()) { < paramset ps = parameter_set; < ps.insert(i); < result += arg_diff * fderivative(serial, ps, seq); < } < } < return result; < } < < int fderivative::compare_same_type(const basic & other) const < { < GINAC_ASSERT(is_a<fderivative>(other)); < const fderivative & o = static_cast<const fderivative &>(other); < < if (parameter_set != o.parameter_set) < return parameter_set < o.parameter_set ? -1 : 1; < else < return inherited::compare_same_type(o); < } < < bool fderivative::is_equal_same_type(const basic & other) const < { < GINAC_ASSERT(is_a<fderivative>(other)); < const fderivative & o = static_cast<const fderivative &>(other); < < if (parameter_set != o.parameter_set) < return false; < else < return inherited::is_equal_same_type(o); < } < < bool fderivative::match_same_type(const basic & other) const < { < GINAC_ASSERT(is_a<fderivative>(other)); < const fderivative & o = static_cast<const fderivative &>(other); < < return parameter_set == o.parameter_set; < } < < } // namespace GiNaC Index: ginac/fderivative.h =================================================================== RCS file: /home/cvs/GiNaC/ginac/fderivative.h,v retrieving revision 1.11 diff -r1.11 fderivative.h 1,86d0 < /** @file fderivative.h < * < * Interface to abstract derivatives of functions. */ < < /* < * GiNaC Copyright (C) 1999-2006 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #ifndef __GINAC_FDERIVATIVE_H__ < #define __GINAC_FDERIVATIVE_H__ < < #include <set> < < #include "function.h" < < namespace GiNaC { < < < typedef std::multiset<unsigned> paramset; < < /** This class represents the (abstract) derivative of a symbolic function. < * It is used to represent the derivatives of functions that do not have < * a derivative or series expansion procedure defined. */ < class fderivative : public function < { < GINAC_DECLARE_REGISTERED_CLASS(fderivative, function) < < // other constructors < public: < /** Construct derivative with respect to one parameter. < * < * @param ser Serial number of function < * @param param Number of parameter with respect to which to take the derivative < * @param args Arguments of derivative function */ < fderivative(unsigned ser, unsigned param, const exvector & args); < < /** Construct derivative with respect to multiple parameters. < * < * @param ser Serial number of function < * @param params Set of numbers of parameters with respect to which to take the derivative < * @param args Arguments of derivative function */ < fderivative(unsigned ser, const paramset & params, const exvector & args); < < // internal constructors < fderivative(unsigned ser, const paramset & params, std::auto_ptr<exvector> vp); < < // functions overriding virtual functions from base classes < public: < void print(const print_context & c, unsigned level = 0) const; < ex eval(int level = 0) const; < ex evalf(int level = 0) const; < ex series(const relational & r, int order, unsigned options = 0) const; < ex thiscontainer(const exvector & v) const; < ex thiscontainer(std::auto_ptr<exvector> vp) const; < protected: < ex derivative(const symbol & s) const; < bool is_equal_same_type(const basic & other) const; < bool match_same_type(const basic & other) const; < < // non-virtual functions in this class < protected: < void do_print(const print_context & c, unsigned level) const; < void do_print_tree(const print_tree & c, unsigned level) const; < < // member variables < protected: < paramset parameter_set; /**< Set of parameter numbers with respect to which to take the derivative */ < }; < < } // namespace GiNaC < < #endif // ndef __GINAC_DERIVATIVE_H__ Index: ginac/flags.h =================================================================== RCS file: /home/cvs/GiNaC/ginac/flags.h,v retrieving revision 1.38 diff -r1.38 flags.h 6c6 < * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany --- > * GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany 233a234,236 > // answered by class function > function, > 261,272d263 < /** Strategies how to clean up the function remember cache. < * @see remember_table */ < class remember_strategies { < public: < enum { < delete_never, ///< Let table grow undefinitely < delete_lru, ///< Least recently used < delete_lfu, ///< Least frequently used < delete_cyclic ///< First (oldest) one in list < }; < }; < Index: ginac/function.pl =================================================================== RCS file: /home/cvs/GiNaC/ginac/function.pl,v retrieving revision 1.95 diff -r1.95 function.pl 1,1284d0 < # This perl script automatically generates function.h and function.cpp < < # function.pl options: \$maxargs=${maxargs} < # < # GiNaC Copyright (C) 1999-2006 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 < # the Free Software Foundation; either version 2 of the License, or < # (at your option) any later version. < # < # This program is distributed in the hope that it will be useful, < # but WITHOUT ANY WARRANTY; without even the implied warranty of < # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < # GNU General Public License for more details. < # < # You should have received a copy of the GNU General Public License < # along with this program; if not, write to the Free Software < # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < < $maxargs=14; < < sub generate_seq { < my ($seq_template,$n)=@_; < my ($res,$N); < < $res=''; < for ($N=1; $N<=$n; $N++) { < $res .= eval('"' . $seq_template . '"'); < if ($N!=$n) { < $res .= ', '; < } < } < return $res; < } < < sub generate_from_to { < my ($template,$seq_template1,$seq_template2,$seq_template3,$from,$to)=@_; < my ($res,$N,$SEQ); < < $res=''; < for ($N=$from; $N<=$to; $N++) { < $SEQ1=generate_seq($seq_template1,$N); < $SEQ2=generate_seq($seq_template2,$N); < $SEQ3=generate_seq($seq_template3,$N); < $res .= eval('"' . $template . '"'); < $SEQ1=''; # to avoid main::SEQ1 used only once warning < $SEQ2=''; # same as above < $SEQ3=''; # same as above < } < return $res; < } < < sub generate { < my ($template,$seq_template1,$seq_template2,$seq_template3)=@_; < return generate_from_to($template,$seq_template1,$seq_template2,$seq_template3,1,$maxargs); < } < < $declare_function_macro = generate( < <<'END_OF_DECLARE_FUNCTION_MACRO','typename T${N}','const T${N} & p${N}','GiNaC::ex(p${N})'); < #define DECLARE_FUNCTION_${N}P(NAME) \\ < class NAME##_SERIAL { public: static unsigned serial; }; \\ < const unsigned NAME##_NPARAMS = ${N}; \\ < template<${SEQ1}> const GiNaC::function NAME(${SEQ2}) { \\ < return GiNaC::function(NAME##_SERIAL::serial, ${SEQ3}); \\ < } < < END_OF_DECLARE_FUNCTION_MACRO < < $typedef_eval_funcp=generate( < 'typedef ex (* eval_funcp_${N})(${SEQ1});'."\n", < 'const ex &','',''); < < $typedef_evalf_funcp=generate( < 'typedef ex (* evalf_funcp_${N})(${SEQ1});'."\n", < 'const ex &','',''); < < $typedef_conjugate_funcp=generate( < 'typedef ex (* conjugate_funcp_${N})(${SEQ1});'."\n", < 'const ex &','',''); < < $typedef_derivative_funcp=generate( < 'typedef ex (* derivative_funcp_${N})(${SEQ1}, unsigned);'."\n", < 'const ex &','',''); < < $typedef_power_funcp=generate( < 'typedef ex (* power_funcp_${N})(${SEQ1}, const ex &);'."\n", < 'const ex &','',''); < < $typedef_series_funcp=generate( < 'typedef ex (* series_funcp_${N})(${SEQ1}, const relational &, int, unsigned);'."\n", < 'const ex &','',''); < < $typedef_print_funcp=generate( < 'typedef void (* print_funcp_${N})(${SEQ1}, const print_context &);'."\n", < 'const ex &','',''); < < $eval_func_interface=generate(' function_options & eval_func(eval_funcp_${N} e);'."\n",'','',''); < < $evalf_func_interface=generate(' function_options & evalf_func(evalf_funcp_${N} ef);'."\n",'','',''); < < $conjugate_func_interface=generate(' function_options & conjugate_func(conjugate_funcp_${N} d);'."\n",'','',''); < < $derivative_func_interface=generate(' function_options & derivative_func(derivative_funcp_${N} d);'."\n",'','',''); < < $power_func_interface=generate(' function_options & power_func(power_funcp_${N} d);'."\n",'','',''); < < $series_func_interface=generate(' function_options & series_func(series_funcp_${N} s);'."\n",'','',''); < < $print_func_interface=generate( < <<'END_OF_PRINT_FUNC_INTERFACE','','',''); < template <class Ctx> function_options & print_func(print_funcp_${N} p) < { < test_and_set_nparams(${N}); < set_print_func(Ctx::get_class_info_static().options.get_id(), print_funcp(p)); < return *this; < } < END_OF_PRINT_FUNC_INTERFACE < < $constructors_interface=generate( < ' function(unsigned ser, ${SEQ1});'."\n", < 'const ex & param${N}','',''); < < $constructors_implementation=generate( < <<'END_OF_CONSTRUCTORS_IMPLEMENTATION','const ex & param${N}','param${N}',''); < function::function(unsigned ser, ${SEQ1}) < : exprseq(${SEQ2}), serial(ser) < { < tinfo_key = &function::tinfo_static; < } < END_OF_CONSTRUCTORS_IMPLEMENTATION < < $eval_switch_statement=generate( < <<'END_OF_EVAL_SWITCH_STATEMENT','seq[${N}-1]','',''); < case ${N}: < eval_result = ((eval_funcp_${N})(opt.eval_f))(${SEQ1}); < break; < END_OF_EVAL_SWITCH_STATEMENT < < $evalf_switch_statement=generate( < <<'END_OF_EVALF_SWITCH_STATEMENT','eseq[${N}-1]','',''); < case ${N}: < return ((evalf_funcp_${N})(opt.evalf_f))(${SEQ1}); < END_OF_EVALF_SWITCH_STATEMENT < < $conjugate_switch_statement=generate( < <<'END_OF_DIFF_SWITCH_STATEMENT','seq[${N}-1]','',''); < case ${N}: < return ((conjugate_funcp_${N})(opt.conjugate_f))(${SEQ1}); < END_OF_DIFF_SWITCH_STATEMENT < < $diff_switch_statement=generate( < <<'END_OF_DIFF_SWITCH_STATEMENT','seq[${N}-1]','',''); < case ${N}: < return ((derivative_funcp_${N})(opt.derivative_f))(${SEQ1},diff_param); < END_OF_DIFF_SWITCH_STATEMENT < < $power_switch_statement=generate( < <<'END_OF_POWER_SWITCH_STATEMENT','seq[${N}-1]','',''); < case ${N}: < return ((power_funcp_${N})(opt.power_f))(${SEQ1},power_param); < END_OF_POWER_SWITCH_STATEMENT < < $series_switch_statement=generate( < <<'END_OF_SERIES_SWITCH_STATEMENT','seq[${N}-1]','',''); < case ${N}: < try { < res = ((series_funcp_${N})(opt.series_f))(${SEQ1},r,order,options); < } catch (do_taylor) { < res = basic::series(r, order, options); < } < return res; < END_OF_SERIES_SWITCH_STATEMENT < < $print_switch_statement=generate( < <<'END_OF_PRINT_SWITCH_STATEMENT','seq[${N}-1]','',''); < case ${N}: < ((print_funcp_${N})(pdt[id]))(${SEQ1}, c); < break; < END_OF_PRINT_SWITCH_STATEMENT < < $eval_func_implementation=generate( < <<'END_OF_EVAL_FUNC_IMPLEMENTATION','','',''); < function_options & function_options::eval_func(eval_funcp_${N} e) < { < test_and_set_nparams(${N}); < eval_f = eval_funcp(e); < return *this; < } < END_OF_EVAL_FUNC_IMPLEMENTATION < < $evalf_func_implementation=generate( < <<'END_OF_EVALF_FUNC_IMPLEMENTATION','','',''); < function_options & function_options::evalf_func(evalf_funcp_${N} ef) < { < test_and_set_nparams(${N}); < evalf_f = evalf_funcp(ef); < return *this; < } < END_OF_EVALF_FUNC_IMPLEMENTATION < < $conjugate_func_implementation=generate( < <<'END_OF_CONJUGATE_FUNC_IMPLEMENTATION','','',''); < function_options & function_options::conjugate_func(conjugate_funcp_${N} c) < { < test_and_set_nparams(${N}); < conjugate_f = conjugate_funcp(c); < return *this; < } < END_OF_CONJUGATE_FUNC_IMPLEMENTATION < < $derivative_func_implementation=generate( < <<'END_OF_DERIVATIVE_FUNC_IMPLEMENTATION','','',''); < function_options & function_options::derivative_func(derivative_funcp_${N} d) < { < test_and_set_nparams(${N}); < derivative_f = derivative_funcp(d); < return *this; < } < END_OF_DERIVATIVE_FUNC_IMPLEMENTATION < < $power_func_implementation=generate( < <<'END_OF_POWER_FUNC_IMPLEMENTATION','','',''); < function_options & function_options::power_func(power_funcp_${N} d) < { < test_and_set_nparams(${N}); < power_f = power_funcp(d); < return *this; < } < END_OF_POWER_FUNC_IMPLEMENTATION < < $series_func_implementation=generate( < <<'END_OF_SERIES_FUNC_IMPLEMENTATION','','',''); < function_options & function_options::series_func(series_funcp_${N} s) < { < test_and_set_nparams(${N}); < series_f = series_funcp(s); < return *this; < } < END_OF_SERIES_FUNC_IMPLEMENTATION < < $interface=<<END_OF_INTERFACE; < /** \@file function.h < * < * Interface to class of symbolic functions. */ < < /* < * This file was generated automatically by function.pl. < * Please do not modify it directly, edit the perl script instead! < * function.pl options: \$maxargs=${maxargs} < * < * GiNaC Copyright (C) 1999-2005 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #ifndef __GINAC_FUNCTION_H__ < #define __GINAC_FUNCTION_H__ < < #include <string> < #include <vector> < < // CINT needs <algorithm> to work properly with <vector> < #include <algorithm> < < #include "exprseq.h" < < // the following lines have been generated for max. ${maxargs} parameters < $declare_function_macro < // end of generated lines < < #define REGISTER_FUNCTION(NAME,OPT) \\ < unsigned NAME##_SERIAL::serial = \\ < GiNaC::function::register_new(GiNaC::function_options(#NAME, NAME##_NPARAMS).OPT); < < namespace GiNaC { < < class function; < class symmetry; < < typedef ex (* eval_funcp)(); < typedef ex (* evalf_funcp)(); < typedef ex (* conjugate_funcp)(); < typedef ex (* derivative_funcp)(); < typedef ex (* power_funcp)(); < typedef ex (* series_funcp)(); < typedef void (* print_funcp)(); < < // the following lines have been generated for max. ${maxargs} parameters < $typedef_eval_funcp < $typedef_evalf_funcp < $typedef_conjugate_funcp < $typedef_derivative_funcp < $typedef_power_funcp < $typedef_series_funcp < $typedef_print_funcp < // end of generated lines < < // Alternatively, an exvector may be passed into the static function, instead < // of individual ex objects. Then, the number of arguments is not limited. < typedef ex (* eval_funcp_exvector)(const exvector &); < typedef ex (* evalf_funcp_exvector)(const exvector &); < typedef ex (* conjugate_funcp_exvector)(const exvector &); < typedef ex (* derivative_funcp_exvector)(const exvector &, unsigned); < typedef ex (* power_funcp_exvector)(const exvector &, const ex &); < typedef ex (* series_funcp_exvector)(const exvector &, const relational &, int, unsigned); < typedef void (* print_funcp_exvector)(const exvector &, const print_context &); < < < class function_options < { < friend class function; < friend class fderivative; < public: < function_options(); < function_options(std::string const & n, std::string const & tn=std::string()); < function_options(std::string const & n, unsigned np); < ~function_options(); < void initialize(); < < function_options & dummy() { return *this; } < function_options & set_name(std::string const & n, std::string const & tn=std::string()); < function_options & latex_name(std::string const & tn); < // the following lines have been generated for max. ${maxargs} parameters < $eval_func_interface < $evalf_func_interface < $conjugate_func_interface < $derivative_func_interface < $power_func_interface < $series_func_interface < $print_func_interface < // end of generated lines < function_options & eval_func(eval_funcp_exvector e); < function_options & evalf_func(evalf_funcp_exvector ef); < function_options & conjugate_func(conjugate_funcp_exvector d); < function_options & derivative_func(derivative_funcp_exvector d); < function_options & power_func(power_funcp_exvector d); < function_options & series_func(series_funcp_exvector s); < < template <class Ctx> function_options & print_func(print_funcp_exvector p) < { < print_use_exvector_args = true; < set_print_func(Ctx::get_class_info_static().options.get_id(), print_funcp(p)); < return *this; < } < < function_options & set_return_type(unsigned rt, tinfo_t rtt=NULL); < function_options & do_not_evalf_params(); < function_options & remember(unsigned size, unsigned assoc_size=0, < unsigned strategy=remember_strategies::delete_never); < function_options & overloaded(unsigned o); < function_options & set_symmetry(const symmetry & s); < < std::string get_name() const { return name; } < unsigned get_nparams() const { return nparams; } < < protected: < bool has_derivative() const { return derivative_f != NULL; } < bool has_power() const { return power_f != NULL; } < void test_and_set_nparams(unsigned n); < void set_print_func(unsigned id, print_funcp f); < < std::string name; < std::string TeX_name; < < unsigned nparams; < < eval_funcp eval_f; < evalf_funcp evalf_f; < conjugate_funcp conjugate_f; < derivative_funcp derivative_f; < power_funcp power_f; < series_funcp series_f; < std::vector<print_funcp> print_dispatch_table; < < bool evalf_params_first; < < bool use_return_type; < unsigned return_type; < tinfo_t return_type_tinfo; < < bool use_remember; < unsigned remember_size; < unsigned remember_assoc_size; < unsigned remember_strategy; < < bool eval_use_exvector_args; < bool evalf_use_exvector_args; < bool conjugate_use_exvector_args; < bool derivative_use_exvector_args; < bool power_use_exvector_args; < bool series_use_exvector_args; < bool print_use_exvector_args; < < unsigned functions_with_same_name; < < ex symtree; < }; < < < /** Exception class thrown by classes which provide their own series expansion < * to signal that ordinary Taylor expansion is safe. */ < class do_taylor {}; < < < /** The class function is used to implement builtin functions like sin, cos... < and user defined functions */ < class function : public exprseq < { < GINAC_DECLARE_REGISTERED_CLASS(function, exprseq) < < // CINT has a linking problem < #ifndef __MAKECINT__ < friend void ginsh_get_ginac_functions(); < #endif // def __MAKECINT__ < < friend class remember_table_entry; < // friend class remember_table_list; < // friend class remember_table; < < // member functions < < // other constructors < public: < function(unsigned ser); < // the following lines have been generated for max. ${maxargs} parameters < $constructors_interface < // end of generated lines < function(unsigned ser, const exprseq & es); < function(unsigned ser, const exvector & v, bool discardable = false); < function(unsigned ser, std::auto_ptr<exvector> vp); < < // functions overriding virtual functions from base classes < public: < void print(const print_context & c, unsigned level = 0) const; < unsigned precedence() const {return 70;} < ex expand(unsigned options=0) const; < ex eval(int level=0) const; < ex evalf(int level=0) const; < unsigned calchash() const; < ex series(const relational & r, int order, unsigned options = 0) const; < ex thiscontainer(const exvector & v) const; < ex thiscontainer(std::auto_ptr<exvector> vp) const; < ex conjugate() const; < protected: < ex derivative(const symbol & s) const; < bool is_equal_same_type(const basic & other) const; < bool match_same_type(const basic & other) const; < unsigned return_type() const; < tinfo_t return_type_tinfo() const; < < // new virtual functions which can be overridden by derived classes < // none < < // non-virtual functions in this class < protected: < ex pderivative(unsigned diff_param) const; // partial differentiation < static std::vector<function_options> & registered_functions(); < bool lookup_remember_table(ex & result) const; < void store_remember_table(ex const & result) const; < public: < ex power(const ex & exp) const; < static unsigned register_new(function_options const & opt); < static unsigned current_serial; < static unsigned find_function(const std::string &name, unsigned nparams); < unsigned get_serial() const {return serial;} < std::string get_name() const; < < // member variables < < protected: < unsigned serial; < }; < < // utility functions/macros < < template <typename T> < inline bool is_the_function(const ex & x) < { < return is_exactly_a<function>(x) < && ex_to<function>(x).get_serial() == T::serial; < } < < // Check whether OBJ is the specified symbolic function. < #define is_ex_the_function(OBJ, FUNCNAME) (GiNaC::is_the_function<FUNCNAME##_SERIAL>(OBJ)) < < } // namespace GiNaC < < #endif // ndef __GINAC_FUNCTION_H__ < < END_OF_INTERFACE < < $implementation=<<END_OF_IMPLEMENTATION; < /** \@file function.cpp < * < * Implementation of class of symbolic functions. */ < < /* < * This file was generated automatically by function.pl. < * Please do not modify it directly, edit the perl script instead! < * function.pl options: \$maxargs=${maxargs} < * < * GiNaC Copyright (C) 1999-2005 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #include <iostream> < #include <string> < #include <stdexcept> < #include <list> < < #include "function.h" < #include "operators.h" < #include "fderivative.h" < #include "ex.h" < #include "lst.h" < #include "symmetry.h" < #include "print.h" < #include "power.h" < #include "archive.h" < #include "inifcns.h" < #include "tostring.h" < #include "utils.h" < #include "remember.h" < < namespace GiNaC { < < ////////// < // helper class function_options < ////////// < < function_options::function_options() < { < initialize(); < } < < function_options::function_options(std::string const & n, std::string const & tn) < { < initialize(); < set_name(n, tn); < } < < function_options::function_options(std::string const & n, unsigned np) < { < initialize(); < set_name(n, std::string()); < nparams = np; < } < < function_options::~function_options() < { < // nothing to clean up at the moment < } < < void function_options::initialize() < { < set_name("unnamed_function", "\\\\mbox{unnamed}"); < nparams = 0; < eval_f = evalf_f = conjugate_f = derivative_f = power_f = series_f = 0; < evalf_params_first = true; < use_return_type = false; < eval_use_exvector_args = false; < evalf_use_exvector_args = false; < conjugate_use_exvector_args = false; < derivative_use_exvector_args = false; < power_use_exvector_args = false; < series_use_exvector_args = false; < print_use_exvector_args = false; < use_remember = false; < functions_with_same_name = 1; < symtree = 0; < } < < function_options & function_options::set_name(std::string const & n, < std::string const & tn) < { < name = n; < if (tn==std::string()) < TeX_name = "\\\\mbox{"+name+"}"; < else < TeX_name = tn; < return *this; < } < < function_options & function_options::latex_name(std::string const & tn) < { < TeX_name = tn; < return *this; < } < < // the following lines have been generated for max. ${maxargs} parameters < $eval_func_implementation < $evalf_func_implementation < $conjugate_func_implementation < $derivative_func_implementation < $power_func_implementation < $series_func_implementation < // end of generated lines < < function_options& function_options::eval_func(eval_funcp_exvector e) < { < eval_use_exvector_args = true; < eval_f = eval_funcp(e); < return *this; < } < function_options& function_options::evalf_func(evalf_funcp_exvector ef) < { < evalf_use_exvector_args = true; < evalf_f = evalf_funcp(ef); < return *this; < } < function_options& function_options::conjugate_func(conjugate_funcp_exvector c) < { < conjugate_use_exvector_args = true; < conjugate_f = conjugate_funcp(c); < return *this; < } < function_options& function_options::derivative_func(derivative_funcp_exvector d) < { < derivative_use_exvector_args = true; < derivative_f = derivative_funcp(d); < return *this; < } < function_options& function_options::power_func(power_funcp_exvector d) < { < power_use_exvector_args = true; < power_f = power_funcp(d); < return *this; < } < function_options& function_options::series_func(series_funcp_exvector s) < { < series_use_exvector_args = true; < series_f = series_funcp(s); < return *this; < } < < function_options & function_options::set_return_type(unsigned rt, tinfo_t rtt) < { < use_return_type = true; < return_type = rt; < return_type_tinfo = rtt; < return *this; < } < < function_options & function_options::do_not_evalf_params() < { < evalf_params_first = false; < return *this; < } < < function_options & function_options::remember(unsigned size, < unsigned assoc_size, < unsigned strategy) < { < use_remember = true; < remember_size = size; < remember_assoc_size = assoc_size; < remember_strategy = strategy; < return *this; < } < < function_options & function_options::overloaded(unsigned o) < { < functions_with_same_name = o; < return *this; < } < < function_options & function_options::set_symmetry(const symmetry & s) < { < symtree = s; < return *this; < } < < void function_options::test_and_set_nparams(unsigned n) < { < if (nparams==0) { < nparams = n; < } else if (nparams!=n) { < // we do not throw an exception here because this code is < // usually executed before main(), so the exception could not < // be caught anyhow < std::cerr << "WARNING: " << name << "(): number of parameters (" < << n << ") differs from number set before (" < << nparams << ")" << std::endl; < } < } < < void function_options::set_print_func(unsigned id, print_funcp f) < { < if (id >= print_dispatch_table.size()) < print_dispatch_table.resize(id + 1); < print_dispatch_table[id] = f; < } < < /** This can be used as a hook for external applications. */ < unsigned function::current_serial = 0; < < < GINAC_IMPLEMENT_REGISTERED_CLASS(function, exprseq) < < ////////// < // default constructor < ////////// < < // public < < function::function() : serial(0) < { < tinfo_key = &function::tinfo_static; < } < < ////////// < // other constructors < ////////// < < // public < < function::function(unsigned ser) : serial(ser) < { < tinfo_key = &function::tinfo_static; < } < < // the following lines have been generated for max. ${maxargs} parameters < $constructors_implementation < // end of generated lines < < function::function(unsigned ser, const exprseq & es) : exprseq(es), serial(ser) < { < tinfo_key = &function::tinfo_static; < < // Force re-evaluation even if the exprseq was already evaluated < // (the exprseq copy constructor copies the flags) < clearflag(status_flags::evaluated); < } < < function::function(unsigned ser, const exvector & v, bool discardable) < : exprseq(v,discardable), serial(ser) < { < tinfo_key = &function::tinfo_static; < } < < function::function(unsigned ser, std::auto_ptr<exvector> vp) < : exprseq(vp), serial(ser) < { < tinfo_key = &function::tinfo_static; < } < < ////////// < // archiving < ////////// < < /** Construct object from archive_node. */ < function::function(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) < { < // Find serial number by function name < std::string s; < if (n.find_string("name", s)) { < unsigned int ser = 0; < std::vector<function_options>::const_iterator i = registered_functions().begin(), iend = registered_functions().end(); < while (i != iend) { < if (s == i->name) { < serial = ser; < return; < } < ++i; ++ser; < } < throw (std::runtime_error("unknown function '" + s + "' in archive")); < } else < throw (std::runtime_error("unnamed function in archive")); < } < < /** Unarchive the object. */ < ex function::unarchive(const archive_node &n, lst &sym_lst) < { < return (new function(n, sym_lst))->setflag(status_flags::dynallocated); < } < < /** Archive the object. */ < void function::archive(archive_node &n) const < { < inherited::archive(n); < GINAC_ASSERT(serial < registered_functions().size()); < n.add_string("name", registered_functions()[serial].name); < } < < ////////// < // functions overriding virtual functions from base classes < ////////// < < // public < < void function::print(const print_context & c, unsigned level) const < { < GINAC_ASSERT(serial<registered_functions().size()); < const function_options &opt = registered_functions()[serial]; < const std::vector<print_funcp> &pdt = opt.print_dispatch_table; < < // Dynamically dispatch on print_context type < const print_context_class_info *pc_info = &c.get_class_info(); < < next_context: < unsigned id = pc_info->options.get_id(); < if (id >= pdt.size() || pdt[id] == NULL) { < < // Method not found, try parent print_context class < const print_context_class_info *parent_pc_info = pc_info->get_parent(); < if (parent_pc_info) { < pc_info = parent_pc_info; < goto next_context; < } < < // Method still not found, use default output < if (is_a<print_tree>(c)) { < < c.s << std::string(level, ' ') << class_name() << " " < << opt.name << " @" << this < << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec < << ", nops=" << nops() < << std::endl; < unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent; < for (size_t i=0; i<seq.size(); ++i) < seq[i].print(c, level + delta_indent); < c.s << std::string(level + delta_indent, ' ') << "=====" << std::endl; < < } else if (is_a<print_csrc>(c)) { < < // Print function name in lowercase < std::string lname = opt.name; < size_t num = lname.size(); < for (size_t i=0; i<num; i++) < lname[i] = tolower(lname[i]); < c.s << lname; < printseq(c, '(', ',', ')', exprseq::precedence(), function::precedence()); < < } else if (is_a<print_latex>(c)) { < c.s << opt.TeX_name; < printseq(c, '(', ',', ')', exprseq::precedence(), function::precedence()); < } else { < c.s << opt.name; < printseq(c, '(', ',', ')', exprseq::precedence(), function::precedence()); < } < < } else { < < // Method found, call it < current_serial = serial; < if (opt.print_use_exvector_args) < ((print_funcp_exvector)pdt[id])(seq, c); < else switch (opt.nparams) { < // the following lines have been generated for max. ${maxargs} parameters < ${print_switch_statement} < // end of generated lines < default: < throw(std::logic_error("function::print(): invalid nparams")); < } < } < } < < ex function::expand(unsigned options) const < { < // Only expand arguments when asked to do so < if (options & expand_options::expand_function_args) < return inherited::expand(options); < else < return (options == 0) ? setflag(status_flags::expanded) : *this; < } < < ex function::eval(int level) const < { < if (level>1) { < // first evaluate children, then we will end up here again < return function(serial,evalchildren(level)); < } < < GINAC_ASSERT(serial<registered_functions().size()); < const function_options &opt = registered_functions()[serial]; < < // Canonicalize argument order according to the symmetry properties < if (seq.size() > 1 && !(opt.symtree.is_zero())) { < exvector v = seq; < GINAC_ASSERT(is_a<symmetry>(opt.symtree)); < int sig = canonicalize(v.begin(), ex_to<symmetry>(opt.symtree)); < if (sig != INT_MAX) { < // Something has changed while sorting arguments, more evaluations later < if (sig == 0) < return _ex0; < return ex(sig) * thiscontainer(v); < } < } < < if (opt.eval_f==0) { < return this->hold(); < } < < bool use_remember = opt.use_remember; < ex eval_result; < if (use_remember && lookup_remember_table(eval_result)) { < return eval_result; < } < current_serial = serial; < if (opt.eval_use_exvector_args) < eval_result = ((eval_funcp_exvector)(opt.eval_f))(seq); < else < switch (opt.nparams) { < // the following lines have been generated for max. ${maxargs} parameters < ${eval_switch_statement} < // end of generated lines < default: < throw(std::logic_error("function::eval(): invalid nparams")); < } < if (use_remember) { < store_remember_table(eval_result); < } < return eval_result; < } < < ex function::evalf(int level) const < { < GINAC_ASSERT(serial<registered_functions().size()); < const function_options &opt = registered_functions()[serial]; < < // Evaluate children first < exvector eseq; < if (level == 1 || !(opt.evalf_params_first)) < eseq = seq; < else if (level == -max_recursion_level) < throw(std::runtime_error("max recursion level reached")); < else { < eseq.reserve(seq.size()); < --level; < exvector::const_iterator it = seq.begin(), itend = seq.end(); < while (it != itend) { < eseq.push_back(it->evalf(level)); < ++it; < } < } < < if (opt.evalf_f==0) { < return function(serial,eseq).hold(); < } < current_serial = serial; < if (opt.evalf_use_exvector_args) < return ((evalf_funcp_exvector)(opt.evalf_f))(seq); < switch (opt.nparams) { < // the following lines have been generated for max. ${maxargs} parameters < ${evalf_switch_statement} < // end of generated lines < } < throw(std::logic_error("function::evalf(): invalid nparams")); < } < < unsigned function::calchash() const < { < unsigned v = golden_ratio_hash(golden_ratio_hash((p_int)tinfo()) ^ serial); < for (size_t i=0; i<nops(); i++) { < v = rotate_left(v); < v ^= this->op(i).gethash(); < } < < if (flags & status_flags::evaluated) { < setflag(status_flags::hash_calculated); < hashvalue = v; < } < return v; < } < < ex function::thiscontainer(const exvector & v) const < { < return function(serial, v); < } < < ex function::thiscontainer(std::auto_ptr<exvector> vp) const < { < return function(serial, vp); < } < < /** Implementation of ex::series for functions. < * \@see ex::series */ < ex function::series(const relational & r, int order, unsigned options) const < { < GINAC_ASSERT(serial<registered_functions().size()); < const function_options &opt = registered_functions()[serial]; < < if (opt.series_f==0) { < return basic::series(r, order); < } < ex res; < current_serial = serial; < if (opt.series_use_exvector_args) { < try { < res = ((series_funcp_exvector)(opt.series_f))(seq, r, order, options); < } catch (do_taylor) { < res = basic::series(r, order, options); < } < return res; < } < switch (opt.nparams) { < // the following lines have been generated for max. ${maxargs} parameters < ${series_switch_statement} < // end of generated lines < } < throw(std::logic_error("function::series(): invalid nparams")); < } < < /** Implementation of ex::conjugate for functions. */ < ex function::conjugate() const < { < GINAC_ASSERT(serial<registered_functions().size()); < const function_options & opt = registered_functions()[serial]; < < if (opt.conjugate_f==0) { < return exprseq::conjugate(); < } < < if (opt.conjugate_use_exvector_args) { < return ((conjugate_funcp_exvector)(opt.conjugate_f))(seq); < } < < switch (opt.nparams) { < // the following lines have been generated for max. ${maxargs} parameters < ${conjugate_switch_statement} < // end of generated lines < } < throw(std::logic_error("function::conjugate(): invalid nparams")); < } < < // protected < < /** Implementation of ex::diff() for functions. It applies the chain rule, < * except for the Order term function. < * \@see ex::diff */ < ex function::derivative(const symbol & s) const < { < ex result; < < if (serial == Order_SERIAL::serial) { < // Order Term function only differentiates the argument < return Order(seq[0].diff(s)); < } else { < // Chain rule < ex arg_diff; < size_t num = seq.size(); < for (size_t i=0; i<num; i++) { < arg_diff = seq[i].diff(s); < // We apply the chain rule only when it makes sense. This is not < // just for performance reasons but also to allow functions to < // throw when differentiated with respect to one of its arguments < // without running into trouble with our automatic full < // differentiation: < if (!arg_diff.is_zero()) < result += pderivative(i)*arg_diff; < } < } < return result; < } < < int function::compare_same_type(const basic & other) const < { < GINAC_ASSERT(is_a<function>(other)); < const function & o = static_cast<const function &>(other); < < if (serial != o.serial) < return serial < o.serial ? -1 : 1; < else < return exprseq::compare_same_type(o); < } < < bool function::is_equal_same_type(const basic & other) const < { < GINAC_ASSERT(is_a<function>(other)); < const function & o = static_cast<const function &>(other); < < if (serial != o.serial) < return false; < else < return exprseq::is_equal_same_type(o); < } < < bool function::match_same_type(const basic & other) const < { < GINAC_ASSERT(is_a<function>(other)); < const function & o = static_cast<const function &>(other); < < return serial == o.serial; < } < < unsigned function::return_type() const < { < GINAC_ASSERT(serial<registered_functions().size()); < const function_options &opt = registered_functions()[serial]; < < if (opt.use_return_type) { < // Return type was explicitly specified < return opt.return_type; < } else { < // Default behavior is to use the return type of the first < // argument. Thus, exp() of a matrix behaves like a matrix, etc. < if (seq.empty()) < return return_types::commutative; < else < return seq.begin()->return_type(); < } < } < < tinfo_t function::return_type_tinfo() const < { < GINAC_ASSERT(serial<registered_functions().size()); < const function_options &opt = registered_functions()[serial]; < < if (opt.use_return_type) { < // Return type was explicitly specified < return opt.return_type_tinfo; < } else { < // Default behavior is to use the return type of the first < // argument. Thus, exp() of a matrix behaves like a matrix, etc. < if (seq.empty()) < return this; < else < return seq.begin()->return_type_tinfo(); < } < } < < ////////// < // new virtual functions which can be overridden by derived classes < ////////// < < // none < < ////////// < // non-virtual functions in this class < ////////// < < // protected < < ex function::pderivative(unsigned diff_param) const // partial differentiation < { < GINAC_ASSERT(serial<registered_functions().size()); < const function_options &opt = registered_functions()[serial]; < < // No derivative defined? Then return abstract derivative object < if (opt.derivative_f == NULL) < return fderivative(serial, diff_param, seq); < < current_serial = serial; < if (opt.derivative_use_exvector_args) < return ((derivative_funcp_exvector)(opt.derivative_f))(seq, diff_param); < switch (opt.nparams) { < // the following lines have been generated for max. ${maxargs} parameters < ${diff_switch_statement} < // end of generated lines < } < throw(std::logic_error("function::pderivative(): no diff function defined")); < } < < ex function::power(const ex & power_param) const // power of function < { < GINAC_ASSERT(serial<registered_functions().size()); < const function_options &opt = registered_functions()[serial]; < < // No derivative defined? Then return abstract derivative object < if (opt.power_f == NULL) < return (new power::power(*this, power_param))->setflag(status_flags::dynallocated | < status_flags::evaluated); < < current_serial = serial; < if (opt.power_use_exvector_args) < return ((power_funcp_exvector)(opt.power_f))(seq, power_param); < switch (opt.nparams) { < // the following lines have been generated for max. ${maxargs} parameters < ${power_switch_statement} < // end of generated lines < } < throw(std::logic_error("function::power(): no power function defined")); < } < < std::vector<function_options> & function::registered_functions() < { < static std::vector<function_options> * rf = new std::vector<function_options>; < return *rf; < } < < bool function::lookup_remember_table(ex & result) const < { < return remember_table::remember_tables()[this->serial].lookup_entry(*this,result); < } < < void function::store_remember_table(ex const & result) const < { < remember_table::remember_tables()[this->serial].add_entry(*this,result); < } < < // public < < unsigned function::register_new(function_options const & opt) < { < size_t same_name = 0; < for (size_t i=0; i<registered_functions().size(); ++i) { < if (registered_functions()[i].name==opt.name) { < ++same_name; < } < } < if (same_name>=opt.functions_with_same_name) { < // we do not throw an exception here because this code is < // usually executed before main(), so the exception could not < // caught anyhow < std::cerr << "WARNING: function name " << opt.name < << " already in use!" << std::endl; < } < registered_functions().push_back(opt); < if (opt.use_remember) { < remember_table::remember_tables(). < push_back(remember_table(opt.remember_size, < opt.remember_assoc_size, < opt.remember_strategy)); < } else { < remember_table::remember_tables().push_back(remember_table()); < } < return registered_functions().size()-1; < } < < /** Find serial number of function by name and number of parameters. < * Throws exception if function was not found. */ < unsigned function::find_function(const std::string &name, unsigned nparams) < { < std::vector<function_options>::const_iterator i = function::registered_functions().begin(), end = function::registered_functions().end(); < unsigned serial = 0; < while (i != end) { < if (i->get_name() == name && i->get_nparams() == nparams) < return serial; < ++i; < ++serial; < } < throw (std::runtime_error("no function '" + name + "' with " + ToString(nparams) + " parameters defined")); < } < < /** Return the print name of the function. */ < std::string function::get_name() const < { < GINAC_ASSERT(serial<registered_functions().size()); < return registered_functions()[serial].name; < } < < } // namespace GiNaC < < END_OF_IMPLEMENTATION < < print "Creating interface file function.h..."; < open OUT,">function.h" or die "cannot open function.h"; < print OUT $interface; < close OUT; < print "ok.\n"; < < print "Creating implementation file function.cpp..."; < open OUT,">function.cpp" or die "cannot open function.cpp"; < print OUT $implementation; < close OUT; < print "ok.\n"; < < print "done.\n"; Index: ginac/ginac.h =================================================================== RCS file: /home/cvs/GiNaC/ginac/ginac.h,v retrieving revision 1.26 diff -r1.26 ginac.h 6c6 < * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany --- > * GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany 59c59,61 < #include "fderivative.h" --- > #include "inifcns_trig.h" > #include "inifcns_exp.h" > #include "inifcns_polylog.h" Index: ginac/inifcns.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/inifcns.cpp,v retrieving revision 1.88 diff -r1.88 inifcns.cpp 6c6 < * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany --- > * GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany 23,25d22 < #include <vector> < #include <stdexcept> < 26a24 > 39a38,40 > #include <stdexcept> > #include <vector> > 43c44 < // complex conjugate --- > // absolute value 46,57c47,49 < static ex conjugate_evalf(const ex & arg) < { < if (is_exactly_a<numeric>(arg)) { < return ex_to<numeric>(arg).conjugate(); < } < return conjugate_function(arg).hold(); < } < < static ex conjugate_eval(const ex & arg) < { < return arg.conjugate(); < } --- > GINAC_IMPLEMENT_FUNCTION_OPT(abs_function, > print_func<print_csrc_float>(&abs_function::do_print_csrc_float). > print_func<print_latex>(&abs_function::do_print_latex)) 59c51 < static void conjugate_print_latex(const ex & arg, const print_context & c) --- > ex abs_function::conjugate() const 61c53 < c.s << "\\bar{"; arg.print(c); c.s << "}"; --- > return *this; 64,87c56 < static ex conjugate_conjugate(const ex & arg) < { < return arg; < } < < REGISTER_FUNCTION(conjugate_function, eval_func(conjugate_eval). < evalf_func(conjugate_evalf). < print_func<print_latex>(conjugate_print_latex). < conjugate_func(conjugate_conjugate). < set_name("conjugate","conjugate")); < < ////////// < // absolute value < ////////// < < static ex abs_evalf(const ex & arg) < { < if (is_exactly_a<numeric>(arg)) < return abs(ex_to<numeric>(arg)); < < return abs(arg).hold(); < } < < static ex abs_eval(const ex & arg) --- > ex abs_function::eval(int level) const 88a58 > const ex& arg = seq[0]; 92c62 < return abs(arg).hold(); --- > return this->hold(); 95c65 < static void abs_print_latex(const ex & arg, const print_context & c) --- > ex abs_function::evalf(int level) const 97c67,72 < c.s << "{|"; arg.print(c); c.s << "|}"; --- > const ex& arg = seq[0]; > if (is_exactly_a<numeric>(arg)) { > return abs(ex_to<numeric>(arg)); > } > > return this->hold(); 100c75 < static void abs_print_csrc_float(const ex & arg, const print_context & c) --- > ex abs_function::power_law(const ex& exp) const 102c77,81 < c.s << "fabs("; arg.print(c); c.s << ")"; --- > const ex& arg = seq[0]; > if (arg.is_equal(arg.conjugate()) && is_a<numeric>(exp) && ex_to<numeric>(exp).is_even()) > return power(arg, exp); > else > return power(abs(arg), exp).hold(); 105c84 < static ex abs_conjugate(const ex & arg) --- > void abs_function::do_print_csrc_float(const print_context& c, unsigned level) const 107c86 < return abs(arg); --- > c.s << "fabs("; seq[0].print(c); c.s << ")"; 110c89 < static ex abs_power(const ex & arg, const ex & exp) --- > void abs_function::do_print_latex(const print_context& c, unsigned level) const 112,115c91 < if (arg.is_equal(arg.conjugate()) && is_a<numeric>(exp) && ex_to<numeric>(exp).is_even()) < return power(arg, exp); < else < return power(abs(arg), exp).hold(); --- > c.s << "{|"; seq[0].print(c); c.s << "|}"; 118,125d93 < REGISTER_FUNCTION(abs, eval_func(abs_eval). < evalf_func(abs_evalf). < print_func<print_latex>(abs_print_latex). < print_func<print_csrc_float>(abs_print_csrc_float). < print_func<print_csrc_double>(abs_print_csrc_float). < conjugate_func(abs_conjugate). < power_func(abs_power)); < 127c95 < // Step function --- > // complex conjugate 130c98,101 < static ex step_evalf(const ex & arg) --- > GINAC_IMPLEMENT_FUNCTION_OPT(conjugate_function, > print_func<print_latex>(&conjugate_function::do_print_latex)) > > ex conjugate_function::conjugate() const 132,135c103 < if (is_exactly_a<numeric>(arg)) < return step(ex_to<numeric>(arg)); < < return step(arg).hold(); --- > return seq[0]; 138c106 < static ex step_eval(const ex & arg) --- > ex conjugate_function::eval(int level) const 140,164c108 < if (is_exactly_a<numeric>(arg)) < return step(ex_to<numeric>(arg)); < < else if (is_exactly_a<mul>(arg) && < is_exactly_a<numeric>(arg.op(arg.nops()-1))) { < numeric oc = ex_to<numeric>(arg.op(arg.nops()-1)); < if (oc.is_real()) { < if (oc > 0) < // step(42*x) -> step(x) < return step(arg/oc).hold(); < else < // step(-42*x) -> step(-x) < return step(-arg/oc).hold(); < } < if (oc.real().is_zero()) { < if (oc.imag() > 0) < // step(42*I*x) -> step(I*x) < return step(I*arg/oc).hold(); < else < // step(-42*I*x) -> step(-I*x) < return step(-I*arg/oc).hold(); < } < } < < return step(arg).hold(); --- > return seq[0].conjugate(); 167,170c111 < static ex step_series(const ex & arg, < const relational & rel, < int order, < unsigned options) --- > ex conjugate_function::evalf(int level) const 172,180c113,117 < const ex arg_pt = arg.subs(rel, subs_options::no_pattern); < if (arg_pt.info(info_flags::numeric) < && ex_to<numeric>(arg_pt).real().is_zero() < && !(options & series_options::suppress_branchcut)) < throw (std::domain_error("step_series(): on imaginary axis")); < < epvector seq; < seq.push_back(expair(step(arg_pt), _ex0)); < return pseries(rel,seq); --- > const ex& arg = seq[0]; > if (is_exactly_a<numeric>(arg)) { > return ex_to<numeric>(arg).conjugate(); > } > return this->hold(); 183c120 < static ex step_conjugate(const ex& arg) --- > void conjugate_function::do_print_latex(const print_context& c, unsigned level) const 185c122,123 < return step(arg); --- > const ex& arg = seq[0]; > c.s << "\\bar{"; arg.print(c); c.s << "}"; 188,192d125 < REGISTER_FUNCTION(step, eval_func(step_eval). < evalf_func(step_evalf). < series_func(step_series). < conjugate_func(step_conjugate)); < 197c130,132 < static ex csgn_evalf(const ex & arg) --- > GINAC_IMPLEMENT_FUNCTION(csgn_function) > > ex csgn_function::conjugate() const 199,202c134 < if (is_exactly_a<numeric>(arg)) < return csgn(ex_to<numeric>(arg)); < < return csgn(arg).hold(); --- > return *this; 205c137 < static ex csgn_eval(const ex & arg) --- > ex csgn_function::eval(int level) const 206a139,140 > const ex& arg = seq[0]; > 231c165 < return csgn(arg).hold(); --- > return this->hold(); 234,237c168 < static ex csgn_series(const ex & arg, < const relational & rel, < int order, < unsigned options) --- > ex csgn_function::evalf(int level) const 239,243c170,172 < const ex arg_pt = arg.subs(rel, subs_options::no_pattern); < if (arg_pt.info(info_flags::numeric) < && ex_to<numeric>(arg_pt).real().is_zero() < && !(options & series_options::suppress_branchcut)) < throw (std::domain_error("csgn_series(): on imaginary axis")); --- > const ex& arg = seq[0]; > if (is_exactly_a<numeric>(arg)) > return csgn(ex_to<numeric>(arg)); 245,247c174 < epvector seq; < seq.push_back(expair(csgn(arg_pt), _ex0)); < return pseries(rel,seq); --- > return this->hold(); 250,255c177 < static ex csgn_conjugate(const ex& arg) < { < return csgn(arg); < } < < static ex csgn_power(const ex & arg, const ex & exp) --- > ex csgn_function::power_law(const ex& exp) const 256a179 > const ex& arg = seq[0]; 259c182 < return csgn(arg); --- > return *this; 261c184 < return power(csgn(arg), _ex2).hold(); --- > return power(*this, _ex2).hold(); 263,302c186 < return power(csgn(arg), exp).hold(); < } < < < REGISTER_FUNCTION(csgn, eval_func(csgn_eval). < evalf_func(csgn_evalf). < series_func(csgn_series). < conjugate_func(csgn_conjugate). < power_func(csgn_power)); < < < ////////// < // Eta function: eta(x,y) == log(x*y) - log(x) - log(y). < // This function is closely related to the unwinding number K, sometimes found < // in modern literature: K(z) == (z-log(exp(z)))/(2*Pi*I). < ////////// < < static ex eta_evalf(const ex &x, const ex &y) < { < // It seems like we basically have to replicate the eval function here, < // since the expression might not be fully evaluated yet. < if (x.info(info_flags::positive) || y.info(info_flags::positive)) < return _ex0; < < if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { < const numeric nx = ex_to<numeric>(x); < const numeric ny = ex_to<numeric>(y); < const numeric nxy = ex_to<numeric>(x*y); < int cut = 0; < if (nx.is_real() && nx.is_negative()) < cut -= 4; < if (ny.is_real() && ny.is_negative()) < cut -= 4; < if (nxy.is_real() && nxy.is_negative()) < cut += 4; < return evalf(I/4*Pi)*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)- < (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut); < } < < return eta(x,y).hold(); --- > return power(*this, exp).hold(); 305c189 < static ex eta_eval(const ex &x, const ex &y) --- > ex csgn_function::series(const relational& r, int order, unsigned options) const 307,325c191,196 < // trivial: eta(x,c) -> 0 if c is real and positive < if (x.info(info_flags::positive) || y.info(info_flags::positive)) < return _ex0; < < if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { < // don't call eta_evalf here because it would call Pi.evalf()! < const numeric nx = ex_to<numeric>(x); < const numeric ny = ex_to<numeric>(y); < const numeric nxy = ex_to<numeric>(x*y); < int cut = 0; < if (nx.is_real() && nx.is_negative()) < cut -= 4; < if (ny.is_real() && ny.is_negative()) < cut -= 4; < if (nxy.is_real() && nxy.is_negative()) < cut += 4; < return (I/4)*Pi*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)- < (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut); < } --- > const ex& arg = seq[0]; > const ex arg_pt = arg.subs(r, subs_options::no_pattern); > if (arg_pt.info(info_flags::numeric) > && ex_to<numeric>(arg_pt).real().is_zero() > && !(options & series_options::suppress_branchcut)) > throw (std::domain_error("csgn_series(): on imaginary axis")); 327,340d197 < return eta(x,y).hold(); < } < < static ex eta_series(const ex & x, const ex & y, < const relational & rel, < int order, < unsigned options) < { < const ex x_pt = x.subs(rel, subs_options::no_pattern); < const ex y_pt = y.subs(rel, subs_options::no_pattern); < if ((x_pt.info(info_flags::numeric) && x_pt.info(info_flags::negative)) || < (y_pt.info(info_flags::numeric) && y_pt.info(info_flags::negative)) || < ((x_pt*y_pt).info(info_flags::numeric) && (x_pt*y_pt).info(info_flags::negative))) < throw (std::domain_error("eta_series(): on discontinuity")); 342,348c199,200 < seq.push_back(expair(eta(x_pt,y_pt), _ex0)); < return pseries(rel,seq); < } < < static ex eta_conjugate(const ex & x, const ex & y) < { < return -eta(x,y); --- > seq.push_back(expair(csgn(arg_pt), _ex0)); > return pseries(r, seq); 351,358d202 < REGISTER_FUNCTION(eta, eval_func(eta_eval). < evalf_func(eta_evalf). < series_func(eta_series). < latex_name("\\eta"). < set_symmetry(sy_symm(0, 1)). < conjugate_func(eta_conjugate)); < < 360c204 < // dilogarithm --- > // Step function 363c207,209 < static ex Li2_evalf(const ex & x) --- > GINAC_IMPLEMENT_FUNCTION(step_function) > > ex step_function::conjugate() const 365,368c211 < if (is_exactly_a<numeric>(x)) < return Li2(ex_to<numeric>(x)); < < return Li2(x).hold(); --- > return *this; 371c214 < static ex Li2_eval(const ex & x) --- > ex step_function::eval(int level) const 373,463c216,230 < if (x.info(info_flags::numeric)) { < // Li2(0) -> 0 < if (x.is_zero()) < return _ex0; < // Li2(1) -> Pi^2/6 < if (x.is_equal(_ex1)) < return power(Pi,_ex2)/_ex6; < // Li2(1/2) -> Pi^2/12 - log(2)^2/2 < if (x.is_equal(_ex1_2)) < return power(Pi,_ex2)/_ex12 + power(log(_ex2),_ex2)*_ex_1_2; < // Li2(-1) -> -Pi^2/12 < if (x.is_equal(_ex_1)) < return -power(Pi,_ex2)/_ex12; < // Li2(I) -> -Pi^2/48+Catalan*I < if (x.is_equal(I)) < return power(Pi,_ex2)/_ex_48 + Catalan*I; < // Li2(-I) -> -Pi^2/48-Catalan*I < if (x.is_equal(-I)) < return power(Pi,_ex2)/_ex_48 - Catalan*I; < // Li2(float) < if (!x.info(info_flags::crational)) < return Li2(ex_to<numeric>(x)); < } < < return Li2(x).hold(); < } < < static ex Li2_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx Li2(x) -> -log(1-x)/x < return -log(_ex1-x)/x; < } < < static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options) < { < const ex x_pt = x.subs(rel, subs_options::no_pattern); < if (x_pt.info(info_flags::numeric)) { < // First special case: x==0 (derivatives have poles) < if (x_pt.is_zero()) { < // method: < // The problem is that in d/dx Li2(x==0) == -log(1-x)/x we cannot < // simply substitute x==0. The limit, however, exists: it is 1. < // We also know all higher derivatives' limits: < // (d/dx)^n Li2(x) == n!/n^2. < // So the primitive series expansion is < // Li2(x==0) == x + x^2/4 + x^3/9 + ... < // and so on. < // We first construct such a primitive series expansion manually in < // a dummy symbol s and then insert the argument's series expansion < // for s. Reexpanding the resulting series returns the desired < // result. < const symbol s; < ex ser; < // manually construct the primitive expansion < for (int i=1; i<order; ++i) < ser += pow(s,i) / pow(numeric(i), *_num2_p); < // 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); < // reexpanding it will collapse the series again < return ser.series(rel, order); < // NB: Of course, this still does not allow us to compute anything < // like sin(Li2(x)).series(x==0,2), since then this code here is < // not reached and the derivative of sin(Li2(x)) doesn't allow the < // substitution x==0. Probably limits *are* needed for the general < // cases. In case L'Hospital's rule is implemented for limits and < // basic::series() takes care of this, this whole block is probably < // obsolete! < } < // second special case: x==1 (branch point) < if (x_pt.is_equal(_ex1)) { < // method: < // construct series manually in a dummy symbol s < const symbol s; < ex ser = zeta(_ex2); < // manually construct the primitive expansion < for (int i=1; i<order; ++i) < ser += pow(1-s,i) * (numeric(1,i)*(I*Pi+log(s-1)) - numeric(1,i*i)); < // 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); < // reexpanding it will collapse the series again < return ser.series(rel, order); --- > const ex& arg = seq[0]; > > if (is_exactly_a<numeric>(arg)) > return step(ex_to<numeric>(arg)); > > else if (is_exactly_a<mul>(arg) && > is_exactly_a<numeric>(arg.op(arg.nops()-1))) { > numeric oc = ex_to<numeric>(arg.op(arg.nops()-1)); > if (oc.is_real()) { > if (oc > 0) > // step(42*x) -> step(x) > return step(arg/oc).hold(); > else > // step(-42*x) -> step(-x) > return step(-arg/oc).hold(); 465,483c232,238 < // third special case: x real, >=1 (branch cut) < if (!(options & series_options::suppress_branchcut) && < ex_to<numeric>(x_pt).is_real() && ex_to<numeric>(x_pt)>1) { < // method: < // This is the branch cut: assemble the primitive series manually < // and then add the corresponding complex step function. < const symbol &s = ex_to<symbol>(rel.lhs()); < const ex point = rel.rhs(); < const symbol foo; < epvector seq; < // zeroth order term: < seq.push_back(expair(Li2(x_pt), _ex0)); < // compute the intermediate terms: < ex replarg = series(Li2(x), s==foo, order); < for (size_t i=1; i<replarg.nops()-1; ++i) < seq.push_back(expair((replarg.op(i)/power(s-foo,i)).series(foo==point,1,options).op(0).subs(foo==s, subs_options::no_pattern),i)); < // append an order term: < seq.push_back(expair(Order(_ex1), replarg.nops()-1)); < return pseries(rel, seq); --- > if (oc.real().is_zero()) { > if (oc.imag() > 0) > // step(42*I*x) -> step(I*x) > return step(I*arg/oc).hold(); > else > // step(-42*I*x) -> step(-I*x) > return step(-I*arg/oc).hold(); 486,504c241,242 < // all other cases should be safe, by now: < throw do_taylor(); // caught by function::series() < } < < REGISTER_FUNCTION(Li2, eval_func(Li2_eval). < evalf_func(Li2_evalf). < derivative_func(Li2_deriv). < series_func(Li2_series). < latex_name("\\mbox{Li}_2")); < < ////////// < // trilogarithm < ////////// < < static ex Li3_eval(const ex & x) < { < if (x.is_zero()) < return x; < return Li3(x).hold(); --- > > return this->hold(); 507,514c245 < REGISTER_FUNCTION(Li3, eval_func(Li3_eval). < latex_name("\\mbox{Li}_3")); < < ////////// < // Derivatives of Riemann's Zeta-function zetaderiv(0,x)==zeta(x) < ////////// < < static ex zetaderiv_eval(const ex & n, const ex & x) --- > ex step_function::evalf(int level) const 516,520c247,249 < if (n.info(info_flags::numeric)) { < // zetaderiv(0,x) -> zeta(x) < if (n.is_zero()) < return zeta(x); < } --- > const ex& arg = seq[0]; > if (is_exactly_a<numeric>(arg)) > return step(ex_to<numeric>(arg)); 522c251 < return zetaderiv(n, x).hold(); --- > return this->hold(); 525c254 < static ex zetaderiv_deriv(const ex & n, const ex & x, unsigned deriv_param) --- > ex step_function::series(const relational& rel, int order, unsigned options) const 527c256,261 < GINAC_ASSERT(deriv_param<2); --- > const ex& arg = seq[0]; > const ex arg_pt = arg.subs(rel, subs_options::no_pattern); > if (arg_pt.info(info_flags::numeric) > && ex_to<numeric>(arg_pt).real().is_zero() > && !(options & series_options::suppress_branchcut)) > throw (std::domain_error("step_series(): on imaginary axis")); 529,534c263,265 < if (deriv_param==0) { < // d/dn zeta(n,x) < throw(std::logic_error("cannot diff zetaderiv(n,x) with respect to n")); < } < // d/dx psi(n,x) < return zetaderiv(n+1,x); --- > epvector seq; > seq.push_back(expair(step(arg_pt), _ex0)); > return pseries(rel,seq); 537,540d267 < REGISTER_FUNCTION(zetaderiv, eval_func(zetaderiv_eval). < derivative_func(zetaderiv_deriv). < latex_name("\\zeta^\\prime")); < 545c272,276 < static ex factorial_evalf(const ex & x) --- > GINAC_IMPLEMENT_FUNCTION_OPT(factorial_function, > print_func<print_dflt>(&factorial_function::do_print_dflt_latex). > print_func<print_latex>(&factorial_function::do_print_dflt_latex)) > > ex factorial_function::conjugate() const 547c278 < return factorial(x).hold(); --- > return *this; 550c281 < static ex factorial_eval(const ex & x) --- > ex factorial_function::eval(int level) const 551a283 > const ex& x = seq[0]; 555c287 < return factorial(x).hold(); --- > return this->hold(); 558c290 < static void factorial_print_dflt_latex(const ex & x, const print_context & c) --- > ex factorial_function::evalf(int level) const 560,562c292,298 < if (is_exactly_a<symbol>(x) || < is_exactly_a<constant>(x) || < is_exactly_a<function>(x)) { --- > return this->hold(); > } > > void factorial_function::do_print_dflt_latex(const print_context& c, unsigned level) const > { > const ex& x = seq[0]; > if (is_exactly_a<symbol>(x) || is_exactly_a<constant>(x) || is_exactly_a<function>(x)) { 569,579d304 < static ex factorial_conjugate(const ex & x) < { < return factorial(x); < } < < REGISTER_FUNCTION(factorial, eval_func(factorial_eval). < evalf_func(factorial_evalf). < print_func<print_dflt>(factorial_print_dflt_latex). < print_func<print_latex>(factorial_print_dflt_latex). < conjugate_func(factorial_conjugate)); < 584c309,314 < static ex binomial_evalf(const ex & x, const ex & y) --- > GINAC_IMPLEMENT_FUNCTION(binomial_function) > > // At the moment the numeric evaluation of a binomail function always > // gives a real number, but if this would be implemented using the gamma > // function, also complex conjugation should be changed (or rather, deleted). > ex binomial_function::conjugate() const 586c316 < return binomial(x, y).hold(); --- > return *this; 589c319 < static ex binomial_sym(const ex & x, const numeric & y) --- > ex binomial_function::sym(const ex& x, const numeric& y) const 604c334 < return binomial(x, y).hold(); --- > return this->hold(); 607c337 < static ex binomial_eval(const ex & x, const ex &y) --- > ex binomial_function::eval(int level) const 608a339,340 > const ex& x = seq[0]; > const ex& y = seq[1]; 613c345 < return binomial_sym(x, ex_to<numeric>(y)); --- > return sym(x, ex_to<numeric>(y)); 615c347 < return binomial(x, y).hold(); --- > return this->hold(); 618,621c350 < // At the moment the numeric evaluation of a binomail function always < // gives a real number, but if this would be implemented using the gamma < // function, also complex conjugation should be changed (or rather, deleted). < static ex binomial_conjugate(const ex & x, const ex & y) --- > ex binomial_function::evalf(int level) const 623c352 < return binomial(x,y); --- > return this->hold(); 626,629d354 < REGISTER_FUNCTION(binomial, eval_func(binomial_eval). < evalf_func(binomial_evalf). < conjugate_func(binomial_conjugate)); < 634c359,367 < static ex Order_eval(const ex & x) --- > GINAC_IMPLEMENT_FUNCTION_OPT(Order_function, > print_func<print_latex>(&Order_function::do_print_latex)) > > ex Order_function::conjugate() const > { > return *this; > } > > ex Order_function::derivative(const symbol& s) const 635a369,374 > return Order(seq[0].diff(s)); > } > > ex Order_function::eval(int level) const > { > const ex& x = op(0); 643c382 < const mul &m = ex_to<mul>(x); --- > const mul& m = ex_to<mul>(x); 648c387 < return Order(x).hold(); --- > return this->hold(); 651c390 < static ex Order_series(const ex & x, const relational & r, int order, unsigned options) --- > ex Order_function::series(const relational& r, int order, unsigned options) const 652a392 > const ex& x = op(0); 656c396 < const symbol &s = ex_to<symbol>(r.lhs()); --- > const symbol& s = ex_to<symbol>(r.lhs()); 661c401 < static ex Order_conjugate(const ex & x) --- > void Order_function::do_print_latex(const print_context& c, unsigned level) const 663c403,404 < return Order(x); --- > c.s << "\\mathcal{O}"; > inherited::do_print(c,level); 666c407,413 < // Differentiation is handled in function::derivative because of its special requirements --- > ////////// > // Abstract derivative of functions > ////////// > > GINAC_IMPLEMENT_FUNCTION_OPT(function_derivative_function, > print_func<print_dflt>(&function_derivative_function::do_print_dflt). > print_func<print_tree>(&function_derivative_function::do_print_tree)) 668,671c415,479 < REGISTER_FUNCTION(Order, eval_func(Order_eval). < series_func(Order_series). < latex_name("\\mathcal{O}"). < conjugate_func(Order_conjugate)); --- > /** Implementation of ex::diff() for derivatives. It applies the chain rule. > * @see ex::diff */ > ex function_derivative_function::derivative(const symbol& s) const > { > GINAC_ASSERT(seq[0].size() == 2); > GINAC_ASSERT(is_a<lst>(seq[0])); > GINAC_ASSERT(is_a<function>(seq[1])); > > const ex& func = seq[1]; > ex result; > for (size_t i=0; i<func.nops(); ++i) { > ex arg_diff = func.op(i).diff(s); > if (!arg_diff.is_zero()) { > lst params = ex_to<lst>(seq[0]); > params.append(i); > params.sort(); > result += arg_diff * function_derivative(params, func); > } > } > return result; > } > > ex function_derivative_function::eval(int level) const > { > if (level > 1) { > // first evaluate children, then we will end up here again > return function_derivative_function(evalchildren(level)); > } > > const ex& params = seq[0]; > // No parameters specified? Then return the function itself > if (params.nops() == 0) { > return seq[1]; > } > > if (params.nops() == 1) { > return ex_to<function>(seq[1]).pderivative((ex_to<numeric>(params.op(0))).to_int()); > } > > return this->hold(); > } > > void function_derivative_function::do_print_dflt(const print_context& c, unsigned level) const > { > c.s << "D["; > const lst& params = ex_to<lst>(seq[0]); > lst::const_iterator i = params.begin(), end = params.end(); > --end; > while (i != end) { > c.s << *i++ << ","; > } > c.s << *i << "](" << seq[1] << ")"; > } > > void function_derivative_function::do_print_tree(const print_tree& c, unsigned level) const > { > c.s << std::string(level, ' ') << class_name() << " " > << ex_to<function>(seq[1]).class_name() << " @" << this > << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec > << ", nops=" << nops() << std::endl; > seq[0].print(c, level + c.delta_indent); > for (size_t i=0; i<seq[1].nops(); ++i) > seq[1].op(i).print(c, level + c.delta_indent); > c.s << std::string(level + c.delta_indent, ' ') << "=====" << std::endl; > } 677c485 < ex lsolve(const ex &eqns, const ex &symbols, unsigned options) --- > ex lsolve(const ex& eqns, const ex& symbols, unsigned options) 758,759c566 < const numeric < fsolve(const ex& f_in, const symbol& x, const numeric& x1, const numeric& x2) --- > const numeric fsolve(const ex& f_in, const symbol& x, const numeric& x1, const numeric& x2) 853,858d659 < < /* Force inclusion of functions from inifcns_gamma and inifcns_zeta < * for static lib (so ginsh will see them). */ < unsigned force_include_tgamma = tgamma_SERIAL::serial; < unsigned force_include_zeta1 = zeta1_SERIAL::serial; < Index: ginac/inifcns.h =================================================================== RCS file: /home/cvs/GiNaC/ginac/inifcns.h,v retrieving revision 1.53 diff -r1.53 inifcns.h 6c6 < * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany --- > * GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany 32,34d31 < /** Complex conjugate. */ < DECLARE_FUNCTION_1P(conjugate_function) < 36,66c33,48 < DECLARE_FUNCTION_1P(abs) < < /** Step function. */ < DECLARE_FUNCTION_1P(step) < < /** Complex sign. */ < DECLARE_FUNCTION_1P(csgn) < < /** Eta function: log(a*b) == log(a) + log(b) + eta(a, b). */ < DECLARE_FUNCTION_2P(eta) < < /** Sine. */ < DECLARE_FUNCTION_1P(sin) < < /** Cosine. */ < DECLARE_FUNCTION_1P(cos) < < /** Tangent. */ < DECLARE_FUNCTION_1P(tan) < < /** Exponential function. */ < DECLARE_FUNCTION_1P(exp) < < /** Natural logarithm. */ < DECLARE_FUNCTION_1P(log) < < /** Inverse sine (arc sine). */ < DECLARE_FUNCTION_1P(asin) < < /** Inverse cosine (arc cosine). */ < DECLARE_FUNCTION_1P(acos) --- > class abs_function : public function > { > GINAC_DECLARE_FUNCTION_1P(abs_function) > public: > virtual ex conjugate() const; > virtual ex eval(int level = 0) const; > virtual ex evalf(int level = 0) const; > virtual ex power_law(const ex& exp) const; > protected: > void do_print_csrc_float(const print_context& c, unsigned level) const; > void do_print_latex(const print_context& c, unsigned level) const; > }; > > template<typename T1> inline abs_function abs(const T1& x1) { return abs_function(x1); } > inline abs_function abs(double x1) { return abs_function(x1); } > inline abs_function abs(float x1) { return abs_function(x1); } 68,69c50,55 < /** Inverse tangent (arc tangent). */ < DECLARE_FUNCTION_1P(atan) --- > /** Complex conjugate. */ > GINAC_FUNCTION_1P(conjugate, > GINAC_FUNCTION_conjugate > GINAC_FUNCTION_eval > GINAC_FUNCTION_evalf > GINAC_FUNCTION_print_latex) 71,72c57,63 < /** Inverse tangent with two arguments. */ < DECLARE_FUNCTION_2P(atan2) --- > /** Complex sign. */ > GINAC_FUNCTION_1P(csgn, > GINAC_FUNCTION_conjugate > GINAC_FUNCTION_eval > GINAC_FUNCTION_evalf > GINAC_FUNCTION_power_law > GINAC_FUNCTION_series) 74,75c65,74 < /** Hyperbolic Sine. */ < DECLARE_FUNCTION_1P(sinh) --- > /** Step function. */ > class step_function : public function > { > GINAC_DECLARE_FUNCTION_1P(step_function) > public: > virtual ex conjugate() const; > virtual ex eval(int level = 0) const; > virtual ex evalf(int level = 0) const; > virtual ex series(const relational& r, int order, unsigned options = 0) const; > }; 77,78c76 < /** Hyperbolic Cosine. */ < DECLARE_FUNCTION_1P(cosh) --- > template<typename T1> inline step_function step(const T1& x1) { return step_function(x1); } 80,81c78,88 < /** Hyperbolic Tangent. */ < DECLARE_FUNCTION_1P(tanh) --- > /** Factorial function. */ > class factorial_function : public function > { > GINAC_DECLARE_FUNCTION_1P(factorial_function) > public: > virtual ex conjugate() const; > virtual ex eval(int level = 0) const; > virtual ex evalf(int level = 0) const; > protected: > void do_print_dflt_latex(const print_context& c, unsigned level) const; > }; 83,84c90 < /** Inverse hyperbolic Sine (area hyperbolic sine). */ < DECLARE_FUNCTION_1P(asinh) --- > template<typename T1> inline factorial_function factorial(const T1& x1) { return factorial_function(x1); } 86,87c92,102 < /** Inverse hyperbolic Cosine (area hyperbolic cosine). */ < DECLARE_FUNCTION_1P(acosh) --- > /** Binomial function. */ > class binomial_function : public function > { > GINAC_DECLARE_FUNCTION_2P(binomial_function) > public: > virtual ex conjugate() const; > virtual ex eval(int level = 0) const; > virtual ex evalf(int level = 0) const; > protected: > ex sym(const ex& x, const numeric& y) const; > }; 89,90c104 < /** Inverse hyperbolic Tangent (area hyperbolic tangent). */ < DECLARE_FUNCTION_1P(atanh) --- > template<typename T1, typename T2> inline binomial_function binomial(const T1& x1, const T2& x2) { return binomial_function(x1, x2); } 92,93c106,117 < /** Dilogarithm. */ < DECLARE_FUNCTION_1P(Li2) --- > /** Order term function (for truncated power series). */ > class Order_function : public function > { > GINAC_DECLARE_FUNCTION_1P(Order_function) > public: > virtual ex conjugate() const; > virtual ex derivative(const symbol& s) const; > virtual ex eval(int level = 0) const; > virtual ex series(const relational& r, int order, unsigned options = 0) const; > protected: > void do_print_latex(const print_context& c, unsigned level) const; > }; 95,96c119 < /** Trilogarithm. */ < DECLARE_FUNCTION_1P(Li3) --- > template<typename T1> inline Order_function Order(const T1& x1) { return Order_function(x1); } 98,99c121,131 < /** Derivatives of Riemann's Zeta-function. */ < DECLARE_FUNCTION_2P(zetaderiv) --- > /** Abstract derivative of functions. */ > class function_derivative_function : public function > { > GINAC_DECLARE_FUNCTION_2P(function_derivative_function) > public: > virtual ex derivative(const symbol& s) const; > virtual ex eval(int level = 0) const; > protected: > void do_print_dflt(const print_context& c, unsigned level) const; > void do_print_tree(const print_tree& c, unsigned level) const; > }; 101,122d132 < // overloading at work: we cannot use the macros here < /** Multiple zeta value including Riemann's zeta-function. */ < class zeta1_SERIAL { public: static unsigned serial; }; < template<typename T1> < inline function zeta(const T1& p1) { < return function(zeta1_SERIAL::serial, ex(p1)); < } < /** Alternating Euler sum or colored MZV. */ < class zeta2_SERIAL { public: static unsigned serial; }; < template<typename T1, typename T2> < inline function zeta(const T1& p1, const T2& p2) { < return function(zeta2_SERIAL::serial, ex(p1), ex(p2)); < } < class zeta_SERIAL; < template<> inline bool is_the_function<zeta_SERIAL>(const ex& x) < { < return is_the_function<zeta1_SERIAL>(x) || is_the_function<zeta2_SERIAL>(x); < } < < // overloading at work: we cannot use the macros here < /** Generalized multiple polylogarithm. */ < class G2_SERIAL { public: static unsigned serial; }; 124,175c134 < inline function G(const T1& x, const T2& y) { < return function(G2_SERIAL::serial, ex(x), ex(y)); < } < /** Generalized multiple polylogarithm with explicit imaginary parts. */ < class G3_SERIAL { public: static unsigned serial; }; < template<typename T1, typename T2, typename T3> < inline function G(const T1& x, const T2& s, const T3& y) { < return function(G3_SERIAL::serial, ex(x), ex(s), ex(y)); < } < class G_SERIAL; < template<> inline bool is_the_function<G_SERIAL>(const ex& x) < { < return is_the_function<G2_SERIAL>(x) || is_the_function<G3_SERIAL>(x); < } < < /** Polylogarithm and multiple polylogarithm. */ < DECLARE_FUNCTION_2P(Li) < < /** Nielsen's generalized polylogarithm. */ < DECLARE_FUNCTION_3P(S) < < /** Harmonic polylogarithm. */ < DECLARE_FUNCTION_2P(H) < < /** Gamma-function. */ < DECLARE_FUNCTION_1P(lgamma) < DECLARE_FUNCTION_1P(tgamma) < < /** Beta-function. */ < DECLARE_FUNCTION_2P(beta) < < // overloading at work: we cannot use the macros here < /** Psi-function (aka digamma-function). */ < class psi1_SERIAL { public: static unsigned serial; }; < template<typename T1> < inline function psi(const T1 & p1) { < return function(psi1_SERIAL::serial, ex(p1)); < } < /** Derivatives of Psi-function (aka polygamma-functions). */ < class psi2_SERIAL { public: static unsigned serial; }; < template<typename T1, typename T2> < inline function psi(const T1 & p1, const T2 & p2) { < return function(psi2_SERIAL::serial, ex(p1), ex(p2)); < } < class psi_SERIAL; < template<> inline bool is_the_function<psi_SERIAL>(const ex & x) < { < return is_the_function<psi1_SERIAL>(x) || is_the_function<psi2_SERIAL>(x); < } < < /** Factorial function. */ < DECLARE_FUNCTION_1P(factorial) --- > inline function_derivative_function function_derivative(const T1& x1, const T2& x2) { return function_derivative_function(x1, x2); } 177,183c136 < /** Binomial function. */ < DECLARE_FUNCTION_2P(binomial) < < /** Order term function (for truncated power series). */ < DECLARE_FUNCTION_1P(Order) < < ex lsolve(const ex &eqns, const ex &symbols, unsigned options = solve_algo::automatic); --- > ex lsolve(const ex& eqns, const ex& symbols, unsigned options = solve_algo::automatic); 196,206d148 < /** Check whether a function is the Order (O(n)) function. */ < inline bool is_order_function(const ex & e) < { < return is_ex_the_function(e, Order); < } < < /** Converts a given list containing parameters for H in Remiddi/Vermaseren notation into < * the corresponding GiNaC functions. < */ < ex convert_H_to_Li(const ex& parameterlst, const ex& arg); < Index: ginac/inifcns_gamma.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/inifcns_gamma.cpp,v retrieving revision 1.64 diff -r1.64 inifcns_gamma.cpp 1,563d0 < /** @file inifcns_gamma.cpp < * < * Implementation of Gamma-function, Beta-function, Polygamma-functions, and < * some related stuff. */ < < /* < * GiNaC Copyright (C) 1999-2005 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #include <vector> < #include <stdexcept> < < #include "inifcns.h" < #include "constant.h" < #include "pseries.h" < #include "numeric.h" < #include "power.h" < #include "relational.h" < #include "operators.h" < #include "symbol.h" < #include "symmetry.h" < #include "utils.h" < < namespace GiNaC { < < ////////// < // Logarithm of Gamma function < ////////// < < static ex lgamma_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) { < try { < return lgamma(ex_to<numeric>(x)); < } catch (const dunno &e) { } < } < < return lgamma(x).hold(); < } < < < /** Evaluation of lgamma(x), the natural logarithm of the Gamma function. < * Knows about integer arguments and that's it. Somebody ought to provide < * some good numerical evaluation some day... < * < * @exception GiNaC::pole_error("lgamma_eval(): logarithmic pole",0) */ < static ex lgamma_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < // trap integer arguments: < if (x.info(info_flags::integer)) { < // lgamma(n) -> log((n-1)!) for postitive n < if (x.info(info_flags::posint)) < return log(factorial(x + _ex_1)); < else < throw (pole_error("lgamma_eval(): logarithmic pole",0)); < } < // lgamma_evalf should be called here once it becomes available < } < < return lgamma(x).hold(); < } < < < static ex lgamma_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx lgamma(x) -> psi(x) < return psi(x); < } < < < static ex lgamma_series(const ex & arg, < const relational & rel, < int order, < unsigned options) < { < // method: < // Taylor series where there is no pole falls back to psi function < // evaluation. < // On a pole at -m we could use the recurrence relation < // lgamma(x) == lgamma(x+1)-log(x) < // from which follows < // series(lgamma(x),x==-m,order) == < // series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order); < const ex arg_pt = arg.subs(rel, subs_options::no_pattern); < if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive)) < throw do_taylor(); // caught by function::series() < // if we got here we have to care for a simple pole of tgamma(-m): < numeric m = -ex_to<numeric>(arg_pt); < ex recur; < for (numeric p = 0; p<=m; ++p) < recur += log(arg+p); < return (lgamma(arg+m+_ex1)-recur).series(rel, order, options); < } < < < REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval). < evalf_func(lgamma_evalf). < derivative_func(lgamma_deriv). < series_func(lgamma_series). < latex_name("\\log \\Gamma")); < < < ////////// < // true Gamma function < ////////// < < static ex tgamma_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) { < try { < return tgamma(ex_to<numeric>(x)); < } catch (const dunno &e) { } < } < < return tgamma(x).hold(); < } < < < /** Evaluation of tgamma(x), the true Gamma function. Knows about integer < * arguments, half-integer arguments and that's it. Somebody ought to provide < * some good numerical evaluation some day... < * < * @exception pole_error("tgamma_eval(): simple pole",0) */ < static ex tgamma_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < // trap integer arguments: < const numeric two_x = (*_num2_p)*ex_to<numeric>(x); < if (two_x.is_even()) { < // tgamma(n) -> (n-1)! for postitive n < if (two_x.is_positive()) { < return factorial(ex_to<numeric>(x).sub(*_num1_p)); < } else { < throw (pole_error("tgamma_eval(): simple pole",1)); < } < } < // trap half integer arguments: < if (two_x.is_integer()) { < // trap positive x==(n+1/2) < // tgamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n) < if (two_x.is_positive()) { < const numeric n = ex_to<numeric>(x).sub(*_num1_2_p); < return (doublefactorial(n.mul(*_num2_p).sub(*_num1_p)).div(pow(*_num2_p,n))) * sqrt(Pi); < } else { < // trap negative x==(-n+1/2) < // tgamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1)) < const numeric n = abs(ex_to<numeric>(x).sub(*_num1_2_p)); < return (pow(*_num_2_p, n).div(doublefactorial(n.mul(*_num2_p).sub(*_num1_p))))*sqrt(Pi); < } < } < // tgamma_evalf should be called here once it becomes available < } < < return tgamma(x).hold(); < } < < < static ex tgamma_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx tgamma(x) -> psi(x)*tgamma(x) < return psi(x)*tgamma(x); < } < < < static ex tgamma_series(const ex & arg, < const relational & rel, < int order, < unsigned options) < { < // method: < // Taylor series where there is no pole falls back to psi function < // evaluation. < // On a pole at -m use the recurrence relation < // tgamma(x) == tgamma(x+1) / x < // from which follows < // series(tgamma(x),x==-m,order) == < // series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order); < const ex arg_pt = arg.subs(rel, subs_options::no_pattern); < if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive)) < throw do_taylor(); // caught by function::series() < // if we got here we have to care for a simple pole at -m: < const numeric m = -ex_to<numeric>(arg_pt); < ex ser_denom = _ex1; < for (numeric p; p<=m; ++p) < ser_denom *= arg+p; < return (tgamma(arg+m+_ex1)/ser_denom).series(rel, order, options); < } < < < REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval). < evalf_func(tgamma_evalf). < derivative_func(tgamma_deriv). < series_func(tgamma_series). < latex_name("\\Gamma")); < < < ////////// < // beta-function < ////////// < < static ex beta_evalf(const ex & x, const ex & y) < { < if (is_exactly_a<numeric>(x) && is_exactly_a<numeric>(y)) { < try { < return tgamma(ex_to<numeric>(x))*tgamma(ex_to<numeric>(y))/tgamma(ex_to<numeric>(x+y)); < } catch (const dunno &e) { } < } < < return beta(x,y).hold(); < } < < < static ex beta_eval(const ex & x, const ex & y) < { < if (x.is_equal(_ex1)) < return 1/y; < if (y.is_equal(_ex1)) < return 1/x; < if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { < // treat all problematic x and y that may not be passed into tgamma, < // because they would throw there although beta(x,y) is well-defined < // using the formula beta(x,y) == (-1)^y * beta(1-x-y, y) < const numeric &nx = ex_to<numeric>(x); < const numeric &ny = ex_to<numeric>(y); < if (nx.is_real() && nx.is_integer() && < ny.is_real() && ny.is_integer()) { < if (nx.is_negative()) { < if (nx<=-ny) < return pow(*_num_1_p, ny)*beta(1-x-y, y); < else < throw (pole_error("beta_eval(): simple pole",1)); < } < if (ny.is_negative()) { < if (ny<=-nx) < return pow(*_num_1_p, nx)*beta(1-y-x, x); < else < throw (pole_error("beta_eval(): simple pole",1)); < } < return tgamma(x)*tgamma(y)/tgamma(x+y); < } < // no problem in numerator, but denominator has pole: < if ((nx+ny).is_real() && < (nx+ny).is_integer() && < !(nx+ny).is_positive()) < return _ex0; < // beta_evalf should be called here once it becomes available < } < < return beta(x,y).hold(); < } < < < static ex beta_deriv(const ex & x, const ex & y, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param<2); < ex retval; < < // d/dx beta(x,y) -> (psi(x)-psi(x+y)) * beta(x,y) < if (deriv_param==0) < retval = (psi(x)-psi(x+y))*beta(x,y); < // d/dy beta(x,y) -> (psi(y)-psi(x+y)) * beta(x,y) < if (deriv_param==1) < retval = (psi(y)-psi(x+y))*beta(x,y); < return retval; < } < < < static ex beta_series(const ex & arg1, < const ex & arg2, < const relational & rel, < int order, < unsigned options) < { < // method: < // Taylor series where there is no pole of one of the tgamma functions < // falls back to beta function evaluation. Otherwise, fall back to < // tgamma series directly. < const ex arg1_pt = arg1.subs(rel, subs_options::no_pattern); < const ex arg2_pt = arg2.subs(rel, subs_options::no_pattern); < GINAC_ASSERT(is_a<symbol>(rel.lhs())); < const symbol &s = ex_to<symbol>(rel.lhs()); < ex arg1_ser, arg2_ser, arg1arg2_ser; < if ((!arg1_pt.info(info_flags::integer) || arg1_pt.info(info_flags::positive)) && < (!arg2_pt.info(info_flags::integer) || arg2_pt.info(info_flags::positive))) < throw do_taylor(); // caught by function::series() < // trap the case where arg1 is on a pole: < if (arg1.info(info_flags::integer) && !arg1.info(info_flags::positive)) < arg1_ser = tgamma(arg1+s); < else < arg1_ser = tgamma(arg1); < // trap the case where arg2 is on a pole: < if (arg2.info(info_flags::integer) && !arg2.info(info_flags::positive)) < arg2_ser = tgamma(arg2+s); < else < arg2_ser = tgamma(arg2); < // trap the case where arg1+arg2 is on a pole: < if ((arg1+arg2).info(info_flags::integer) && !(arg1+arg2).info(info_flags::positive)) < arg1arg2_ser = tgamma(arg2+arg1+s); < else < arg1arg2_ser = tgamma(arg2+arg1); < // compose the result (expanding all the terms): < return (arg1_ser*arg2_ser/arg1arg2_ser).series(rel, order, options).expand(); < } < < < REGISTER_FUNCTION(beta, eval_func(beta_eval). < evalf_func(beta_evalf). < derivative_func(beta_deriv). < series_func(beta_series). < latex_name("\\mbox{B}"). < set_symmetry(sy_symm(0, 1))); < < < ////////// < // Psi-function (aka digamma-function) < ////////// < < static ex psi1_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) { < try { < return psi(ex_to<numeric>(x)); < } catch (const dunno &e) { } < } < < return psi(x).hold(); < } < < /** Evaluation of digamma-function psi(x). < * Somebody ought to provide some good numerical evaluation some day... */ < static ex psi1_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < const numeric &nx = ex_to<numeric>(x); < if (nx.is_integer()) { < // integer case < if (nx.is_positive()) { < // psi(n) -> 1 + 1/2 +...+ 1/(n-1) - Euler < numeric rat = 0; < for (numeric i(nx+(*_num_1_p)); i>0; --i) < rat += i.inverse(); < return rat-Euler; < } else { < // for non-positive integers there is a pole: < throw (pole_error("psi_eval(): simple pole",1)); < } < } < if (((*_num2_p)*nx).is_integer()) { < // half integer case < if (nx.is_positive()) { < // psi((2m+1)/2) -> 2/(2m+1) + 2/2m +...+ 2/1 - Euler - 2log(2) < numeric rat = 0; < for (numeric i = (nx+(*_num_1_p))*(*_num2_p); i>0; i-=(*_num2_p)) < rat += (*_num2_p)*i.inverse(); < return rat-Euler-_ex2*log(_ex2); < } else { < // use the recurrence relation < // psi(-m-1/2) == psi(-m-1/2+1) - 1 / (-m-1/2) < // to relate psi(-m-1/2) to psi(1/2): < // psi(-m-1/2) == psi(1/2) + r < // where r == ((-1/2)^(-1) + ... + (-m-1/2)^(-1)) < numeric recur = 0; < for (numeric p = nx; p<0; ++p) < recur -= pow(p, *_num_1_p); < return recur+psi(_ex1_2); < } < } < // psi1_evalf should be called here once it becomes available < } < < return psi(x).hold(); < } < < static ex psi1_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx psi(x) -> psi(1,x) < return psi(_ex1, x); < } < < static ex psi1_series(const ex & arg, < const relational & rel, < int order, < unsigned options) < { < // method: < // Taylor series where there is no pole falls back to polygamma function < // evaluation. < // On a pole at -m use the recurrence relation < // psi(x) == psi(x+1) - 1/z < // from which follows < // series(psi(x),x==-m,order) == < // series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x==-m,order); < const ex arg_pt = arg.subs(rel, subs_options::no_pattern); < if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive)) < throw do_taylor(); // caught by function::series() < // if we got here we have to care for a simple pole at -m: < const numeric m = -ex_to<numeric>(arg_pt); < ex recur; < for (numeric p; p<=m; ++p) < recur += power(arg+p,_ex_1); < return (psi(arg+m+_ex1)-recur).series(rel, order, options); < } < < unsigned psi1_SERIAL::serial = < function::register_new(function_options("psi", 1). < eval_func(psi1_eval). < evalf_func(psi1_evalf). < derivative_func(psi1_deriv). < series_func(psi1_series). < latex_name("\\psi"). < overloaded(2)); < < ////////// < // Psi-functions (aka polygamma-functions) psi(0,x)==psi(x) < ////////// < < static ex psi2_evalf(const ex & n, const ex & x) < { < if (is_exactly_a<numeric>(n) && is_exactly_a<numeric>(x)) { < try { < return psi(ex_to<numeric>(n),ex_to<numeric>(x)); < } catch (const dunno &e) { } < } < < return psi(n,x).hold(); < } < < /** Evaluation of polygamma-function psi(n,x). < * Somebody ought to provide some good numerical evaluation some day... */ < static ex psi2_eval(const ex & n, const ex & x) < { < // psi(0,x) -> psi(x) < if (n.is_zero()) < return psi(x); < // psi(-1,x) -> log(tgamma(x)) < if (n.is_equal(_ex_1)) < return log(tgamma(x)); < if (n.info(info_flags::numeric) && n.info(info_flags::posint) && < x.info(info_flags::numeric)) { < const numeric &nn = ex_to<numeric>(n); < const numeric &nx = ex_to<numeric>(x); < if (nx.is_integer()) { < // integer case < if (nx.is_equal(*_num1_p)) < // use psi(n,1) == (-)^(n+1) * n! * zeta(n+1) < return pow(*_num_1_p,nn+(*_num1_p))*factorial(nn)*zeta(ex(nn+(*_num1_p))); < if (nx.is_positive()) { < // use the recurrence relation < // psi(n,m) == psi(n,m+1) - (-)^n * n! / m^(n+1) < // to relate psi(n,m) to psi(n,1): < // psi(n,m) == psi(n,1) + r < // where r == (-)^n * n! * (1^(-n-1) + ... + (m-1)^(-n-1)) < numeric recur = 0; < for (numeric p = 1; p<nx; ++p) < recur += pow(p, -nn+(*_num_1_p)); < recur *= factorial(nn)*pow((*_num_1_p), nn); < return recur+psi(n,_ex1); < } else { < // for non-positive integers there is a pole: < throw (pole_error("psi2_eval(): pole",1)); < } < } < if (((*_num2_p)*nx).is_integer()) { < // half integer case < if (nx.is_equal(*_num1_2_p)) < // use psi(n,1/2) == (-)^(n+1) * n! * (2^(n+1)-1) * zeta(n+1) < return pow(*_num_1_p,nn+(*_num1_p))*factorial(nn)*(pow(*_num2_p,nn+(*_num1_p)) + (*_num_1_p))*zeta(ex(nn+(*_num1_p))); < if (nx.is_positive()) { < const numeric m = nx - (*_num1_2_p); < // use the multiplication formula < // psi(n,2*m) == (psi(n,m) + psi(n,m+1/2)) / 2^(n+1) < // to revert to positive integer case < return psi(n,(*_num2_p)*m)*pow((*_num2_p),nn+(*_num1_p))-psi(n,m); < } else { < // use the recurrence relation < // psi(n,-m-1/2) == psi(n,-m-1/2+1) - (-)^n * n! / (-m-1/2)^(n+1) < // to relate psi(n,-m-1/2) to psi(n,1/2): < // psi(n,-m-1/2) == psi(n,1/2) + r < // where r == (-)^(n+1) * n! * ((-1/2)^(-n-1) + ... + (-m-1/2)^(-n-1)) < numeric recur = 0; < for (numeric p = nx; p<0; ++p) < recur += pow(p, -nn+(*_num_1_p)); < recur *= factorial(nn)*pow(*_num_1_p, nn+(*_num_1_p)); < return recur+psi(n,_ex1_2); < } < } < // psi2_evalf should be called here once it becomes available < } < < return psi(n, x).hold(); < } < < static ex psi2_deriv(const ex & n, const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param<2); < < if (deriv_param==0) { < // d/dn psi(n,x) < throw(std::logic_error("cannot diff psi(n,x) with respect to n")); < } < // d/dx psi(n,x) -> psi(n+1,x) < return psi(n+_ex1, x); < } < < static ex psi2_series(const ex & n, < const ex & arg, < const relational & rel, < int order, < unsigned options) < { < // method: < // Taylor series where there is no pole falls back to polygamma function < // evaluation. < // On a pole at -m use the recurrence relation < // psi(n,x) == psi(n,x+1) - (-)^n * n! / x^(n+1) < // from which follows < // series(psi(x),x==-m,order) == < // series(psi(x+m+1) - (-1)^n * n! * ((x)^(-n-1) + (x+1)^(-n-1) + ... < // ... + (x+m)^(-n-1))),x==-m,order); < const ex arg_pt = arg.subs(rel, subs_options::no_pattern); < if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive)) < throw do_taylor(); // caught by function::series() < // if we got here we have to care for a pole of order n+1 at -m: < const numeric m = -ex_to<numeric>(arg_pt); < ex recur; < for (numeric p; p<=m; ++p) < recur += power(arg+p,-n+_ex_1); < recur *= factorial(n)*power(_ex_1,n); < return (psi(n, arg+m+_ex1)-recur).series(rel, order, options); < } < < unsigned psi2_SERIAL::serial = < function::register_new(function_options("psi", 2). < eval_func(psi2_eval). < evalf_func(psi2_evalf). < derivative_func(psi2_deriv). < series_func(psi2_series). < latex_name("\\psi"). < overloaded(2)); < < < } // namespace GiNaC Index: ginac/inifcns_nstdsums.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/inifcns_nstdsums.cpp,v retrieving revision 1.28 diff -r1.28 inifcns_nstdsums.cpp 1,3948d0 < /** @file inifcns_nstdsums.cpp < * < * Implementation of some special functions that have a representation as nested sums. < * < * 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) < * 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)) < * < * Some remarks: < * < * - All formulae used can be looked up in the following publications: < * [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258. < * [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172. < * [ReV] Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754 < * [BBB] Special Values of Multiple Polylogarithms, J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941 < * [VSW] Numerical evaluation of multiple polylogarithms, J.Vollinga, S.Weinzierl, hep-ph/0410259 < * < * - The order of parameters and arguments of Li and zeta is defined according to the nested sums < * representation. The parameters for H are understood as in [ReV]. They can be in expanded --- only < * 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 < * 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. < * < * - The calculation of classical polylogarithms is speeded up by using Bernoulli numbers and < * look-up tables. S uses look-up tables as well. The zeta function applies the algorithms in < * [Cra] and [BBB] for speed up. Multiple polylogarithms use Hoelder convolution [BBB]. < * < * - The functions have no means to do a series expansion into nested sums. To do this, you have to convert < * these functions into the appropriate objects from the nestedsums library, do the expansion and convert < * the result back. < * < * - Numerical testing of this implementation has been performed by doing a comparison of results < * between this software and the commercial M.......... 4.1. Multiple zeta values have been checked < * by means of evaluations into simple zeta values. Harmonic polylogarithms have been checked by < * comparison to S(n,p,x) for corresponding parameter combinations and by continuity checks < * around |x|=1 along with comparisons to corresponding zeta functions. Multiple polylogarithms were < * checked against H and zeta and by means of shuffle and quasi-shuffle relations. < * < */ < < /* < * GiNaC Copyright (C) 1999-2005 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #include <sstream> < #include <stdexcept> < #include <vector> < #include <cln/cln.h> < < #include "inifcns.h" < < #include "add.h" < #include "constant.h" < #include "lst.h" < #include "mul.h" < #include "numeric.h" < #include "operators.h" < #include "power.h" < #include "pseries.h" < #include "relational.h" < #include "symbol.h" < #include "utils.h" < #include "wildcard.h" < < < namespace GiNaC { < < < ////////////////////////////////////////////////////////////////////// < // < // Classical polylogarithm Li(n,x) < // < // helper functions < // < ////////////////////////////////////////////////////////////////////// < < < // anonymous namespace for helper functions < namespace { < < < // lookup table for factors built from Bernoulli numbers < // see fill_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; < int xnsize = 0; < < < // This function calculates the X_n. The X_n are needed for speed up of classical polylogarithms. < // With these numbers the polylogs can be calculated as follows: < // Li_p (x) = \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with u = -log(1-x) < // X_0(n) = B_n (Bernoulli numbers) < // X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k) < // The calculation of Xn depends on X0 and X{n-1}. < // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater. < // This results in a slightly more complicated algorithm for the X_n. < // The first index in Xn corresponds to the index of the polylog minus 2. < // The second index in Xn corresponds to the index from the actual sum. < void fill_Xn(int n) < { < if (n>1) { < // calculate X_2 and higher (corresponding to Li_4 and higher) < std::vector<cln::cl_N> buf(xninitsize); < std::vector<cln::cl_N>::iterator it = buf.begin(); < cln::cl_N result; < *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1 < it++; < for (int i=2; i<=xninitsize; i++) { < if (i&1) { < result = 0; // k == 0 < } else { < result = Xn[0][i/2-1]; // k == 0 < } < for (int k=1; k<i-1; k++) { < if ( !(((i-k) & 1) && ((i-k) > 1)) ) { < result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1); < } < } < result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1 < result = result + Xn[n-1][i-1] / (i+1); // k == i < < *it = result; < it++; < } < Xn.push_back(buf); < } else if (n==1) { < // special case to handle the X_0 correct < std::vector<cln::cl_N> buf(xninitsize); < std::vector<cln::cl_N>::iterator it = buf.begin(); < cln::cl_N result; < *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1 < it++; < *it = cln::cl_I(17)/cln::cl_I(36); // i == 2 < it++; < for (int i=3; i<=xninitsize; i++) { < if (i & 1) { < result = -Xn[0][(i-3)/2]/2; < *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result; < it++; < } else { < result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1); < for (int k=1; k<i/2; k++) { < result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1); < } < *it = result; < it++; < } < } < Xn.push_back(buf); < } else { < // calculate X_0 < std::vector<cln::cl_N> buf(xninitsize/2); < std::vector<cln::cl_N>::iterator it = buf.begin(); < for (int i=1; i<=xninitsize/2; i++) { < *it = bernoulli(i*2).to_cl_N(); < it++; < } < Xn.push_back(buf); < } < < xnsize++; < } < < // doubles the number of entries in each Xn[] < void double_Xn() < { < const int pos0 = xninitsize / 2; < // X_0 < for (int i=1; i<=xninitsizestep/2; ++i) { < Xn[0].push_back(bernoulli((i+pos0)*2).to_cl_N()); < } < if (Xn.size() > 1) { < int xend = xninitsize + xninitsizestep; < cln::cl_N result; < // X_1 < for (int i=xninitsize+1; i<=xend; ++i) { < if (i & 1) { < result = -Xn[0][(i-3)/2]/2; < Xn[1].push_back((cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result); < } else { < result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1); < for (int k=1; k<i/2; k++) { < result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1); < } < Xn[1].push_back(result); < } < } < // X_n < for (int n=2; n<Xn.size(); ++n) { < for (int i=xninitsize+1; i<=xend; ++i) { < if (i & 1) { < result = 0; // k == 0 < } else { < result = Xn[0][i/2-1]; // k == 0 < } < for (int k=1; k<i-1; ++k) { < if ( !(((i-k) & 1) && ((i-k) > 1)) ) { < result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1); < } < } < result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1 < result = result + Xn[n-1][i-1] / (i+1); // k == i < Xn[n].push_back(result); < } < } < } < xninitsize += xninitsizestep; < } < < < // calculates Li(2,x) without Xn < cln::cl_N Li2_do_sum(const cln::cl_N& x) < { < cln::cl_N res = x; < cln::cl_N resbuf; < cln::cl_N num = x * cln::cl_float(1, cln::float_format(Digits)); < cln::cl_I den = 1; // n^2 = 1 < unsigned i = 3; < do { < resbuf = res; < num = num * x; < den = den + i; // n^2 = 4, 9, 16, ... < i += 2; < res = res + num / den; < } while (res != resbuf); < return res; < } < < < // calculates Li(2,x) with Xn < cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x) < { < std::vector<cln::cl_N>::const_iterator it = Xn[0].begin(); < std::vector<cln::cl_N>::const_iterator xend = Xn[0].end(); < cln::cl_N u = -cln::log(1-x); < cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits)); < cln::cl_N uu = cln::square(u); < cln::cl_N res = u - uu/4; < cln::cl_N resbuf; < unsigned i = 1; < do { < resbuf = res; < factor = factor * uu / (2*i * (2*i+1)); < res = res + (*it) * factor; < i++; < if (++it == xend) { < double_Xn(); < it = Xn[0].begin() + (i-1); < xend = Xn[0].end(); < } < } while (res != resbuf); < return res; < } < < < // calculates Li(n,x), n>2 without Xn < cln::cl_N Lin_do_sum(int n, const cln::cl_N& x) < { < cln::cl_N factor = x * cln::cl_float(1, cln::float_format(Digits)); < cln::cl_N res = x; < cln::cl_N resbuf; < int i=2; < do { < resbuf = res; < factor = factor * x; < res = res + factor / cln::expt(cln::cl_I(i),n); < i++; < } while (res != resbuf); < return res; < } < < < // calculates Li(n,x), n>2 with Xn < cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x) < { < std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin(); < std::vector<cln::cl_N>::const_iterator xend = Xn[n-2].end(); < cln::cl_N u = -cln::log(1-x); < cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits)); < cln::cl_N res = u; < cln::cl_N resbuf; < unsigned i=2; < do { < resbuf = res; < factor = factor * u / i; < res = res + (*it) * factor; < i++; < if (++it == xend) { < double_Xn(); < it = Xn[n-2].begin() + (i-2); < xend = Xn[n-2].end(); < } < } while (res != resbuf); < return res; < } < < < // forward declaration needed by function Li_projection and C below < numeric S_num(int n, int p, const numeric& x); < < < // helper function for classical polylog Li < cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec) < { < // treat n=2 as special case < if (n == 2) { < // check if precalculated X0 exists < if (xnsize == 0) { < fill_Xn(0); < } < < if (cln::realpart(x) < 0.5) { < // choose the faster algorithm < // the switching point was empirically determined. the optimal point < // depends on hardware, Digits, ... so an approx value is okay. < // it solves also the problem with precision due to the u=-log(1-x) transformation < if (cln::abs(cln::realpart(x)) < 0.25) { < < return Li2_do_sum(x); < } else { < return Li2_do_sum_Xn(x); < } < } else { < // choose the faster algorithm < if (cln::abs(cln::realpart(x)) > 0.75) { < return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2); < } else { < return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2); < } < } < } else { < // check if precalculated Xn exist < if (n > xnsize+1) { < for (int i=xnsize; i<n-1; i++) { < fill_Xn(i); < } < } < < if (cln::realpart(x) < 0.5) { < // choose the faster algorithm < // with n>=12 the "normal" summation always wins against the method with Xn < if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) { < return Lin_do_sum(n, x); < } else { < return Lin_do_sum_Xn(n, x); < } < } 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()) < * 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) < { < if (n == 1) { < // just a log < return -cln::log(1-x.to_cl_N()); < } < if (x.is_zero()) { < return 0; < } < if (x == 1) { < // [Kol] (2.22) < return cln::zeta(n); < } < else 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); < 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); < } < 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(); < // 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))); < < // [Kol] (5.15) < if (cln::abs(value) > 1) { < cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n); < // check if argument is complex. if it is real, the new polylog has to be conjugated. < if (cln::zerop(cln::imagpart(value))) { < if (n & 1) { < result = result + conjugate(Li_projection(n, cln::recip(value), prec)); < } < else { < result = result - conjugate(Li_projection(n, cln::recip(value), prec)); < } < } < else { < if (n & 1) { < result = result + Li_projection(n, cln::recip(value), prec); < } < else { < result = result - Li_projection(n, cln::recip(value), prec); < } < } < 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); < } < result = result - add; < return result; < } < else { < return Li_projection(n, value, prec); < } < } < < < } // end of anonymous namespace < < < ////////////////////////////////////////////////////////////////////// < // < // Multiple polylogarithm Li(n,x) < // < // helper function < // < ////////////////////////////////////////////////////////////////////// < < < // anonymous namespace for helper function < namespace { < < < // performs the actual series summation for multiple polylogarithms < cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x) < { < const int j = s.size(); < < std::vector<cln::cl_N> t(j); < cln::cl_F one = cln::cl_float(1, cln::float_format(Digits)); < < cln::cl_N t0buf; < int q = 0; < do { < t0buf = t[0]; < // do it once ... < 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]); < } < // ... and do it again (to avoid premature drop out due to special arguments) < 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]); < } < } while (t[0] != t0buf); < < return t[0]; < } < < < // converts parameter types and calls multipleLi_do_sum (convenience function for G_numeric) < cln::cl_N mLi_do_summation(const lst& m, const lst& x) < { < std::vector<int> m_int; < std::vector<cln::cl_N> x_cln; < for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) { < m_int.push_back(ex_to<numeric>(*itm).to_int()); < x_cln.push_back(ex_to<numeric>(*itx).to_cl_N()); < } < return multipleLi_do_sum(m_int, x_cln); < } < < < // forward declaration for Li_eval() < lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf); < < < // holding dummy-symbols for the G/Li transformations < std::vector<ex> gsyms; < < < // type used by the transformation functions for G < typedef std::vector<int> Gparameter; < < < // G_eval1-function for G transformations < ex G_eval1(int a, int scale) < { < if (a != 0) { < const ex& scs = gsyms[std::abs(scale)]; < const ex& as = gsyms[std::abs(a)]; < if (as != scs) { < return -log(1 - scs/as); < } else { < return -zeta(1); < } < } else { < return log(gsyms[std::abs(scale)]); < } < } < < < // G_eval-function for G transformations < ex G_eval(const Gparameter& a, int scale) < { < // check for properties of G < ex sc = gsyms[std::abs(scale)]; < lst newa; < bool all_zero = true; < bool all_ones = true; < int count_ones = 0; < for (Gparameter::const_iterator it = a.begin(); it != a.end(); ++it) { < if (*it != 0) { < const ex sym = gsyms[std::abs(*it)]; < newa.append(sym); < all_zero = false; < if (sym != sc) { < all_ones = false; < } < if (all_ones) { < ++count_ones; < } < } else { < all_ones = false; < } < } < < // care about divergent G: shuffle to separate divergencies that will be canceled < // later on in the transformation < if (newa.nops() > 1 && newa.op(0) == sc && !all_ones && a.front()!=0) { < // do shuffle < Gparameter short_a; < Gparameter::const_iterator it = a.begin(); < ++it; < for (; it != a.end(); ++it) { < short_a.push_back(*it); < } < ex result = G_eval1(a.front(), scale) * G_eval(short_a, scale); < it = short_a.begin(); < for (int i=1; i<count_ones; ++i) { < ++it; < } < for (; it != short_a.end(); ++it) { < < Gparameter newa; < Gparameter::const_iterator it2 = short_a.begin(); < for (--it2; it2 != it;) { < ++it2; < newa.push_back(*it2); < } < newa.push_back(a[0]); < ++it2; < for (; it2 != short_a.end(); ++it2) { < newa.push_back(*it2); < } < result -= G_eval(newa, scale); < } < return result / count_ones; < } < < // G({1,...,1};y) -> G({1};y)^k / k! < if (all_ones && a.size() > 1) { < return pow(G_eval1(a.front(),scale), count_ones) / factorial(count_ones); < } < < // G({0,...,0};y) -> log(y)^k / k! < if (all_zero) { < return pow(log(gsyms[std::abs(scale)]), a.size()) / factorial(a.size()); < } < < // no special cases anymore -> convert it into Li < lst m; < lst x; < ex argbuf = gsyms[std::abs(scale)]; < ex mval = _ex1; < for (Gparameter::const_iterator it=a.begin(); it!=a.end(); ++it) { < if (*it != 0) { < const ex& sym = gsyms[std::abs(*it)]; < x.append(argbuf / sym); < m.append(mval); < mval = _ex1; < argbuf = sym; < } else { < ++mval; < } < } < return pow(-1, x.nops()) * Li(m, x); < } < < < // converts data for G: pending_integrals -> a < Gparameter convert_pending_integrals_G(const Gparameter& pending_integrals) < { < GINAC_ASSERT(pending_integrals.size() != 1); < < if (pending_integrals.size() > 0) { < // get rid of the first element, which would stand for the new upper limit < Gparameter new_a(pending_integrals.begin()+1, pending_integrals.end()); < return new_a; < } else { < // just return empty parameter list < Gparameter new_a; < return new_a; < } < } < < < // check the parameters a and scale for G and return information about convergence, depth, etc. < // convergent : true if G(a,scale) is convergent < // depth : depth of G(a,scale) < // trailing_zeros : number of trailing zeros of a < // min_it : iterator of a pointing on the smallest element in a < Gparameter::const_iterator check_parameter_G(const Gparameter& a, int scale, < bool& convergent, int& depth, int& trailing_zeros, Gparameter::const_iterator& min_it) < { < convergent = true; < depth = 0; < trailing_zeros = 0; < min_it = a.end(); < Gparameter::const_iterator lastnonzero = a.end(); < for (Gparameter::const_iterator it = a.begin(); it != a.end(); ++it) { < if (std::abs(*it) > 0) { < ++depth; < trailing_zeros = 0; < lastnonzero = it; < if (std::abs(*it) < scale) { < convergent = false; < if ((min_it == a.end()) || (std::abs(*it) < std::abs(*min_it))) { < min_it = it; < } < } < } else { < ++trailing_zeros; < } < } < return ++lastnonzero; < } < < < // add scale to pending_integrals if pending_integrals is empty < Gparameter prepare_pending_integrals(const Gparameter& pending_integrals, int scale) < { < GINAC_ASSERT(pending_integrals.size() != 1); < < if (pending_integrals.size() > 0) { < return pending_integrals; < } else { < Gparameter new_pending_integrals; < new_pending_integrals.push_back(scale); < return new_pending_integrals; < } < } < < < // handles trailing zeroes for an otherwise convergent integral < ex trailing_zeros_G(const Gparameter& a, int scale) < { < bool convergent; < int depth, trailing_zeros; < Gparameter::const_iterator last, dummyit; < last = check_parameter_G(a, scale, convergent, depth, trailing_zeros, dummyit); < < GINAC_ASSERT(convergent); < < if ((trailing_zeros > 0) && (depth > 0)) { < ex result; < Gparameter new_a(a.begin(), a.end()-1); < result += G_eval1(0, scale) * trailing_zeros_G(new_a, scale); < for (Gparameter::const_iterator it = a.begin(); it != last; ++it) { < Gparameter new_a(a.begin(), it); < new_a.push_back(0); < new_a.insert(new_a.end(), it, a.end()-1); < result -= trailing_zeros_G(new_a, scale); < } < < return result / trailing_zeros; < } else { < return G_eval(a, scale); < } < } < < < // G transformation [VSW] (57),(58) < ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, int scale) < { < // pendint = ( y1, b1, ..., br ) < // a = ( 0, ..., 0, amin ) < // scale = y2 < // < // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(0, ..., 0, sr; y2) < // where sr replaces amin < < GINAC_ASSERT(a.back() != 0); < GINAC_ASSERT(a.size() > 0); < < ex result; < Gparameter new_pending_integrals = prepare_pending_integrals(pending_integrals, std::abs(a.back())); < const int psize = pending_integrals.size(); < < // length == 1 < // G(sr_{+-}; y2 ) = G(y2_{-+}; sr) - G(0; sr) + ln(-y2_{-+}) < < if (a.size() == 1) { < < // ln(-y2_{-+}) < result += log(gsyms[ex_to<numeric>(scale).to_int()]); < if (a.back() > 0) { < new_pending_integrals.push_back(-scale); < result += I*Pi; < } else { < new_pending_integrals.push_back(scale); < result -= I*Pi; < } < if (psize) { < result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front()); < } < < // G(y2_{-+}; sr) < result += trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front()); < < // G(0; sr) < new_pending_integrals.back() = 0; < result -= trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front()); < < return result; < } < < // length > 1 < // G_m(sr_{+-}; y2) = -zeta_m + int_0^y2 dt/t G_{m-1}( (1/y2)_{+-}; 1/t ) < // - int_0^sr dt/t G_{m-1}( (1/y2)_{+-}; 1/t ) < < //term zeta_m < result -= zeta(a.size()); < if (psize) { < result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front()); < } < < // term int_0^sr dt/t G_{m-1}( (1/y2)_{+-}; 1/t ) < // = int_0^sr dt/t G_{m-1}( t_{+-}; y2 ) < Gparameter new_a(a.begin()+1, a.end()); < new_pending_integrals.push_back(0); < result -= depth_one_trafo_G(new_pending_integrals, new_a, scale); < < // term int_0^y2 dt/t G_{m-1}( (1/y2)_{+-}; 1/t ) < // = int_0^y2 dt/t G_{m-1}( t_{+-}; y2 ) < Gparameter new_pending_integrals_2; < new_pending_integrals_2.push_back(scale); < new_pending_integrals_2.push_back(0); < if (psize) { < result += trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front()) < * depth_one_trafo_G(new_pending_integrals_2, new_a, scale); < } else { < result += depth_one_trafo_G(new_pending_integrals_2, new_a, scale); < } < < return result; < } < < < // forward declaration < ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2, < const Gparameter& pendint, const Gparameter& a_old, int scale); < < < // G transformation [VSW] < ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) < { < // main recursion routine < // < // pendint = ( y1, b1, ..., br ) < // a = ( a1, ..., amin, ..., aw ) < // scale = y2 < // < // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,sr,...,aw,y2) < // where sr replaces amin < < // find smallest alpha, determine depth and trailing zeros, and check for convergence < bool convergent; < int depth, trailing_zeros; < Gparameter::const_iterator min_it; < Gparameter::const_iterator firstzero = < check_parameter_G(a, scale, convergent, depth, trailing_zeros, min_it); < int min_it_pos = min_it - a.begin(); < < // special case: all a's are zero < if (depth == 0) { < ex result; < < if (a.size() == 0) { < result = 1; < } else { < result = G_eval(a, scale); < } < if (pendint.size() > 0) { < result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front()); < } < return result; < } < < // handle trailing zeros < if (trailing_zeros > 0) { < ex result; < Gparameter new_a(a.begin(), a.end()-1); < result += G_eval1(0, scale) * G_transform(pendint, new_a, scale); < for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) { < Gparameter new_a(a.begin(), it); < new_a.push_back(0); < new_a.insert(new_a.end(), it, a.end()-1); < result -= G_transform(pendint, new_a, scale); < } < return result / trailing_zeros; < } < < // convergence case < if (convergent) { < if (pendint.size() > 0) { < return G_eval(convert_pending_integrals_G(pendint), pendint.front()) * G_eval(a, scale); < } else { < return G_eval(a, scale); < } < } < < // call basic transformation for depth equal one < if (depth == 1) { < return depth_one_trafo_G(pendint, a, scale); < } < < // do recursion < // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,sr,...,aw,y2) < // = int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,0,...,aw,y2) < // + int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) int_0^{sr} ds_{r+1} d/ds_{r+1} G(a1,...,s_{r+1},...,aw,y2) < < // smallest element in last place < if (min_it + 1 == a.end()) { < do { --min_it; } while (*min_it == 0); < Gparameter empty; < Gparameter a1(a.begin(),min_it+1); < Gparameter a2(min_it+1,a.end()); < < ex result = G_transform(pendint,a2,scale)*G_transform(empty,a1,scale); < < result -= shuffle_G(empty,a1,a2,pendint,a,scale); < return result; < } < < Gparameter empty; < Gparameter::iterator changeit; < < // first term G(a_1,..,0,...,a_w;a_0) < Gparameter new_pendint = prepare_pending_integrals(pendint, a[min_it_pos]); < Gparameter new_a = a; < new_a[min_it_pos] = 0; < ex result = G_transform(empty, new_a, scale); < if (pendint.size() > 0) { < result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front()); < } < < // other terms < changeit = new_a.begin() + min_it_pos; < changeit = new_a.erase(changeit); < if (changeit != new_a.begin()) { < // smallest in the middle < new_pendint.push_back(*changeit); < result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front()) < * G_transform(empty, new_a, scale); < int buffer = *changeit; < *changeit = *min_it; < result += G_transform(new_pendint, new_a, scale); < *changeit = buffer; < new_pendint.pop_back(); < --changeit; < new_pendint.push_back(*changeit); < result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front()) < * G_transform(empty, new_a, scale); < *changeit = *min_it; < result -= G_transform(new_pendint, new_a, scale); < } else { < // smallest at the front < new_pendint.push_back(scale); < result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front()) < * G_transform(empty, new_a, scale); < new_pendint.back() = *changeit; < result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front()) < * G_transform(empty, new_a, scale); < *changeit = *min_it; < result += G_transform(new_pendint, new_a, scale); < } < return result; < } < < < // shuffles the two parameter list a1 and a2 and calls G_transform for every term except < // for the one that is equal to a_old < ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2, < const Gparameter& pendint, const Gparameter& a_old, int scale) < { < if (a1.size()==0 && a2.size()==0) { < // veto the one configuration we don't want < if ( a0 == a_old ) return 0; < < return G_transform(pendint,a0,scale); < } < < if (a2.size()==0) { < Gparameter empty; < Gparameter aa0 = a0; < aa0.insert(aa0.end(),a1.begin(),a1.end()); < return shuffle_G(aa0,empty,empty,pendint,a_old,scale); < } < < if (a1.size()==0) { < Gparameter empty; < Gparameter aa0 = a0; < aa0.insert(aa0.end(),a2.begin(),a2.end()); < return shuffle_G(aa0,empty,empty,pendint,a_old,scale); < } < < Gparameter a1_removed(a1.begin()+1,a1.end()); < Gparameter a2_removed(a2.begin()+1,a2.end()); < < Gparameter a01 = a0; < Gparameter a02 = a0; < < a01.push_back( a1[0] ); < a02.push_back( a2[0] ); < < return shuffle_G(a01,a1_removed,a2,pendint,a_old,scale) < + shuffle_G(a02,a1,a2_removed,pendint,a_old,scale); < } < < < // handles the transformations and the numerical evaluation of G < // the parameter x, s and y must only contain numerics < ex G_numeric(const lst& x, const lst& s, const ex& y) < { < // check for convergence and necessary accelerations < bool need_trafo = false; < bool need_hoelder = false; < int depth = 0; < for (lst::const_iterator it = x.begin(); it != x.end(); ++it) { < if (!(*it).is_zero()) { < ++depth; < if (abs(*it) - y < -pow(10,-Digits+2)) { < need_trafo = true; < break; < } < if (abs((abs(*it) - y)/y) < 0.01) { < need_hoelder = true; < } < } < } < if (x.op(x.nops()-1).is_zero()) { < need_trafo = true; < } < if (depth == 1 && !need_trafo) { < return -Li(x.nops(), y / x.op(x.nops()-1)).evalf(); < } < < // convergence transformation < if (need_trafo) { < < // sort (|x|<->position) to determine indices < std::multimap<ex,int> sortmap; < int size = 0; < for (int i=0; i<x.nops(); ++i) { < if (!x[i].is_zero()) { < sortmap.insert(std::pair<ex,int>(abs(x[i]), i)); < ++size; < } < } < // include upper limit (scale) < sortmap.insert(std::pair<ex,int>(abs(y), x.nops())); < < // generate missing dummy-symbols < int i = 1; < gsyms.clear(); < gsyms.push_back(symbol("GSYMS_ERROR")); < ex lastentry; < for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) { < if (it != sortmap.begin()) { < if (it->second < x.nops()) { < if (x[it->second] == lastentry) { < gsyms.push_back(gsyms.back()); < continue; < } < } else { < if (y == lastentry) { < gsyms.push_back(gsyms.back()); < continue; < } < } < } < std::ostringstream os; < os << "a" << i; < gsyms.push_back(symbol(os.str())); < ++i; < if (it->second < x.nops()) { < lastentry = x[it->second]; < } else { < lastentry = y; < } < } < < // fill position data according to sorted indices and prepare substitution list < Gparameter a(x.nops()); < lst subslst; < int pos = 1; < int scale; < for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) { < if (it->second < x.nops()) { < if (s[it->second] > 0) { < a[it->second] = pos; < } else { < a[it->second] = -pos; < } < subslst.append(gsyms[pos] == x[it->second]); < } else { < scale = pos; < subslst.append(gsyms[pos] == y); < } < ++pos; < } < < // do transformation < Gparameter pendint; < ex result = G_transform(pendint, a, scale); < // replace dummy symbols with their values < result = result.eval().expand(); < result = result.subs(subslst).evalf(); < < return result; < } < < // do acceleration transformation (hoelder convolution [BBB]) < if (need_hoelder) { < < ex result; < const int size = x.nops(); < lst newx; < for (lst::const_iterator it = x.begin(); it != x.end(); ++it) { < newx.append(*it / y); < } < < for (int r=0; r<=size; ++r) { < ex buffer = pow(-1, r); < ex p = 2; < bool adjustp; < do { < adjustp = false; < for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) { < if (*it == 1/p) { < p += (3-p)/2; < adjustp = true; < continue; < } < } < } while (adjustp); < ex q = p / (p-1); < lst qlstx; < lst qlsts; < for (int j=r; j>=1; --j) { < qlstx.append(1-newx.op(j-1)); < if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) { < qlsts.append( s.op(j-1)); < } else { < qlsts.append( -s.op(j-1)); < } < } < if (qlstx.nops() > 0) { < buffer *= G_numeric(qlstx, qlsts, 1/q); < } < lst plstx; < lst plsts; < for (int j=r+1; j<=size; ++j) { < plstx.append(newx.op(j-1)); < plsts.append(s.op(j-1)); < } < if (plstx.nops() > 0) { < buffer *= G_numeric(plstx, plsts, 1/p); < } < result += buffer; < } < return result; < } < < // do summation < lst newx; < lst m; < int mcount = 1; < ex sign = 1; < ex factor = y; < for (lst::const_iterator it = x.begin(); it != x.end(); ++it) { < if ((*it).is_zero()) { < ++mcount; < } else { < newx.append(factor / (*it)); < factor = *it; < m.append(mcount); < mcount = 1; < sign = -sign; < } < } < < return sign * numeric(mLi_do_summation(m, newx)); < } < < < ex mLi_numeric(const lst& m, const lst& x) < { < // let G_numeric do the transformation < lst newx; < lst s; < ex factor = 1; < for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) { < for (int i = 1; i < *itm; ++i) { < newx.append(0); < s.append(1); < } < newx.append(factor / *itx); < factor /= *itx; < s.append(1); < } < return pow(-1, m.nops()) * G_numeric(newx, s, _ex1); < } < < < } // end of anonymous namespace < < < ////////////////////////////////////////////////////////////////////// < // < // Generalized multiple polylogarithm G(x, y) and G(x, s, y) < // < // GiNaC function < // < ////////////////////////////////////////////////////////////////////// < < < static ex G2_evalf(const ex& x_, const ex& y) < { < if (!y.info(info_flags::positive)) { < return G(x_, y).hold(); < } < lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_); < if (x.nops() == 0) { < return _ex1; < } < if (x.op(0) == y) { < return G(x_, y).hold(); < } < lst s; < bool all_zero = true; < for (lst::const_iterator it = x.begin(); it != x.end(); ++it) { < if (!(*it).info(info_flags::numeric)) { < return G(x_, y).hold(); < } < if (*it != _ex0) { < all_zero = false; < } < s.append(+1); < } < if (all_zero) { < return pow(log(y), x.nops()) / factorial(x.nops()); < } < return G_numeric(x, s, y); < } < < < static ex G2_eval(const ex& x_, const ex& y) < { < //TODO eval to MZV or H or S or Lin < < if (!y.info(info_flags::positive)) { < return G(x_, y).hold(); < } < lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_); < if (x.nops() == 0) { < return _ex1; < } < if (x.op(0) == y) { < return G(x_, y).hold(); < } < lst s; < bool all_zero = true; < bool crational = true; < for (lst::const_iterator it = x.begin(); it != x.end(); ++it) { < if (!(*it).info(info_flags::numeric)) { < return G(x_, y).hold(); < } < if (!(*it).info(info_flags::crational)) { < crational = false; < } < if (*it != _ex0) { < all_zero = false; < } < s.append(+1); < } < if (all_zero) { < return pow(log(y), x.nops()) / factorial(x.nops()); < } < if (!y.info(info_flags::crational)) { < crational = false; < } < if (crational) { < return G(x_, y).hold(); < } < return G_numeric(x, s, y); < } < < < 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). < // print_func<print_latex>(G2_print_latex). < < < static ex G3_evalf(const ex& x_, const ex& s_, const ex& y) < { < if (!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_); < if (x.nops() != s.nops()) { < return G(x_, s_, y).hold(); < } < if (x.nops() == 0) { < return _ex1; < } < if (x.op(0) == y) { < return G(x_, s_, y).hold(); < } < lst sn; < bool all_zero = true; < for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) { < if (!(*itx).info(info_flags::numeric)) { < return G(x_, y).hold(); < } < if (!(*its).info(info_flags::real)) { < return G(x_, y).hold(); < } < if (*itx != _ex0) { < all_zero = false; < } < if (*its >= 0) { < sn.append(+1); < } else { < sn.append(-1); < } < } < if (all_zero) { < return pow(log(y), x.nops()) / factorial(x.nops()); < } < return G_numeric(x, sn, y); < } < < < static ex G3_eval(const ex& x_, const ex& s_, const ex& y) < { < //TODO eval to MZV or H or S or Lin < < if (!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_); < if (x.nops() != s.nops()) { < return G(x_, s_, y).hold(); < } < if (x.nops() == 0) { < return _ex1; < } < if (x.op(0) == y) { < return G(x_, s_, y).hold(); < } < lst sn; < bool all_zero = true; < bool crational = true; < for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) { < if (!(*itx).info(info_flags::numeric)) { < return G(x_, s_, y).hold(); < } < if (!(*its).info(info_flags::real)) { < return G(x_, s_, y).hold(); < } < if (!(*itx).info(info_flags::crational)) { < crational = false; < } < if (*itx != _ex0) { < all_zero = false; < } < if (*its >= 0) { < sn.append(+1); < } else { < sn.append(-1); < } < } < if (all_zero) { < return pow(log(y), x.nops()) / factorial(x.nops()); < } < if (!y.info(info_flags::crational)) { < crational = false; < } < if (crational) { < return G(x_, s_, y).hold(); < } < return G_numeric(x, sn, y); < } < < < 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). < // print_func<print_latex>(G3_print_latex). < < < ////////////////////////////////////////////////////////////////////// < // < // Classical polylogarithm and multiple polylogarithm Li(m,x) < // < // GiNaC function < // < ////////////////////////////////////////////////////////////////////// < < < static ex Li_evalf(const ex& m_, const ex& x_) < { < // 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_)); < } 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)); < } < } < } < // multiple polylogs < if (is_a<lst>(m_) && is_a<lst>(x_)) { < < const lst& m = ex_to<lst>(m_); < const lst& x = ex_to<lst>(x_); < if (m.nops() != x.nops()) { < return Li(m_,x_).hold(); < } < if (x.nops() == 0) { < return _ex1; < } < if ((m.op(0) == _ex1) && (x.op(0) == _ex1)) { < return Li(m_,x_).hold(); < } < < for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) { < if (!(*itm).info(info_flags::posint)) { < return Li(m_, x_).hold(); < } < if (!(*itx).info(info_flags::numeric)) { < return Li(m_, x_).hold(); < } < if (*itx == _ex0) { < return _ex0; < } < } < < return mLi_numeric(m, x); < } < < return Li(m_,x_).hold(); < } < < < static ex Li_eval(const ex& m_, const ex& x_) < { < if (is_a<lst>(m_)) { < if (is_a<lst>(x_)) { < // multiple polylogs < const lst& m = ex_to<lst>(m_); < const lst& x = ex_to<lst>(x_); < if (m.nops() != x.nops()) { < return Li(m_,x_).hold(); < } < if (x.nops() == 0) { < return _ex1; < } < bool is_H = true; < bool is_zeta = true; < bool do_evalf = true; < bool crational = true; < for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) { < if (!(*itm).info(info_flags::posint)) { < return Li(m_,x_).hold(); < } < if ((*itx != _ex1) && (*itx != _ex_1)) { < if (itx != x.begin()) { < is_H = false; < } < is_zeta = false; < } < if (*itx == _ex0) { < return _ex0; < } < if (!(*itx).info(info_flags::numeric)) { < do_evalf = false; < } < if (!(*itx).info(info_flags::crational)) { < crational = false; < } < } < if (is_zeta) { < return zeta(m_,x_); < } < if (is_H) { < ex prefactor; < lst newm = convert_parameter_Li_to_H(m, x, prefactor); < return prefactor * H(newm, x[0]); < } < if (do_evalf && !crational) { < return mLi_numeric(m,x); < } < } < return Li(m_, x_).hold(); < } else if (is_a<lst>(x_)) { < return Li(m_, x_).hold(); < } < < // classical polylogs < if (x_ == _ex0) { < return _ex0; < } < if (x_ == _ex1) { < return zeta(m_); < } < if (x_ == _ex_1) { < return (pow(2,1-m_)-1) * zeta(m_); < } < if (m_ == _ex1) { < return -log(1-x_); < } < if (m_ == _ex2) { < if (x_.is_equal(I)) { < return power(Pi,_ex2)/_ex_48 + Catalan*I; < } < if (x_.is_equal(-I)) { < return power(Pi,_ex2)/_ex_48 - Catalan*I; < } < } < 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_)); < } < < return Li(m_, x_).hold(); < } < < < static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options) < { < 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); < } < < // classical polylog < const ex x_pt = x.subs(rel, subs_options::no_pattern); < if (m.info(info_flags::numeric) && x_pt.info(info_flags::numeric)) { < // First special case: x==0 (derivatives have poles) < if (x_pt.is_zero()) { < const symbol s; < ex ser; < // manually construct the primitive expansion < for (int i=1; i<order; ++i) < ser += pow(s,i) / pow(numeric(i), m); < // 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); < // reexpanding it will collapse the series again < return ser.series(rel, order); < } < // TODO special cases: x==1 (branch point) and x real, >=1 (branch cut) < throw std::runtime_error("Li_series: don't know how to do the series expansion at this point!"); < } < // all other cases should be safe, by now: < throw do_taylor(); // caught by function::series() < } < < < static ex Li_deriv(const ex& m_, const ex& x_, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param < 2); < if (deriv_param == 0) { < return _ex0; < } < if (m_.nops() > 1) { < throw std::runtime_error("don't know how to derivate multiple polylogarithm!"); < } < ex m; < if (is_a<lst>(m_)) { < m = m_.op(0); < } else { < m = m_; < } < ex x; < if (is_a<lst>(x_)) { < x = x_.op(0); < } else { < x = x_; < } < if (m > 0) { < return Li(m-1, x) / x; < } else { < return 1/(1-x); < } < } < < < static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c) < { < lst m; < if (is_a<lst>(m_)) { < m = ex_to<lst>(m_); < } else { < m = lst(m_); < } < lst x; < if (is_a<lst>(x_)) { < x = ex_to<lst>(x_); < } else { < x = lst(x_); < } < c.s << "\\mbox{Li}_{"; < lst::const_iterator itm = m.begin(); < (*itm).print(c); < itm++; < for (; itm != m.end(); itm++) { < c.s << ","; < (*itm).print(c); < } < c.s << "}("; < lst::const_iterator itx = x.begin(); < (*itx).print(c); < itx++; < for (; itx != x.end(); itx++) { < c.s << ","; < (*itx).print(c); < } < c.s << ")"; < } < < < REGISTER_FUNCTION(Li, < evalf_func(Li_evalf). < eval_func(Li_eval). < series_func(Li_series). < derivative_func(Li_deriv). < print_func<print_latex>(Li_print_latex). < do_not_evalf_params()); < < < ////////////////////////////////////////////////////////////////////// < // < // Nielsen's generalized polylogarithm S(n,p,x) < // < // helper functions < // < ////////////////////////////////////////////////////////////////////// < < < // anonymous namespace for helper functions < namespace { < < < // lookup table for special Euler-Zagier-Sums (used for S_n,p(x)) < // see fill_Yn() < std::vector<std::vector<cln::cl_N> > Yn; < int ynsize = 0; // number of Yn[] < int ynlength = 100; // initial length of all Yn[i] < < < // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x). < // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum < // representing S_{n,p}(x). < // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the < // equivalent Z-sum. < // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum < // representing S_{n,p}(x). < // The calculation of Y_n uses the values from Y_{n-1}. < void fill_Yn(int n, const cln::float_format_t& prec) < { < const int initsize = ynlength; < //const int initsize = initsize_Yn; < cln::cl_N one = cln::cl_float(1, prec); < < if (n) { < std::vector<cln::cl_N> buf(initsize); < std::vector<cln::cl_N>::iterator it = buf.begin(); < std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin(); < *it = (*itprev) / cln::cl_N(n+1) * one; < it++; < itprev++; < // sums with an index smaller than the depth are zero and need not to be calculated. < // calculation starts with depth, which is n+2) < for (int i=n+2; i<=initsize+n; i++) { < *it = *(it-1) + (*itprev) / cln::cl_N(i) * one; < it++; < itprev++; < } < Yn.push_back(buf); < } else { < std::vector<cln::cl_N> buf(initsize); < std::vector<cln::cl_N>::iterator it = buf.begin(); < *it = 1 * one; < it++; < for (int i=2; i<=initsize; i++) { < *it = *(it-1) + 1 / cln::cl_N(i) * one; < it++; < } < Yn.push_back(buf); < } < ynsize++; < } < < < // make Yn longer ... < void make_Yn_longer(int newsize, const cln::float_format_t& prec) < { < < cln::cl_N one = cln::cl_float(1, prec); < < Yn[0].resize(newsize); < std::vector<cln::cl_N>::iterator it = Yn[0].begin(); < it += ynlength; < for (int i=ynlength+1; i<=newsize; i++) { < *it = *(it-1) + 1 / cln::cl_N(i) * one; < it++; < } < < for (int n=1; n<ynsize; n++) { < Yn[n].resize(newsize); < std::vector<cln::cl_N>::iterator it = Yn[n].begin(); < std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin(); < it += ynlength; < itprev += ynlength; < for (int i=ynlength+n+1; i<=newsize+n; i++) { < *it = *(it-1) + (*itprev) / cln::cl_N(i) * one; < it++; < itprev++; < } < } < < ynlength = newsize; < } < < < // helper function for S(n,p,x) < // [Kol] (7.2) < cln::cl_N C(int n, int p) < { < cln::cl_N result; < < for (int k=0; k<p; k++) { < for (int j=0; j<=(n+k-1)/2; j++) { < 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); < } < else { < result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j); < } < } < } < else { < 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::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::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() < / (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::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j)); < } < } < } < } < } < int np = n+p; < if ((np-1) & 1) { < if (((np)/2+n) & 1) { < result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p)); < } < else { < result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p)); < } < } < < return result; < } < < < // helper function for S(n,p,x) < // [Kol] remark to (9.1) < cln::cl_N a_k(int k) < { < cln::cl_N result; < < if (k == 0) { < return 1; < } < < result = result; < for (int m=2; m<=k; m++) { < result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m); < } < < return -result / k; < } < < < // helper function for S(n,p,x) < // [Kol] remark to (9.1) < cln::cl_N b_k(int k) < { < cln::cl_N result; < < if (k == 0) { < return 1; < } < < result = result; < for (int m=2; m<=k; m++) { < result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m); < } < < return result / k; < } < < < // 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) < { < if (p==1) { < return Li_projection(n+1, x, prec); < } < < // check if precalculated values are sufficient < if (p > ynsize+1) { < for (int i=ynsize; i<p-1; i++) { < fill_Yn(i, prec); < } < } < < // should be done otherwise < cln::cl_F one = cln::cl_float(1, cln::float_format(Digits)); < cln::cl_N xf = x * one; < //cln::cl_N xf = x * cln::cl_float(1, prec); < < cln::cl_N res; < cln::cl_N resbuf; < cln::cl_N factor = cln::expt(xf, p); < int i = p; < do { < resbuf = res; < if (i-p >= ynlength) { < // make Yn longer < make_Yn_longer(ynlength*2, prec); < } < res = res + factor / cln::expt(cln::cl_I(i),n+1) * Yn[p-2][i-p]; // should we check it? or rely on magic number? ... < //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ... < factor = factor * xf; < i++; < } while (res != resbuf); < < return res; < } < < < // helper function for S(n,p,x) < cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec) < { < // [Kol] (5.3) < if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) { < < cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n) < * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p); < < for (int s=0; s<n; s++) { < 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-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); < } < < return result; < } < < return S_do_sum(n, p, x, prec); < } < < < // helper function for S(n,p,x) < numeric S_num(int n, int p, const numeric& x) < { < if (x == 1) { < if (n == 1) { < // [Kol] (2.22) with (2.21) < return cln::zeta(p+1); < } < < if (p == 1) { < // [Kol] (2.22) < return cln::zeta(n+1); < } < < // [Kol] (9.1) < cln::cl_N result; < for (int nu=0; nu<n; nu++) { < for (int rho=0; rho<=p; rho++) { < result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1) < * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1); < } < } < result = result * cln::expt(cln::cl_I(-1),n+p-1); < < return result; < } < else if (x == -1) { < // [Kol] (2.22) < if (p == 1) { < return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1); < } < // throw std::runtime_error("don't know how to evaluate this function!"); < } < < // 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))); < < // [Kol] (5.3) < if ((cln::realpart(value) < -0.5) || (n == 0)) { < < cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n) < * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p); < < for (int s=0; s<n; s++) { < 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); < } < result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s); < } < < return result; < < } < // [Kol] (5.12) < if (cln::abs(value) > 1) { < < cln::cl_N result; < < for (int s=0; s<p; s++) { < 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(); < } < } < result = result * cln::expt(cln::cl_I(-1),n); < < cln::cl_N res2; < for (int r=0; r<n; r++) { < res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r); < } < res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p); < < result = result + cln::expt(cln::cl_I(-1),p) * res2; < < return result; < } < else { < return S_projection(n, p, value, prec); < } < } < < < } // end of anonymous namespace < < < ////////////////////////////////////////////////////////////////////// < // < // Nielsen's generalized polylogarithm S(n,p,x) < // < // GiNaC function < // < ////////////////////////////////////////////////////////////////////// < < < static ex S_evalf(const ex& n, const ex& p, const ex& x) < { < if (n.info(info_flags::posint) && p.info(info_flags::posint)) { < if (is_a<numeric>(x)) { < return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x)); < } 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)); < } < } < } < return S(n, p, x).hold(); < } < < < static ex S_eval(const ex& n, const ex& p, const ex& x) < { < if (n.info(info_flags::posint) && p.info(info_flags::posint)) { < if (x == 0) { < return _ex0; < } < if (x == 1) { < lst m(n+1); < for (int i=ex_to<numeric>(p).to_int()-1; i>0; i--) { < m.append(1); < } < return zeta(m); < } < if (p == 1) { < 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)); < } < } < if (n.is_zero()) { < // [Kol] (5.3) < return pow(-log(1-x), p) / factorial(p); < } < return S(n, p, x).hold(); < } < < < static ex S_series(const ex& n, const ex& p, const ex& x, const relational& rel, int order, unsigned options) < { < if (p == _ex1) { < return Li(n+1, x).series(rel, order, options); < } < < const ex x_pt = x.subs(rel, subs_options::no_pattern); < if (n.info(info_flags::posint) && p.info(info_flags::posint) && x_pt.info(info_flags::numeric)) { < // First special case: x==0 (derivatives have poles) < if (x_pt.is_zero()) { < const symbol s; < ex ser; < // manually construct the primitive expansion < // subsum = Euler-Zagier-Sum is needed < // dirty hack (slow ...) calculation of subsum: < std::vector<ex> presubsum, subsum; < subsum.push_back(0); < for (int i=1; i<order-1; ++i) { < subsum.push_back(subsum[i-1] + numeric(1, i)); < } < for (int depth=2; depth<p; ++depth) { < presubsum = subsum; < for (int i=1; i<order-1; ++i) { < subsum[i] = subsum[i-1] + numeric(1, i) * presubsum[i-1]; < } < } < < for (int i=1; i<order; ++i) { < ser += pow(s,i) / pow(numeric(i), n+1) * subsum[i-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); < // reexpanding it will collapse the series again < return ser.series(rel, order); < } < // TODO special cases: x==1 (branch point) and x real, >=1 (branch cut) < throw std::runtime_error("S_series: don't know how to do the series expansion at this point!"); < } < // all other cases should be safe, by now: < throw do_taylor(); // caught by function::series() < } < < < static ex S_deriv(const ex& n, const ex& p, const ex& x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param < 3); < if (deriv_param < 2) { < return _ex0; < } < if (n > 0) { < return S(n-1, p, x) / x; < } else { < return S(n, p-1, x) / (1-x); < } < } < < < static void S_print_latex(const ex& n, const ex& p, const ex& x, const print_context& c) < { < c.s << "\\mbox{S}_{"; < n.print(c); < c.s << ","; < p.print(c); < c.s << "}("; < x.print(c); < c.s << ")"; < } < < < REGISTER_FUNCTION(S, < evalf_func(S_evalf). < eval_func(S_eval). < series_func(S_series). < derivative_func(S_deriv). < print_func<print_latex>(S_print_latex). < do_not_evalf_params()); < < < ////////////////////////////////////////////////////////////////////// < // < // Harmonic polylogarithm H(m,x) < // < // helper functions < // < ////////////////////////////////////////////////////////////////////// < < < // anonymous namespace for helper functions < namespace { < < < // regulates the pole (used by 1/x-transformation) < symbol H_polesign("IMSIGN"); < < < // convert parameters from H to Li representation < // parameters are expected to be in expanded form, i.e. only 0, 1 and -1 < // returns true if some parameters are negative < bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf) < { < // expand parameter list < lst mexp; < for (lst::const_iterator it = l.begin(); it != l.end(); it++) { < if (*it > 1) { < for (ex count=*it-1; count > 0; count--) { < mexp.append(0); < } < mexp.append(1); < } else if (*it < -1) { < for (ex count=*it+1; count < 0; count++) { < mexp.append(0); < } < mexp.append(-1); < } else { < mexp.append(*it); < } < } < < ex signum = 1; < pf = 1; < bool has_negative_parameters = false; < ex acc = 1; < for (lst::const_iterator it = mexp.begin(); it != mexp.end(); it++) { < if (*it == 0) { < acc++; < continue; < } < if (*it > 0) { < m.append((*it+acc-1) * signum); < } else { < m.append((*it-acc+1) * signum); < } < acc = 1; < signum = *it; < pf *= *it; < if (pf < 0) { < has_negative_parameters = true; < } < } < if (has_negative_parameters) { < for (int i=0; i<m.nops(); i++) { < if (m.op(i) < 0) { < m.let_op(i) = -m.op(i); < s.append(-1); < } else { < s.append(1); < } < } < } < < return has_negative_parameters; < } < < < // recursivly transforms H to corresponding multiple polylogarithms < struct map_trafo_H_convert_to_Li : 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); < < lst m; < lst s; < ex pf; < if (convert_parameter_H_to_Li(parameter, m, s, pf)) { < s.let_op(0) = s.op(0) * arg; < return pf * Li(m, s).hold(); < } else { < for (int i=0; i<m.nops(); i++) { < s.append(1); < } < s.let_op(0) = s.op(0) * arg; < return Li(m, s).hold(); < } < } < } < return e; < } < }; < < < // recursivly transforms H to corresponding zetas < struct map_trafo_H_convert_to_zeta : 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)); < } < < 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(); < int signum = 1; < pf = _ex1; < res.append(*itm); < itm++; < while (itx != x.end()) { < signum *= (*itx > 0) ? 1 : -1; < 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] (49) < // put parameter 1 in front of existing parameters < ex trafo_H_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()); < } else { < return e * H(lst(1),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_1mx : 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) { < lst newparameter; < for (int i=parameter.nops(); i>0; i--) { < newparameter.append(0); < } < return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold(); < } < } else if (parameter.op(0) == -1) { < throw std::runtime_error("map_trafo_H_1mx: cannot handle weights equal -1!"); < } else { < for (int i=1; i<parameter.nops(); i++) { < if (parameter.op(i) != 1) { < allthesame = false; < break; < } < } < if (allthesame) { < lst newparameter; < for (int i=parameter.nops(); i>0; i--) { < newparameter.append(1); < } < return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold(); < } < } < < lst newparameter = parameter; < newparameter.remove_first(); < < if (parameter.op(0) == 0) { < < // leading zero < ex res = convert_H_to_zeta(parameter); < //ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1))); < map_trafo_H_1mx recursion; < ex buffer = recursion(H(newparameter, arg).hold()); < if (is_a<add>(buffer)) { < for (int i=0; i<buffer.nops(); i++) { < res -= trafo_H_prepend_one(buffer.op(i), arg); < } < } else { < res -= trafo_H_prepend_one(buffer, arg); < } < return res; < < } else { < < // leading one < map_trafo_H_1mx recursion; < map_trafo_H_mult unify; < ex res; < 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(); < } < return (unify((-H(lst(0), 1-arg).hold() * recursion(H(newparameter, arg).hold())).expand()) + < recursion(res)) / firstzero; < < } < < } < } < return e; < } < }; < < < // 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() + H_polesign, 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 < // < ////////////////////////////////////////////////////////////////////// < < < static ex H_evalf(const ex& x1, const ex& x2) < { < if (is_a<lst>(x1)) { < < cln::cl_N x; < if (is_a<numeric>(x2)) { < x = ex_to<numeric>(x2).to_cl_N(); < } else { < ex x2_val = x2.evalf(); < if (is_a<numeric>(x2_val)) { < x = ex_to<numeric>(x2_val).to_cl_N(); < } < } < < for (int i=0; i<x1.nops(); i++) { < if (!x1.op(i).info(info_flags::integer)) { < return H(x1, x2).hold(); < } < } < if (x1.nops() < 1) { < return H(x1, x2).hold(); < } < < const lst& morg = ex_to<lst>(x1); < // remove trailing zeros ... < if (*(--morg.end()) == 0) { < symbol xtemp("xtemp"); < map_trafo_H_reduce_trailing_zeros filter; < return filter(H(x1, xtemp).hold()).subs(xtemp==x2).evalf(); < } < // ... and expand parameter notation < bool has_minus_one = false; < lst m; < for (lst::const_iterator it = morg.begin(); it != morg.end(); it++) { < if (*it > 1) { < for (ex count=*it-1; count > 0; count--) { < m.append(0); < } < m.append(1); < } else if (*it <= -1) { < for (ex count=*it+1; count < 0; count++) { < m.append(0); < } < m.append(-1); < has_minus_one = true; < } else { < m.append(*it); < } < } < < // do summation < if (cln::abs(x) < 0.95) { < lst m_lst; < lst s_lst; < ex pf; < if (convert_parameter_H_to_Li(m, m_lst, s_lst, pf)) { < // negative parameters -> s_lst is filled < std::vector<int> m_int; < std::vector<cln::cl_N> x_cln; < for (lst::const_iterator it_int = m_lst.begin(), it_cln = s_lst.begin(); < it_int != m_lst.end(); it_int++, it_cln++) { < m_int.push_back(ex_to<numeric>(*it_int).to_int()); < x_cln.push_back(ex_to<numeric>(*it_cln).to_cl_N()); < } < x_cln.front() = x_cln.front() * x; < return pf * numeric(multipleLi_do_sum(m_int, x_cln)); < } else { < // only positive parameters < //TODO < if (m_lst.nops() == 1) { < return Li(m_lst.op(0), x2).evalf(); < } < std::vector<int> m_int; < for (lst::const_iterator it = m_lst.begin(); it != m_lst.end(); it++) { < m_int.push_back(ex_to<numeric>(*it).to_int()); < } < return numeric(H_do_sum(m_int, x)); < } < } < < symbol xtemp("xtemp"); < ex res = 1; < < // ensure that the realpart of the argument is positive < if (cln::realpart(x) < 0) { < x = -x; < for (int i=0; i<m.nops(); i++) { < if (m.op(i) != 0) { < m.let_op(i) = -m.op(i); < res *= -1; < } < } < } < < // x -> 1/x < if (cln::abs(x) >= 2.0) { < map_trafo_H_1overx trafo; < res *= trafo(H(m, xtemp)); < if (cln::imagpart(x) <= 0) { < res = res.subs(H_polesign == -I*Pi); < } else { < res = res.subs(H_polesign == I*Pi); < } < return res.subs(xtemp == numeric(x)).evalf(); < } < < // check transformations for 0.95 <= |x| < 2.0 < < // |(1-x)/(1+x)| < 0.9 -> circular area with center=9,53+0i and radius=9.47 < if (cln::abs(x-9.53) <= 9.47) { < // x -> (1-x)/(1+x) < map_trafo_H_1mxt1px trafo; < res *= trafo(H(m, xtemp)); < } else { < // x -> 1-x < if (has_minus_one) { < map_trafo_H_convert_to_Li filter; < return filter(H(m, numeric(x)).hold()).evalf(); < } < map_trafo_H_1mx trafo; < res *= trafo(H(m, xtemp)); < } < < return res.subs(xtemp == numeric(x)).evalf(); < } < < return H(x1,x2).hold(); < } < < < static ex H_eval(const ex& m_, const ex& x) < { < lst m; < if (is_a<lst>(m_)) { < m = ex_to<lst>(m_); < } else { < m = lst(m_); < } < if (m.nops() == 0) { < return _ex1; < } < ex pos1; < ex pos2; < ex n; < ex p; < int step = 0; < if (*m.begin() > _ex1) { < step++; < pos1 = _ex0; < pos2 = _ex1; < n = *m.begin()-1; < p = _ex1; < } else if (*m.begin() < _ex_1) { < step++; < pos1 = _ex0; < pos2 = _ex_1; < n = -*m.begin()-1; < p = _ex1; < } else if (*m.begin() == _ex0) { < pos1 = _ex0; < n = _ex1; < } else { < pos1 = *m.begin(); < p = _ex1; < } < for (lst::const_iterator it = ++m.begin(); it != m.end(); it++) { < if ((*it).info(info_flags::integer)) { < if (step == 0) { < if (*it > _ex1) { < if (pos1 == _ex0) { < step = 1; < pos2 = _ex1; < n += *it-1; < p = _ex1; < } else { < step = 2; < } < } else if (*it < _ex_1) { < if (pos1 == _ex0) { < step = 1; < pos2 = _ex_1; < n += -*it-1; < p = _ex1; < } else { < step = 2; < } < } else { < if (*it != pos1) { < step = 1; < pos2 = *it; < } < if (*it == _ex0) { < n++; < } else { < p++; < } < } < } else if (step == 1) { < if (*it != pos2) { < step = 2; < } else { < if (*it == _ex0) { < n++; < } else { < p++; < } < } < } < } else { < // if some m_i is not an integer < return H(m_, x).hold(); < } < } < if ((x == _ex1) && (*(--m.end()) != _ex0)) { < return convert_H_to_zeta(m); < } < if (step == 0) { < if (pos1 == _ex0) { < // all zero < if (x == _ex0) { < return H(m_, x).hold(); < } < return pow(log(x), m.nops()) / factorial(m.nops()); < } else { < // all (minus) one < return pow(-pos1*log(1-pos1*x), m.nops()) / factorial(m.nops()); < } < } else if ((step == 1) && (pos1 == _ex0)){ < // convertible to S < if (pos2 == _ex1) { < return S(n, p, x); < } else { < return pow(-1, p) * S(n, p, -x); < } < } < if (x == _ex0) { < return _ex0; < } < if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) { < return H(m_, x).evalf(); < } < return H(m_, x).hold(); < } < < < 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); < } < < < static ex H_deriv(const ex& m_, const ex& x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param < 2); < if (deriv_param == 0) { < return _ex0; < } < lst m; < if (is_a<lst>(m_)) { < m = ex_to<lst>(m_); < } else { < m = lst(m_); < } < ex mb = *m.begin(); < if (mb > _ex1) { < m[0]--; < return H(m, x) / x; < } < if (mb < _ex_1) { < m[0]++; < return H(m, x) / x; < } < m.remove_first(); < if (mb == _ex1) { < return 1/(1-x) * H(m, x); < } else if (mb == _ex_1) { < return 1/(1+x) * H(m, x); < } else { < return H(m, x) / x; < } < } < < < static void H_print_latex(const ex& m_, const ex& x, const print_context& c) < { < lst m; < if (is_a<lst>(m_)) { < m = ex_to<lst>(m_); < } else { < m = lst(m_); < } < c.s << "\\mbox{H}_{"; < lst::const_iterator itm = m.begin(); < (*itm).print(c); < itm++; < for (; itm != m.end(); itm++) { < c.s << ","; < (*itm).print(c); < } < c.s << "}("; < x.print(c); < c.s << ")"; < } < < < REGISTER_FUNCTION(H, < evalf_func(H_evalf). < eval_func(H_eval). < series_func(H_series). < derivative_func(H_deriv). < print_func<print_latex>(H_print_latex). < do_not_evalf_params()); < < < // takes a parameter list for H and returns an expression with corresponding multiple polylogarithms < ex convert_H_to_Li(const ex& m, const ex& x) < { < map_trafo_H_reduce_trailing_zeros filter; < map_trafo_H_convert_to_Li filter2; < if (is_a<lst>(m)) { < return filter2(filter(H(m, x).hold())); < } else { < return filter2(filter(H(lst(m), x).hold())); < } < } < < < ////////////////////////////////////////////////////////////////////// < // < // Multiple zeta values zeta(x) and zeta(x,s) < // < // helper functions < // < ////////////////////////////////////////////////////////////////////// < < < // anonymous namespace for helper functions < namespace { < < < // parameters and data for [Cra] algorithm < const cln::cl_N lambda = cln::cl_N("319/320"); < int L1; < int L2; < std::vector<std::vector<cln::cl_N> > f_kj; < std::vector<cln::cl_N> crB; < std::vector<std::vector<cln::cl_N> > crG; < std::vector<cln::cl_N> crX; < < < void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c) < { < const int size = a.size(); < for (int n=0; n<size; n++) { < c[n] = 0; < for (int m=0; m<=n; m++) { < c[n] = c[n] + a[m]*b[n-m]; < } < } < } < < < // [Cra] section 4 < void initcX(const std::vector<int>& s) < { < const int k = s.size(); < < crX.clear(); < crG.clear(); < crB.clear(); < < for (int i=0; i<=L2; i++) { < crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i)); < } < < int Sm = 0; < int Smp1 = 0; < for (int m=0; m<k-1; m++) { < std::vector<cln::cl_N> crGbuf; < Sm = Sm + s[m]; < Smp1 = Sm + s[m+1]; < for (int i=0; i<=L2; i++) { < crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2)); < } < crG.push_back(crGbuf); < } < < crX = crB; < < for (int m=0; m<k-1; m++) { < std::vector<cln::cl_N> Xbuf; < for (int i=0; i<=L2; i++) { < Xbuf.push_back(crX[i] * crG[m][i]); < } < halfcyclic_convolute(Xbuf, crB, crX); < } < } < < < // [Cra] section 4 < cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk) < { < cln::cl_F one = cln::cl_float(1, cln::float_format(Digits)); < cln::cl_N factor = cln::expt(lambda, Sqk); < cln::cl_N res = factor / Sqk * crX[0] * one; < cln::cl_N resbuf; < int N = 0; < do { < resbuf = res; < factor = factor * lambda; < N++; < res = res + crX[N] * factor / (N+Sqk); < } while ((res != resbuf) || cln::zerop(crX[N])); < return res; < } < < < // [Cra] section 4 < void calc_f(int maxr) < { < f_kj.clear(); < f_kj.resize(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(); < cln::cl_F one = cln::cl_float(1, cln::float_format(Digits)); < < t0 = cln::exp(-lambda); < t2 = 1; < for (k=1; k<=L1; k++) { < t1 = k * lambda; < t2 = t0 * t2; < for (j=1; j<=maxr; j++) { < t3 = 1; < t4 = 1; < for (i=2; i<=j; i++) { < t4 = t4 * (j-i+1); < t3 = t1 * t3 + t4; < } < (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one); < } < it++; < } < } < < < // [Cra] (3.1) < cln::cl_N crandall_Z(const std::vector<int>& s) < { < const int j = s.size(); < < if (j == 1) { < cln::cl_N t0; < cln::cl_N t0buf; < int q = 0; < do { < t0buf = t0; < q++; < t0 = t0 + f_kj[q+j-2][s[0]-1]; < } while (t0 != t0buf); < < return t0 / cln::factorial(s[0]-1); < } < < std::vector<cln::cl_N> t(j); < < 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),s[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), s[k]); < } < t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1]; < } while (t[0] != t0buf); < < return t[0] / cln::factorial(s[0]-1); < } < < < // [Cra] (2.4) < cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s) < { < std::vector<int> r = s; < const int j = r.size(); < < // decide on maximal size of f_kj for crandall_Z < if (Digits < 50) { < L1 = 150; < } else { < L1 = Digits * 3 + j*2; < } < < // decide on maximal size of crX for crandall_Y < if (Digits < 38) { < L2 = 63; < } else if (Digits < 86) { < L2 = 127; < } else if (Digits < 192) { < L2 = 255; < } else if (Digits < 394) { < L2 = 511; < } else if (Digits < 808) { < L2 = 1023; < } else { < L2 = 2047; < } < < cln::cl_N res; < < int maxr = 0; < int S = 0; < for (int i=0; i<j; i++) { < S += r[i]; < if (r[i] > maxr) { < maxr = r[i]; < } < } < < calc_f(maxr); < < const cln::cl_N r0factorial = cln::factorial(r[0]-1); < < std::vector<int> rz; < int skp1buf; < int Srun = S; < for (int k=r.size()-1; k>0; k--) { < < rz.insert(rz.begin(), r.back()); < skp1buf = rz.front(); < Srun -= skp1buf; < r.pop_back(); < < initcX(r); < < for (int q=0; q<skp1buf; q++) { < < cln::cl_N pp1 = crandall_Y_loop(Srun+q-k); < cln::cl_N pp2 = crandall_Z(rz); < < rz.front()--; < < if (q & 1) { < res = res - pp1 * pp2 / cln::factorial(q); < } else { < res = res + pp1 * pp2 / cln::factorial(q); < } < } < rz.front() = skp1buf; < } < rz.insert(rz.begin(), r.back()); < < initcX(rz); < < res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz); < < return res; < } < < < cln::cl_N zeta_do_sum_simple(const std::vector<int>& r) < { < const int j = r.size(); < < // buffer for subsums < std::vector<cln::cl_N> t(j); < cln::cl_F one = cln::cl_float(1, cln::float_format(Digits)); < < cln::cl_N t0buf; < int q = 0; < do { < t0buf = t[0]; < q++; < t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]); < for (int k=j-2; k>=0; k--) { < t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]); < } < } while (t[0] != t0buf); < < return t[0]; < } < < < // does Hoelder convolution. see [BBB] (7.0) < cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vector<int>& s_) < { < // prepare parameters < // holds Li arguments in [BBB] notation < std::vector<int> s = s_; < std::vector<int> m_p = m_; < std::vector<int> m_q; < // holds Li arguments in nested sums notation < std::vector<cln::cl_N> s_p(s.size(), cln::cl_N(1)); < s_p[0] = s_p[0] * cln::cl_N("1/2"); < // convert notations < int sig = 1; < for (int i=0; i<s_.size(); i++) { < if (s_[i] < 0) { < sig = -sig; < s_p[i] = -s_p[i]; < } < s[i] = sig * std::abs(s[i]); < } < std::vector<cln::cl_N> s_q; < cln::cl_N signum = 1; < < // first term < cln::cl_N res = multipleLi_do_sum(m_p, s_p); < < // middle terms < do { < < // change parameters < if (s.front() > 0) { < if (m_p.front() == 1) { < m_p.erase(m_p.begin()); < s_p.erase(s_p.begin()); < if (s_p.size() > 0) { < s_p.front() = s_p.front() * cln::cl_N("1/2"); < } < s.erase(s.begin()); < m_q.front()++; < } else { < m_p.front()--; < m_q.insert(m_q.begin(), 1); < if (s_q.size() > 0) { < s_q.front() = s_q.front() * 2; < } < s_q.insert(s_q.begin(), cln::cl_N("1/2")); < } < } else { < if (m_p.front() == 1) { < m_p.erase(m_p.begin()); < cln::cl_N spbuf = s_p.front(); < s_p.erase(s_p.begin()); < if (s_p.size() > 0) { < s_p.front() = s_p.front() * spbuf; < } < s.erase(s.begin()); < m_q.insert(m_q.begin(), 1); < if (s_q.size() > 0) { < s_q.front() = s_q.front() * 4; < } < s_q.insert(s_q.begin(), cln::cl_N("1/4")); < signum = -signum; < } else { < m_p.front()--; < m_q.insert(m_q.begin(), 1); < if (s_q.size() > 0) { < s_q.front() = s_q.front() * 2; < } < s_q.insert(s_q.begin(), cln::cl_N("1/2")); < } < } < < // exiting the loop < if (m_p.size() == 0) break; < < res = res + signum * multipleLi_do_sum(m_p, s_p) * multipleLi_do_sum(m_q, s_q); < < } while (true); < < // last term < res = res + signum * multipleLi_do_sum(m_q, s_q); < < return res; < } < < < } // end of anonymous namespace < < < ////////////////////////////////////////////////////////////////////// < // < // Multiple zeta values zeta(x) < // < // GiNaC function < // < ////////////////////////////////////////////////////////////////////// < < < static ex zeta1_evalf(const ex& x) < { < if (is_exactly_a<lst>(x) && (x.nops()>1)) { < < // multiple zeta value < const int count = x.nops(); < const lst& xlst = ex_to<lst>(x); < std::vector<int> r(count); < < // check parameters and convert them < lst::const_iterator it1 = xlst.begin(); < std::vector<int>::iterator it2 = r.begin(); < do { < if (!(*it1).info(info_flags::posint)) { < return zeta(x).hold(); < } < *it2 = ex_to<numeric>(*it1).to_int(); < it1++; < it2++; < } while (it2 != r.end()); < < // check for divergence < if (r[0] == 1) { < return zeta(x).hold(); < } < < // decide on summation algorithm < // this is still a bit clumsy < int limit = (Digits>17) ? 10 : 6; < if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) { < return numeric(zeta_do_sum_Crandall(r)); < } else { < return numeric(zeta_do_sum_simple(r)); < } < } < < // single zeta value < if (is_exactly_a<numeric>(x) && (x != 1)) { < try { < return zeta(ex_to<numeric>(x)); < } catch (const dunno &e) { } < } < < return zeta(x).hold(); < } < < < static ex zeta1_eval(const ex& m) < { < if (is_exactly_a<lst>(m)) { < if (m.nops() == 1) { < return zeta(m.op(0)); < } < return zeta(m).hold(); < } < < if (m.info(info_flags::numeric)) { < const numeric& y = ex_to<numeric>(m); < // trap integer arguments: < if (y.is_integer()) { < if (y.is_zero()) { < return _ex_1_2; < } < if (y.is_equal(*_num1_p)) { < return zeta(m).hold(); < } < if (y.info(info_flags::posint)) { < if (y.info(info_flags::odd)) { < return zeta(m).hold(); < } else { < return abs(bernoulli(y)) * pow(Pi, y) * pow(*_num2_p, y-(*_num1_p)) / factorial(y); < } < } else { < if (y.info(info_flags::odd)) { < return -bernoulli((*_num1_p)-y) / ((*_num1_p)-y); < } else { < return _ex0; < } < } < } < // zeta(float) < if (y.info(info_flags::numeric) && !y.info(info_flags::crational)) { < return zeta1_evalf(m); < } < } < return zeta(m).hold(); < } < < < static ex zeta1_deriv(const ex& m, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < if (is_exactly_a<lst>(m)) { < return _ex0; < } else { < return zetaderiv(_ex1, m); < } < } < < < static void zeta1_print_latex(const ex& m_, const print_context& c) < { < c.s << "\\zeta("; < if (is_a<lst>(m_)) { < const lst& m = ex_to<lst>(m_); < lst::const_iterator it = m.begin(); < (*it).print(c); < it++; < for (; it != m.end(); it++) { < c.s << ","; < (*it).print(c); < } < } else { < m_.print(c); < } < c.s << ")"; < } < < < unsigned zeta1_SERIAL::serial = function::register_new(function_options("zeta", 1). < evalf_func(zeta1_evalf). < eval_func(zeta1_eval). < derivative_func(zeta1_deriv). < print_func<print_latex>(zeta1_print_latex). < do_not_evalf_params(). < overloaded(2)); < < < ////////////////////////////////////////////////////////////////////// < // < // Alternating Euler sum zeta(x,s) < // < // GiNaC function < // < ////////////////////////////////////////////////////////////////////// < < < static ex zeta2_evalf(const ex& x, const ex& s) < { < if (is_exactly_a<lst>(x)) { < < // alternating Euler sum < const int count = x.nops(); < const lst& xlst = ex_to<lst>(x); < const lst& slst = ex_to<lst>(s); < std::vector<int> xi(count); < std::vector<int> si(count); < < // check parameters and convert them < lst::const_iterator it_xread = xlst.begin(); < lst::const_iterator it_sread = slst.begin(); < std::vector<int>::iterator it_xwrite = xi.begin(); < std::vector<int>::iterator it_swrite = si.begin(); < do { < if (!(*it_xread).info(info_flags::posint)) { < return zeta(x, s).hold(); < } < *it_xwrite = ex_to<numeric>(*it_xread).to_int(); < if (*it_sread > 0) { < *it_swrite = 1; < } else { < *it_swrite = -1; < } < it_xread++; < it_sread++; < it_xwrite++; < it_swrite++; < } while (it_xwrite != xi.end()); < < // check for divergence < if ((xi[0] == 1) && (si[0] == 1)) { < return zeta(x, s).hold(); < } < < // use Hoelder convolution < return numeric(zeta_do_Hoelder_convolution(xi, si)); < } < < return zeta(x, s).hold(); < } < < < static ex zeta2_eval(const ex& m, const ex& s_) < { < if (is_exactly_a<lst>(s_)) { < const lst& s = ex_to<lst>(s_); < for (lst::const_iterator it = s.begin(); it != s.end(); it++) { < if ((*it).info(info_flags::positive)) { < continue; < } < return zeta(m, s_).hold(); < } < return zeta(m); < } else if (s_.info(info_flags::positive)) { < return zeta(m); < } < < return zeta(m, s_).hold(); < } < < < static ex zeta2_deriv(const ex& m, const ex& s, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < if (is_exactly_a<lst>(m)) { < return _ex0; < } else { < if ((is_exactly_a<lst>(s) && s.op(0).info(info_flags::positive)) || s.info(info_flags::positive)) { < return zetaderiv(_ex1, m); < } < return _ex0; < } < } < < < static void zeta2_print_latex(const ex& m_, const ex& s_, const print_context& c) < { < lst m; < if (is_a<lst>(m_)) { < m = ex_to<lst>(m_); < } else { < m = lst(m_); < } < lst s; < if (is_a<lst>(s_)) { < s = ex_to<lst>(s_); < } else { < s = lst(s_); < } < c.s << "\\zeta("; < lst::const_iterator itm = m.begin(); < lst::const_iterator its = s.begin(); < if (*its < 0) { < c.s << "\\overline{"; < (*itm).print(c); < c.s << "}"; < } else { < (*itm).print(c); < } < its++; < itm++; < for (; itm != m.end(); itm++, its++) { < c.s << ","; < if (*its < 0) { < c.s << "\\overline{"; < (*itm).print(c); < c.s << "}"; < } else { < (*itm).print(c); < } < } < c.s << ")"; < } < < < unsigned zeta2_SERIAL::serial = function::register_new(function_options("zeta", 2). < evalf_func(zeta2_evalf). < eval_func(zeta2_eval). < derivative_func(zeta2_deriv). < print_func<print_latex>(zeta2_print_latex). < do_not_evalf_params(). < overloaded(2)); < < < } // namespace GiNaC < Index: ginac/inifcns_trans.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/inifcns_trans.cpp,v retrieving revision 1.62 diff -r1.62 inifcns_trans.cpp 1,1240d0 < /** @file inifcns_trans.cpp < * < * Implementation of transcendental (and trigonometric and hyperbolic) < * functions. */ < < /* < * GiNaC Copyright (C) 1999-2005 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #include <vector> < #include <stdexcept> < < #include "inifcns.h" < #include "ex.h" < #include "constant.h" < #include "numeric.h" < #include "power.h" < #include "operators.h" < #include "relational.h" < #include "symbol.h" < #include "pseries.h" < #include "utils.h" < < namespace GiNaC { < < ////////// < // exponential function < ////////// < < static ex exp_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return exp(ex_to<numeric>(x)); < < return exp(x).hold(); < } < < static ex exp_eval(const ex & x) < { < // exp(0) -> 1 < if (x.is_zero()) { < return _ex1; < } < < // exp(n*Pi*I/2) -> {+1|+I|-1|-I} < const ex TwoExOverPiI=(_ex2*x)/(Pi*I); < if (TwoExOverPiI.info(info_flags::integer)) { < const numeric z = mod(ex_to<numeric>(TwoExOverPiI),*_num4_p); < if (z.is_equal(*_num0_p)) < return _ex1; < if (z.is_equal(*_num1_p)) < return ex(I); < if (z.is_equal(*_num2_p)) < return _ex_1; < if (z.is_equal(*_num3_p)) < return ex(-I); < } < < // exp(log(x)) -> x < if (is_ex_the_function(x, log)) < return x.op(0); < < // exp(float) -> float < if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) < return exp(ex_to<numeric>(x)); < < return exp(x).hold(); < } < < static ex exp_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx exp(x) -> exp(x) < return exp(x); < } < < REGISTER_FUNCTION(exp, eval_func(exp_eval). < evalf_func(exp_evalf). < derivative_func(exp_deriv). < latex_name("\\exp")); < < ////////// < // natural logarithm < ////////// < < static ex log_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return log(ex_to<numeric>(x)); < < return log(x).hold(); < } < < static ex log_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < if (x.is_zero()) // log(0) -> infinity < throw(pole_error("log_eval(): log(0)",0)); < if (x.info(info_flags::rational) && x.info(info_flags::negative)) < return (log(-x)+I*Pi); < if (x.is_equal(_ex1)) // log(1) -> 0 < return _ex0; < if (x.is_equal(I)) // log(I) -> Pi*I/2 < return (Pi*I*_ex1_2); < if (x.is_equal(-I)) // log(-I) -> -Pi*I/2 < return (Pi*I*_ex_1_2); < < // log(float) -> float < if (!x.info(info_flags::crational)) < return log(ex_to<numeric>(x)); < } < < // log(exp(t)) -> t (if -Pi < t.imag() <= Pi): < if (is_ex_the_function(x, exp)) { < const ex &t = x.op(0); < if (t.info(info_flags::real)) < return t; < } < < return log(x).hold(); < } < < static ex log_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx log(x) -> 1/x < return power(x, _ex_1); < } < < static ex log_series(const ex &arg, < const relational &rel, < int order, < unsigned options) < { < GINAC_ASSERT(is_a<symbol>(rel.lhs())); < ex arg_pt; < bool must_expand_arg = false; < // maybe substitution of rel into arg fails because of a pole < try { < arg_pt = arg.subs(rel, subs_options::no_pattern); < } catch (pole_error) { < must_expand_arg = true; < } < // or we are at the branch point anyways < if (arg_pt.is_zero()) < must_expand_arg = true; < < if (must_expand_arg) { < // method: < // This is the branch point: Series expand the argument first, then < // trivially factorize it to isolate that part which has constant < // leading coefficient in this fashion: < // x^n + x^(n+1) +...+ Order(x^(n+m)) -> x^n * (1 + x +...+ Order(x^m)). < // Return a plain n*log(x) for the x^n part and series expand the < // other part. Add them together and reexpand again in order to have < // one unnested pseries object. All this also works for negative n. < pseries argser; // series expansion of log's argument < unsigned extra_ord = 0; // extra expansion order < do { < // oops, the argument expanded to a pure Order(x^something)... < argser = ex_to<pseries>(arg.series(rel, order+extra_ord, options)); < ++extra_ord; < } while (!argser.is_terminating() && argser.nops()==1); < < const symbol &s = ex_to<symbol>(rel.lhs()); < const ex &point = rel.rhs(); < const int n = argser.ldegree(s); < epvector seq; < // construct what we carelessly called the n*log(x) term above < const ex coeff = argser.coeff(s, n); < // expand the log, but only if coeff is real and > 0, since otherwise < // it would make the branch cut run into the wrong direction < if (coeff.info(info_flags::positive)) < seq.push_back(expair(n*log(s-point)+log(coeff), _ex0)); < else < seq.push_back(expair(log(coeff*pow(s-point, n)), _ex0)); < < if (!argser.is_terminating() || argser.nops()!=1) { < // in this case n more (or less) terms are needed < // (sadly, to generate them, we have to start from the beginning) < if (n == 0 && coeff == 1) { < epvector epv; < ex acc = (new pseries(rel, epv))->setflag(status_flags::dynallocated); < epv.reserve(2); < epv.push_back(expair(-1, _ex0)); < epv.push_back(expair(Order(_ex1), order)); < ex rest = pseries(rel, epv).add_series(argser); < for (int i = order-1; i>0; --i) { < epvector cterm; < cterm.reserve(1); < cterm.push_back(expair(i%2 ? _ex1/i : _ex_1/i, _ex0)); < acc = pseries(rel, cterm).add_series(ex_to<pseries>(acc)); < acc = (ex_to<pseries>(rest)).mul_series(ex_to<pseries>(acc)); < } < return acc; < } < const ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true); < return pseries(rel, seq).add_series(ex_to<pseries>(log(newarg).series(rel, order, options))); < } else // it was a monomial < return pseries(rel, seq); < } < if (!(options & series_options::suppress_branchcut) && < arg_pt.info(info_flags::negative)) { < // method: < // This is the branch cut: assemble the primitive series manually and < // then add the corresponding complex step function. < const symbol &s = ex_to<symbol>(rel.lhs()); < const ex &point = rel.rhs(); < const symbol foo; < const ex replarg = series(log(arg), s==foo, order).subs(foo==point, subs_options::no_pattern); < epvector seq; < seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0)); < seq.push_back(expair(Order(_ex1), order)); < return series(replarg - I*Pi + pseries(rel, seq), rel, order); < } < throw do_taylor(); // caught by function::series() < } < < REGISTER_FUNCTION(log, eval_func(log_eval). < evalf_func(log_evalf). < derivative_func(log_deriv). < series_func(log_series). < latex_name("\\ln")); < < ////////// < // sine (trigonometric function) < ////////// < < static ex sin_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return sin(ex_to<numeric>(x)); < < return sin(x).hold(); < } < < static ex sin_eval(const ex & x) < { < // sin(n/d*Pi) -> { all known non-nested radicals } < const ex SixtyExOverPi = _ex60*x/Pi; < ex sign = _ex1; < if (SixtyExOverPi.info(info_flags::integer)) { < numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p); < if (z>=*_num60_p) { < // wrap to interval [0, Pi) < z -= *_num60_p; < sign = _ex_1; < } < if (z>*_num30_p) { < // wrap to interval [0, Pi/2) < z = *_num60_p-z; < } < if (z.is_equal(*_num0_p)) // sin(0) -> 0 < return _ex0; < if (z.is_equal(*_num5_p)) // sin(Pi/12) -> sqrt(6)/4*(1-sqrt(3)/3) < return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3)); < if (z.is_equal(*_num6_p)) // sin(Pi/10) -> sqrt(5)/4-1/4 < return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4); < if (z.is_equal(*_num10_p)) // sin(Pi/6) -> 1/2 < return sign*_ex1_2; < if (z.is_equal(*_num15_p)) // sin(Pi/4) -> sqrt(2)/2 < return sign*_ex1_2*sqrt(_ex2); < if (z.is_equal(*_num18_p)) // sin(3/10*Pi) -> sqrt(5)/4+1/4 < return sign*(_ex1_4*sqrt(_ex5)+_ex1_4); < if (z.is_equal(*_num20_p)) // sin(Pi/3) -> sqrt(3)/2 < return sign*_ex1_2*sqrt(_ex3); < if (z.is_equal(*_num25_p)) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3) < return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3)); < if (z.is_equal(*_num30_p)) // sin(Pi/2) -> 1 < return sign; < } < < if (is_exactly_a<function>(x)) { < const ex &t = x.op(0); < < // sin(asin(x)) -> x < if (is_ex_the_function(x, asin)) < return t; < < // sin(acos(x)) -> sqrt(1-x^2) < if (is_ex_the_function(x, acos)) < return sqrt(_ex1-power(t,_ex2)); < < // sin(atan(x)) -> x/sqrt(1+x^2) < if (is_ex_the_function(x, atan)) < return t*power(_ex1+power(t,_ex2),_ex_1_2); < } < < // sin(float) -> float < if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) < return sin(ex_to<numeric>(x)); < < // sin() is odd < if (x.info(info_flags::negative)) < return -sin(-x); < < return sin(x).hold(); < } < < static ex sin_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx sin(x) -> cos(x) < return cos(x); < } < < REGISTER_FUNCTION(sin, eval_func(sin_eval). < evalf_func(sin_evalf). < derivative_func(sin_deriv). < latex_name("\\sin")); < < ////////// < // cosine (trigonometric function) < ////////// < < static ex cos_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return cos(ex_to<numeric>(x)); < < return cos(x).hold(); < } < < static ex cos_eval(const ex & x) < { < // cos(n/d*Pi) -> { all known non-nested radicals } < const ex SixtyExOverPi = _ex60*x/Pi; < ex sign = _ex1; < if (SixtyExOverPi.info(info_flags::integer)) { < numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p); < if (z>=*_num60_p) { < // wrap to interval [0, Pi) < z = *_num120_p-z; < } < if (z>=*_num30_p) { < // wrap to interval [0, Pi/2) < z = *_num60_p-z; < sign = _ex_1; < } < if (z.is_equal(*_num0_p)) // cos(0) -> 1 < return sign; < if (z.is_equal(*_num5_p)) // cos(Pi/12) -> sqrt(6)/4*(1+sqrt(3)/3) < return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3)); < if (z.is_equal(*_num10_p)) // cos(Pi/6) -> sqrt(3)/2 < return sign*_ex1_2*sqrt(_ex3); < if (z.is_equal(*_num12_p)) // cos(Pi/5) -> sqrt(5)/4+1/4 < return sign*(_ex1_4*sqrt(_ex5)+_ex1_4); < if (z.is_equal(*_num15_p)) // cos(Pi/4) -> sqrt(2)/2 < return sign*_ex1_2*sqrt(_ex2); < if (z.is_equal(*_num20_p)) // cos(Pi/3) -> 1/2 < return sign*_ex1_2; < if (z.is_equal(*_num24_p)) // cos(2/5*Pi) -> sqrt(5)/4-1/4x < return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4); < if (z.is_equal(*_num25_p)) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3) < return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3)); < if (z.is_equal(*_num30_p)) // cos(Pi/2) -> 0 < return _ex0; < } < < if (is_exactly_a<function>(x)) { < const ex &t = x.op(0); < < // cos(acos(x)) -> x < if (is_ex_the_function(x, acos)) < return t; < < // cos(asin(x)) -> sqrt(1-x^2) < if (is_ex_the_function(x, asin)) < return sqrt(_ex1-power(t,_ex2)); < < // cos(atan(x)) -> 1/sqrt(1+x^2) < if (is_ex_the_function(x, atan)) < return power(_ex1+power(t,_ex2),_ex_1_2); < } < < // cos(float) -> float < if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) < return cos(ex_to<numeric>(x)); < < // cos() is even < if (x.info(info_flags::negative)) < return cos(-x); < < return cos(x).hold(); < } < < static ex cos_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx cos(x) -> -sin(x) < return -sin(x); < } < < REGISTER_FUNCTION(cos, eval_func(cos_eval). < evalf_func(cos_evalf). < derivative_func(cos_deriv). < latex_name("\\cos")); < < ////////// < // tangent (trigonometric function) < ////////// < < static ex tan_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return tan(ex_to<numeric>(x)); < < return tan(x).hold(); < } < < static ex tan_eval(const ex & x) < { < // tan(n/d*Pi) -> { all known non-nested radicals } < const ex SixtyExOverPi = _ex60*x/Pi; < ex sign = _ex1; < if (SixtyExOverPi.info(info_flags::integer)) { < numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num60_p); < if (z>=*_num60_p) { < // wrap to interval [0, Pi) < z -= *_num60_p; < } < if (z>=*_num30_p) { < // wrap to interval [0, Pi/2) < z = *_num60_p-z; < sign = _ex_1; < } < if (z.is_equal(*_num0_p)) // tan(0) -> 0 < return _ex0; < if (z.is_equal(*_num5_p)) // tan(Pi/12) -> 2-sqrt(3) < return sign*(_ex2-sqrt(_ex3)); < if (z.is_equal(*_num10_p)) // tan(Pi/6) -> sqrt(3)/3 < return sign*_ex1_3*sqrt(_ex3); < if (z.is_equal(*_num15_p)) // tan(Pi/4) -> 1 < return sign; < if (z.is_equal(*_num20_p)) // tan(Pi/3) -> sqrt(3) < return sign*sqrt(_ex3); < if (z.is_equal(*_num25_p)) // tan(5/12*Pi) -> 2+sqrt(3) < return sign*(sqrt(_ex3)+_ex2); < if (z.is_equal(*_num30_p)) // tan(Pi/2) -> infinity < throw (pole_error("tan_eval(): simple pole",1)); < } < < if (is_exactly_a<function>(x)) { < const ex &t = x.op(0); < < // tan(atan(x)) -> x < if (is_ex_the_function(x, atan)) < return t; < < // tan(asin(x)) -> x/sqrt(1+x^2) < if (is_ex_the_function(x, asin)) < return t*power(_ex1-power(t,_ex2),_ex_1_2); < < // tan(acos(x)) -> sqrt(1-x^2)/x < if (is_ex_the_function(x, acos)) < return power(t,_ex_1)*sqrt(_ex1-power(t,_ex2)); < } < < // tan(float) -> float < if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) { < return tan(ex_to<numeric>(x)); < } < < // tan() is odd < if (x.info(info_flags::negative)) < return -tan(-x); < < return tan(x).hold(); < } < < static ex tan_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx tan(x) -> 1+tan(x)^2; < return (_ex1+power(tan(x),_ex2)); < } < < static ex tan_series(const ex &x, < const relational &rel, < int order, < unsigned options) < { < GINAC_ASSERT(is_a<symbol>(rel.lhs())); < // method: < // Taylor series where there is no pole falls back to tan_deriv. < // On a pole simply expand sin(x)/cos(x). < const ex x_pt = x.subs(rel, subs_options::no_pattern); < if (!(2*x_pt/Pi).info(info_flags::odd)) < throw do_taylor(); // caught by function::series() < // if we got here we have to care for a simple pole < return (sin(x)/cos(x)).series(rel, order, options); < } < < REGISTER_FUNCTION(tan, eval_func(tan_eval). < evalf_func(tan_evalf). < derivative_func(tan_deriv). < series_func(tan_series). < latex_name("\\tan")); < < ////////// < // inverse sine (arc sine) < ////////// < < static ex asin_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return asin(ex_to<numeric>(x)); < < return asin(x).hold(); < } < < static ex asin_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // asin(0) -> 0 < if (x.is_zero()) < return x; < < // asin(1/2) -> Pi/6 < if (x.is_equal(_ex1_2)) < return numeric(1,6)*Pi; < < // asin(1) -> Pi/2 < if (x.is_equal(_ex1)) < return _ex1_2*Pi; < < // asin(-1/2) -> -Pi/6 < if (x.is_equal(_ex_1_2)) < return numeric(-1,6)*Pi; < < // asin(-1) -> -Pi/2 < if (x.is_equal(_ex_1)) < return _ex_1_2*Pi; < < // asin(float) -> float < if (!x.info(info_flags::crational)) < return asin(ex_to<numeric>(x)); < < // asin() is odd < if (x.info(info_flags::negative)) < return -asin(-x); < } < < return asin(x).hold(); < } < < static ex asin_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx asin(x) -> 1/sqrt(1-x^2) < return power(1-power(x,_ex2),_ex_1_2); < } < < REGISTER_FUNCTION(asin, eval_func(asin_eval). < evalf_func(asin_evalf). < derivative_func(asin_deriv). < latex_name("\\arcsin")); < < ////////// < // inverse cosine (arc cosine) < ////////// < < static ex acos_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return acos(ex_to<numeric>(x)); < < return acos(x).hold(); < } < < static ex acos_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // acos(1) -> 0 < if (x.is_equal(_ex1)) < return _ex0; < < // acos(1/2) -> Pi/3 < if (x.is_equal(_ex1_2)) < return _ex1_3*Pi; < < // acos(0) -> Pi/2 < if (x.is_zero()) < return _ex1_2*Pi; < < // acos(-1/2) -> 2/3*Pi < if (x.is_equal(_ex_1_2)) < return numeric(2,3)*Pi; < < // acos(-1) -> Pi < if (x.is_equal(_ex_1)) < return Pi; < < // acos(float) -> float < if (!x.info(info_flags::crational)) < return acos(ex_to<numeric>(x)); < < // acos(-x) -> Pi-acos(x) < if (x.info(info_flags::negative)) < return Pi-acos(-x); < } < < return acos(x).hold(); < } < < static ex acos_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx acos(x) -> -1/sqrt(1-x^2) < return -power(1-power(x,_ex2),_ex_1_2); < } < < REGISTER_FUNCTION(acos, eval_func(acos_eval). < evalf_func(acos_evalf). < derivative_func(acos_deriv). < latex_name("\\arccos")); < < ////////// < // inverse tangent (arc tangent) < ////////// < < static ex atan_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return atan(ex_to<numeric>(x)); < < return atan(x).hold(); < } < < static ex atan_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // atan(0) -> 0 < if (x.is_zero()) < return _ex0; < < // atan(1) -> Pi/4 < if (x.is_equal(_ex1)) < return _ex1_4*Pi; < < // atan(-1) -> -Pi/4 < if (x.is_equal(_ex_1)) < return _ex_1_4*Pi; < < if (x.is_equal(I) || x.is_equal(-I)) < throw (pole_error("atan_eval(): logarithmic pole",0)); < < // atan(float) -> float < if (!x.info(info_flags::crational)) < return atan(ex_to<numeric>(x)); < < // atan() is odd < if (x.info(info_flags::negative)) < return -atan(-x); < } < < return atan(x).hold(); < } < < static ex atan_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx atan(x) -> 1/(1+x^2) < return power(_ex1+power(x,_ex2), _ex_1); < } < < static ex atan_series(const ex &arg, < const relational &rel, < int order, < unsigned options) < { < GINAC_ASSERT(is_a<symbol>(rel.lhs())); < // method: < // Taylor series where there is no pole or cut falls back to atan_deriv. < // There are two branch cuts, one runnig from I up the imaginary axis and < // one running from -I down the imaginary axis. The points I and -I are < // poles. < // On the branch cuts and the poles series expand < // (log(1+I*x)-log(1-I*x))/(2*I) < // instead. < const ex arg_pt = arg.subs(rel, subs_options::no_pattern); < if (!(I*arg_pt).info(info_flags::real)) < throw do_taylor(); // Re(x) != 0 < if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1) < throw do_taylor(); // Re(x) == 0, but abs(x)<1 < // care for the poles, using the defining formula for atan()... < if (arg_pt.is_equal(I) || arg_pt.is_equal(-I)) < return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options); < if (!(options & series_options::suppress_branchcut)) { < // method: < // This is the branch cut: assemble the primitive series manually and < // then add the corresponding complex step function. < const symbol &s = ex_to<symbol>(rel.lhs()); < const ex &point = rel.rhs(); < const symbol foo; < const ex replarg = series(atan(arg), s==foo, order).subs(foo==point, subs_options::no_pattern); < ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2; < if ((I*arg_pt)<_ex0) < Order0correction += log((I*arg_pt+_ex_1)/(I*arg_pt+_ex1))*I*_ex_1_2; < else < Order0correction += log((I*arg_pt+_ex1)/(I*arg_pt+_ex_1))*I*_ex1_2; < epvector seq; < seq.push_back(expair(Order0correction, _ex0)); < seq.push_back(expair(Order(_ex1), order)); < return series(replarg - pseries(rel, seq), rel, order); < } < throw do_taylor(); < } < < REGISTER_FUNCTION(atan, eval_func(atan_eval). < evalf_func(atan_evalf). < derivative_func(atan_deriv). < series_func(atan_series). < latex_name("\\arctan")); < < ////////// < // inverse tangent (atan2(y,x)) < ////////// < < static ex atan2_evalf(const ex &y, const ex &x) < { < if (is_exactly_a<numeric>(y) && is_exactly_a<numeric>(x)) < return atan(ex_to<numeric>(y), ex_to<numeric>(x)); < < return atan2(y, x).hold(); < } < < static ex atan2_eval(const ex & y, const ex & x) < { < if (y.info(info_flags::numeric) && x.info(info_flags::numeric)) { < < if (y.is_zero()) { < < // atan(0, 0) -> 0 < if (x.is_zero()) < return _ex0; < < // atan(0, x), x real and positive -> 0 < if (x.info(info_flags::positive)) < return _ex0; < < // atan(0, x), x real and negative -> -Pi < if (x.info(info_flags::negative)) < return _ex_1*Pi; < } < < if (x.is_zero()) { < < // atan(y, 0), y real and positive -> Pi/2 < if (y.info(info_flags::positive)) < return _ex1_2*Pi; < < // atan(y, 0), y real and negative -> -Pi/2 < if (y.info(info_flags::negative)) < return _ex_1_2*Pi; < } < < if (y.is_equal(x)) { < < // atan(y, y), y real and positive -> Pi/4 < if (y.info(info_flags::positive)) < return _ex1_4*Pi; < < // atan(y, y), y real and negative -> -3/4*Pi < if (y.info(info_flags::negative)) < return numeric(-3, 4)*Pi; < } < < if (y.is_equal(-x)) { < < // atan(y, -y), y real and positive -> 3*Pi/4 < if (y.info(info_flags::positive)) < return numeric(3, 4)*Pi; < < // atan(y, -y), y real and negative -> -Pi/4 < if (y.info(info_flags::negative)) < return _ex_1_4*Pi; < } < < // atan(float, float) -> float < if (!y.info(info_flags::crational) && !x.info(info_flags::crational)) < return atan(ex_to<numeric>(y), ex_to<numeric>(x)); < < // atan(real, real) -> atan(y/x) +/- Pi < if (y.info(info_flags::real) && x.info(info_flags::real)) { < if (x.info(info_flags::positive)) < return atan(y/x); < else if(y.info(info_flags::positive)) < return atan(y/x)+Pi; < else < return atan(y/x)-Pi; < } < } < < return atan2(y, x).hold(); < } < < static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param<2); < < if (deriv_param==0) { < // d/dy atan(y,x) < return x*power(power(x,_ex2)+power(y,_ex2),_ex_1); < } < // d/dx atan(y,x) < return -y*power(power(x,_ex2)+power(y,_ex2),_ex_1); < } < < REGISTER_FUNCTION(atan2, eval_func(atan2_eval). < evalf_func(atan2_evalf). < derivative_func(atan2_deriv)); < < ////////// < // hyperbolic sine (trigonometric function) < ////////// < < static ex sinh_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return sinh(ex_to<numeric>(x)); < < return sinh(x).hold(); < } < < static ex sinh_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // sinh(0) -> 0 < if (x.is_zero()) < return _ex0; < < // sinh(float) -> float < if (!x.info(info_flags::crational)) < return sinh(ex_to<numeric>(x)); < < // sinh() is odd < if (x.info(info_flags::negative)) < return -sinh(-x); < } < < if ((x/Pi).info(info_flags::numeric) && < ex_to<numeric>(x/Pi).real().is_zero()) // sinh(I*x) -> I*sin(x) < return I*sin(x/I); < < if (is_exactly_a<function>(x)) { < const ex &t = x.op(0); < < // sinh(asinh(x)) -> x < if (is_ex_the_function(x, asinh)) < return t; < < // sinh(acosh(x)) -> sqrt(x-1) * sqrt(x+1) < if (is_ex_the_function(x, acosh)) < return sqrt(t-_ex1)*sqrt(t+_ex1); < < // sinh(atanh(x)) -> x/sqrt(1-x^2) < if (is_ex_the_function(x, atanh)) < return t*power(_ex1-power(t,_ex2),_ex_1_2); < } < < return sinh(x).hold(); < } < < static ex sinh_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx sinh(x) -> cosh(x) < return cosh(x); < } < < REGISTER_FUNCTION(sinh, eval_func(sinh_eval). < evalf_func(sinh_evalf). < derivative_func(sinh_deriv). < latex_name("\\sinh")); < < ////////// < // hyperbolic cosine (trigonometric function) < ////////// < < static ex cosh_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return cosh(ex_to<numeric>(x)); < < return cosh(x).hold(); < } < < static ex cosh_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // cosh(0) -> 1 < if (x.is_zero()) < return _ex1; < < // cosh(float) -> float < if (!x.info(info_flags::crational)) < return cosh(ex_to<numeric>(x)); < < // cosh() is even < if (x.info(info_flags::negative)) < return cosh(-x); < } < < if ((x/Pi).info(info_flags::numeric) && < ex_to<numeric>(x/Pi).real().is_zero()) // cosh(I*x) -> cos(x) < return cos(x/I); < < if (is_exactly_a<function>(x)) { < const ex &t = x.op(0); < < // cosh(acosh(x)) -> x < if (is_ex_the_function(x, acosh)) < return t; < < // cosh(asinh(x)) -> sqrt(1+x^2) < if (is_ex_the_function(x, asinh)) < return sqrt(_ex1+power(t,_ex2)); < < // cosh(atanh(x)) -> 1/sqrt(1-x^2) < if (is_ex_the_function(x, atanh)) < return power(_ex1-power(t,_ex2),_ex_1_2); < } < < return cosh(x).hold(); < } < < static ex cosh_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx cosh(x) -> sinh(x) < return sinh(x); < } < < REGISTER_FUNCTION(cosh, eval_func(cosh_eval). < evalf_func(cosh_evalf). < derivative_func(cosh_deriv). < latex_name("\\cosh")); < < ////////// < // hyperbolic tangent (trigonometric function) < ////////// < < static ex tanh_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return tanh(ex_to<numeric>(x)); < < return tanh(x).hold(); < } < < static ex tanh_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // tanh(0) -> 0 < if (x.is_zero()) < return _ex0; < < // tanh(float) -> float < if (!x.info(info_flags::crational)) < return tanh(ex_to<numeric>(x)); < < // tanh() is odd < if (x.info(info_flags::negative)) < return -tanh(-x); < } < < if ((x/Pi).info(info_flags::numeric) && < ex_to<numeric>(x/Pi).real().is_zero()) // tanh(I*x) -> I*tan(x); < return I*tan(x/I); < < if (is_exactly_a<function>(x)) { < const ex &t = x.op(0); < < // tanh(atanh(x)) -> x < if (is_ex_the_function(x, atanh)) < return t; < < // tanh(asinh(x)) -> x/sqrt(1+x^2) < if (is_ex_the_function(x, asinh)) < return t*power(_ex1+power(t,_ex2),_ex_1_2); < < // tanh(acosh(x)) -> sqrt(x-1)*sqrt(x+1)/x < if (is_ex_the_function(x, acosh)) < return sqrt(t-_ex1)*sqrt(t+_ex1)*power(t,_ex_1); < } < < return tanh(x).hold(); < } < < static ex tanh_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx tanh(x) -> 1-tanh(x)^2 < return _ex1-power(tanh(x),_ex2); < } < < static ex tanh_series(const ex &x, < const relational &rel, < int order, < unsigned options) < { < GINAC_ASSERT(is_a<symbol>(rel.lhs())); < // method: < // Taylor series where there is no pole falls back to tanh_deriv. < // On a pole simply expand sinh(x)/cosh(x). < const ex x_pt = x.subs(rel, subs_options::no_pattern); < if (!(2*I*x_pt/Pi).info(info_flags::odd)) < throw do_taylor(); // caught by function::series() < // if we got here we have to care for a simple pole < return (sinh(x)/cosh(x)).series(rel, order, options); < } < < REGISTER_FUNCTION(tanh, eval_func(tanh_eval). < evalf_func(tanh_evalf). < derivative_func(tanh_deriv). < series_func(tanh_series). < latex_name("\\tanh")); < < ////////// < // inverse hyperbolic sine (trigonometric function) < ////////// < < static ex asinh_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return asinh(ex_to<numeric>(x)); < < return asinh(x).hold(); < } < < static ex asinh_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // asinh(0) -> 0 < if (x.is_zero()) < return _ex0; < < // asinh(float) -> float < if (!x.info(info_flags::crational)) < return asinh(ex_to<numeric>(x)); < < // asinh() is odd < if (x.info(info_flags::negative)) < return -asinh(-x); < } < < return asinh(x).hold(); < } < < static ex asinh_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx asinh(x) -> 1/sqrt(1+x^2) < return power(_ex1+power(x,_ex2),_ex_1_2); < } < < REGISTER_FUNCTION(asinh, eval_func(asinh_eval). < evalf_func(asinh_evalf). < derivative_func(asinh_deriv)); < < ////////// < // inverse hyperbolic cosine (trigonometric function) < ////////// < < static ex acosh_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return acosh(ex_to<numeric>(x)); < < return acosh(x).hold(); < } < < static ex acosh_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // acosh(0) -> Pi*I/2 < if (x.is_zero()) < return Pi*I*numeric(1,2); < < // acosh(1) -> 0 < if (x.is_equal(_ex1)) < return _ex0; < < // acosh(-1) -> Pi*I < if (x.is_equal(_ex_1)) < return Pi*I; < < // acosh(float) -> float < if (!x.info(info_flags::crational)) < return acosh(ex_to<numeric>(x)); < < // acosh(-x) -> Pi*I-acosh(x) < if (x.info(info_flags::negative)) < return Pi*I-acosh(-x); < } < < return acosh(x).hold(); < } < < static ex acosh_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1)) < return power(x+_ex_1,_ex_1_2)*power(x+_ex1,_ex_1_2); < } < < REGISTER_FUNCTION(acosh, eval_func(acosh_eval). < evalf_func(acosh_evalf). < derivative_func(acosh_deriv)); < < ////////// < // inverse hyperbolic tangent (trigonometric function) < ////////// < < static ex atanh_evalf(const ex & x) < { < if (is_exactly_a<numeric>(x)) < return atanh(ex_to<numeric>(x)); < < return atanh(x).hold(); < } < < static ex atanh_eval(const ex & x) < { < if (x.info(info_flags::numeric)) { < < // atanh(0) -> 0 < if (x.is_zero()) < return _ex0; < < // atanh({+|-}1) -> throw < if (x.is_equal(_ex1) || x.is_equal(_ex_1)) < throw (pole_error("atanh_eval(): logarithmic pole",0)); < < // atanh(float) -> float < if (!x.info(info_flags::crational)) < return atanh(ex_to<numeric>(x)); < < // atanh() is odd < if (x.info(info_flags::negative)) < return -atanh(-x); < } < < return atanh(x).hold(); < } < < static ex atanh_deriv(const ex & x, unsigned deriv_param) < { < GINAC_ASSERT(deriv_param==0); < < // d/dx atanh(x) -> 1/(1-x^2) < return power(_ex1-power(x,_ex2),_ex_1); < } < < static ex atanh_series(const ex &arg, < const relational &rel, < int order, < unsigned options) < { < GINAC_ASSERT(is_a<symbol>(rel.lhs())); < // method: < // Taylor series where there is no pole or cut falls back to atanh_deriv. < // There are two branch cuts, one runnig from 1 up the real axis and one < // one running from -1 down the real axis. The points 1 and -1 are poles < // On the branch cuts and the poles series expand < // (log(1+x)-log(1-x))/2 < // instead. < const ex arg_pt = arg.subs(rel, subs_options::no_pattern); < if (!(arg_pt).info(info_flags::real)) < throw do_taylor(); // Im(x) != 0 < if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1) < throw do_taylor(); // Im(x) == 0, but abs(x)<1 < // care for the poles, using the defining formula for atanh()... < if (arg_pt.is_equal(_ex1) || arg_pt.is_equal(_ex_1)) < return ((log(_ex1+arg)-log(_ex1-arg))*_ex1_2).series(rel, order, options); < // ...and the branch cuts (the discontinuity at the cut being just I*Pi) < if (!(options & series_options::suppress_branchcut)) { < // method: < // This is the branch cut: assemble the primitive series manually and < // then add the corresponding complex step function. < const symbol &s = ex_to<symbol>(rel.lhs()); < const ex &point = rel.rhs(); < const symbol foo; < const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point, subs_options::no_pattern); < ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2; < if (arg_pt<_ex0) < Order0correction += log((arg_pt+_ex_1)/(arg_pt+_ex1))*_ex1_2; < else < Order0correction += log((arg_pt+_ex1)/(arg_pt+_ex_1))*_ex_1_2; < epvector seq; < seq.push_back(expair(Order0correction, _ex0)); < seq.push_back(expair(Order(_ex1), order)); < return series(replarg - pseries(rel, seq), rel, order); < } < throw do_taylor(); < } < < REGISTER_FUNCTION(atanh, eval_func(atanh_eval). < evalf_func(atanh_evalf). < derivative_func(atanh_deriv). < series_func(atanh_series)); < < < } // namespace GiNaC Index: ginac/input_parser.yy =================================================================== RCS file: /home/cvs/GiNaC/ginac/input_parser.yy,v retrieving revision 1.22 diff -r1.22 input_parser.yy 7c7 < * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany --- > * GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany 109,110c109,120 < unsigned i = function::find_function(n, $3.nops()); < $$ = function(i, ex_to<exprseq>($3)).eval(1); --- > factory_p p = find_func_factory(n+"_function"); > // exprseq -> exvector > const exprseq& in = ex_to<exprseq>($3); > exvector out; > for (exprseq::const_iterator it = in.begin(); it != in.end(); ++it) { > out.push_back(*it); > } > $$ = (*p)(out); > // $$ = n(ex_to<exprseq>($3)); > //TODO > // unsigned i = function::find_function(n, $3.nops()); > // $$ = function(i, ex_to<exprseq>($3)).eval(1); 128c138 < | exp '!' {$$ = factorial($1);} --- > | exp '!' {$$ = factorial($1);} Index: ginac/integral.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/integral.cpp,v retrieving revision 1.9 diff -r1.9 integral.cpp 29a30 > #include "inifcns_exp.h" Index: ginac/matrix.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/matrix.cpp,v retrieving revision 1.107 diff -r1.107 matrix.cpp 30a31,32 > > #include "relational.h" Index: ginac/normal.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/normal.cpp,v retrieving revision 1.103 diff -r1.103 normal.cpp 29a30 > Index: ginac/power.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/power.cpp,v retrieving revision 1.108 diff -r1.108 power.cpp 36c36 < #include "inifcns.h" // for log() in power::derivative() --- > #include "inifcns_exp.h" // for log() in power::derivative() 394,395c394,395 < if (is_exactly_a<function>(ebasis)) < return ex_to<function>(ebasis).power(eexponent); --- > if (is_a<function>(ebasis)) > return ex_to<function>(ebasis).power_law(eexponent); Index: ginac/pseries.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/pseries.cpp,v retrieving revision 1.87 diff -r1.87 pseries.cpp 133c133 < if (!is_order_function(i->rest)) { --- > if (!is_exactly_a<Order_function>(i->rest)) { 270c270 < if (is_order_function(seq[i].rest)) --- > if (is_exactly_a<Order_function>(seq[i].rest)) 502c502 < if (is_order_function(it->rest)) { --- > if (is_exactly_a<Order_function>(it->rest)) { 515c515 < if (is_order_function(it->rest)) { --- > if (is_exactly_a<Order_function>(it->rest)) { 535c535 < if (is_order_function(it->rest)) { --- > if (is_exactly_a<Order_function>(it->rest)) { 547c547 < return seq.empty() || !is_order_function((seq.end()-1)->rest); --- > return seq.empty() || !is_exactly_a<Order_function>((seq.end()-1)->rest); 682c682 < if (is_order_function((*a).rest)) --- > if (is_exactly_a<Order_function>((*a).rest)) 688c688 < if (is_order_function((*b).rest)) --- > if (is_exactly_a<Order_function>((*b).rest)) 693c693 < if (is_order_function((*a).rest) || is_order_function((*b).rest)) { --- > if (is_exactly_a<Order_function>((*a).rest) || is_exactly_a<Order_function>((*b).rest)) { 750c750 < if (!is_order_function(it->rest)) --- > if (!is_exactly_a<Order_function>(it->rest)) 791c791 < if (is_order_function(coeff(var, a_max))) --- > if (is_exactly_a<Order_function>(coeff(var, a_max))) 793c793 < if (is_order_function(other.coeff(var, b_max))) --- > if (is_exactly_a<Order_function>(other.coeff(var, b_max))) 805c805 < if (!is_order_function(a_coeff) && !is_order_function(b_coeff)) --- > if (!is_exactly_a<Order_function>(a_coeff) && !is_exactly_a<Order_function>(b_coeff)) 940c940 < if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative()) --- > if (seq.size() == 1 && is_exactly_a<Order_function>(seq[0].rest) && p.real().is_negative()) 951c951 < if (is_order_function(c)) { --- > if (is_exactly_a<Order_function>(c)) { 966c966 < if (is_order_function(co[i])) { --- > if (is_exactly_a<Order_function>(co[i])) { 1117c1117 < if (is_order_function(currcoeff)) --- > if (is_exactly_a<Order_function>(currcoeff)) 1133c1133 < if (is_order_function(currcoeff)) --- > if (is_exactly_a<Order_function>(currcoeff)) Index: ginac/registrar.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/registrar.cpp,v retrieving revision 1.13 diff -r1.13 registrar.cpp 40a41,45 > factory_p find_func_factory(const std::string& class_name) > { > return registered_class_info::find(class_name)->options.get_func_factory(); > } > Index: ginac/registrar.h =================================================================== RCS file: /home/cvs/GiNaC/ginac/registrar.h,v retrieving revision 1.26 diff -r1.26 registrar.h 47a48,49 > /** Definitions for the function registration mechanism. */ > typedef ex (*factory_p)(const std::vector<ex>& v); 58a61 > factory_p get_func_factory() const { return function_factory_pointer; } 88a92,97 > registered_class_options & func_factory(factory_p p) > { > function_factory_pointer = p; > return *this; > } > 94a104 > factory_p function_factory_pointer; /**< Pointer to function factory. */ 168a179,181 > /** Find function factory by class name. */ > extern factory_p find_func_factory(const std::string& class_name); > Index: ginac/remember.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/remember.cpp,v retrieving revision 1.14 diff -r1.14 remember.cpp 1,190d0 < /** @file remember.cpp < * < * Implementation of helper classes for using the remember option < * in GiNaC functions */ < < /* < * GiNaC Copyright (C) 1999-2005 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #include <stdexcept> < < #include "function.h" < #include "utils.h" < #include "remember.h" < < namespace GiNaC { < < ////////// < // class remember_table_entry < ////////// < < remember_table_entry::remember_table_entry(function const & f, ex const & r) < : hashvalue(f.gethash()), seq(f.seq), result(r) < { < ++last_access = access_counter; < successful_hits = 0; < } < < bool remember_table_entry::is_equal(function const & f) const < { < GINAC_ASSERT(f.seq.size()==seq.size()); < if (f.gethash()!=hashvalue) return false; < size_t num = seq.size(); < for (size_t i=0; i<num; ++i) < if (!seq[i].is_equal(f.seq[i])) return false; < ++last_access = access_counter; < ++successful_hits; < return true; < } < < unsigned long remember_table_entry::access_counter = 0; < < ////////// < // class remember_table_list < ////////// < < remember_table_list::remember_table_list(unsigned as, unsigned strat) < { < max_assoc_size = as; < remember_strategy = strat; < } < < < void remember_table_list::add_entry(function const & f, ex const & result) < { < if ((max_assoc_size!=0) && < (remember_strategy!=remember_strategies::delete_never) && < (size()>=max_assoc_size)) { < // table is full, we must delete an older entry < GINAC_ASSERT(size()>0); // there must be at least one entry < < switch (remember_strategy) { < case remember_strategies::delete_cyclic: { < // delete oldest entry (first in list) < erase(begin()); < break; < } < case remember_strategies::delete_lru: { < // delete least recently used entry < iterator it = begin(); < iterator lowest_access_it = it; < unsigned long lowest_access = (*it).get_last_access(); < ++it; < while (it!=end()) { < if ((*it).get_last_access()<lowest_access) { < lowest_access = (*it).get_last_access(); < lowest_access_it = it; < } < ++it; < } < erase(lowest_access_it); < break; < } < case remember_strategies::delete_lfu: { < // delete least frequently used entry < iterator it = begin(); < iterator lowest_hits_it = it; < unsigned lowest_hits = (*it).get_successful_hits(); < ++it; < while (it!=end()) { < if ((*it).get_successful_hits()<lowest_hits) { < lowest_hits = (*it).get_successful_hits(); < lowest_hits_it = it; < } < ++it; < } < erase(lowest_hits_it); < break; < } < default: < throw(std::logic_error("remember_table_list::add_entry(): invalid remember_strategy")); < } < GINAC_ASSERT(size()==max_assoc_size-1); < } < push_back(remember_table_entry(f,result)); < } < < bool remember_table_list::lookup_entry(function const & f, ex & result) const < { < const_iterator i = begin(), iend = end(); < while (i != iend) { < if (i->is_equal(f)) { < result = i->get_result(); < return true; < } < ++i; < } < return false; < } < < ////////// < // class remember_table < ////////// < < remember_table::remember_table() < { < table_size=0; < max_assoc_size=0; < remember_strategy=remember_strategies::delete_never; < } < < remember_table::remember_table(unsigned s, unsigned as, unsigned strat) < : max_assoc_size(as), remember_strategy(strat) < { < // we keep max_assoc_size and remember_strategy if we need to clear < // all entries < < // use some power of 2 next to s < table_size = 1 << log2(s); < init_table(); < } < < bool remember_table::lookup_entry(function const & f, ex & result) const < { < unsigned entry = f.gethash() & (table_size-1); < GINAC_ASSERT(entry<size()); < return operator[](entry).lookup_entry(f,result); < } < < void remember_table::add_entry(function const & f, ex const & result) < { < unsigned entry = f.gethash() & (table_size-1); < GINAC_ASSERT(entry<size()); < operator[](entry).add_entry(f,result); < } < < void remember_table::clear_all_entries() < { < clear(); < init_table(); < } < < void remember_table::init_table() < { < reserve(table_size); < for (unsigned i=0; i<table_size; ++i) < push_back(remember_table_list(max_assoc_size,remember_strategy)); < } < < std::vector<remember_table> & remember_table::remember_tables() < { < static std::vector<remember_table> * rt = new std::vector<remember_table>; < return *rt; < } < < } // namespace GiNaC Index: ginac/remember.h =================================================================== RCS file: /home/cvs/GiNaC/ginac/remember.h,v retrieving revision 1.14 diff -r1.14 remember.h 1,101d0 < /** @file remember.h < * < * Interface to helper classes for using the remember option < * in GiNaC functions */ < < /* < * GiNaC Copyright (C) 1999-2005 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 < * the Free Software Foundation; either version 2 of the License, or < * (at your option) any later version. < * < * This program is distributed in the hope that it will be useful, < * but WITHOUT ANY WARRANTY; without even the implied warranty of < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the < * GNU General Public License for more details. < * < * You should have received a copy of the GNU General Public License < * along with this program; if not, write to the Free Software < * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA < */ < < #ifndef __GINAC_REMEMBER_H__ < #define __GINAC_REMEMBER_H__ < < #include <iosfwd> < #include <vector> < #include <list> < < namespace GiNaC { < < class function; < class ex; < < /** A single entry in the remember table of a function. < * Needs to be a friend of class function to access 'seq'. < * 'last_access' and 'successful_hits' are updated at each successful < * 'is_equal'. */ < class remember_table_entry { < public: < remember_table_entry(function const & f, ex const & r); < bool is_equal(function const & f) const; < ex get_result() const { return result; } < unsigned long get_last_access() const { return last_access; } < unsigned long get_successful_hits() const { return successful_hits; }; < < protected: < unsigned hashvalue; < exvector seq; < ex result; < mutable unsigned long last_access; < mutable unsigned successful_hits; < static unsigned long access_counter; < }; < < /** A list of entries in the remember table having some least < * significant bits of the hashvalue in common. */ < class remember_table_list : public std::list<remember_table_entry> { < public: < remember_table_list(unsigned as, unsigned strat); < void add_entry(function const & f, ex const & result); < bool lookup_entry(function const & f, ex & result) const; < protected: < unsigned max_assoc_size; < unsigned remember_strategy; < }; < < /** The remember table is organized like an n-fold associative cache < * in a microprocessor. The table has a width of 's' (which is rounded < * to table_size, some power of 2 near 's', internally) and a depth of 'as' < * (unless you choose that entries are never discarded). The place where < * an entry is stored depends on the hashvalue of the parameters of the < * function (this corresponds to the address of byte to be cached). < * The 'log_2(table_size)' least significant bits of this hashvalue < * give the slot in which the entry will be stored or looked up. < * Each slot can take up to 'as' entries. If a slot is full, an older < * entry is removed by one of the following strategies: < * - oldest entry (the first one in the list) < * - least recently used (the one with the lowest 'last_access') < * - least frequently used (the one with the lowest 'successful_hits') < * or all entries are kept which means that the table grows indefinitely. */ < class remember_table : public std::vector<remember_table_list> { < public: < remember_table(); < remember_table(unsigned s, unsigned as, unsigned strat); < bool lookup_entry(function const & f, ex & result) const; < void add_entry(function const & f, ex const & result); < void clear_all_entries(); < void show_statistics(std::ostream & os, unsigned level) const; < static std::vector<remember_table> & remember_tables(); < protected: < void init_table(); < unsigned table_size; < unsigned max_assoc_size; < unsigned remember_strategy; < }; < < } // namespace GiNaC < < #endif // ndef __GINAC_REMEMBER_H__ Index: ginac/symbol.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/symbol.cpp,v retrieving revision 1.62 diff -r1.62 symbol.cpp 223c223 < return GiNaC::conjugate_function(*this).hold(); --- > return GiNaC::conjugate(*this); Index: ginac/symmetry.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/symmetry.cpp,v retrieving revision 1.26 diff -r1.26 symmetry.cpp 29c29 < #include "numeric.h" // for factorial() --- > #include "inifcns.h" // for factorial() Index: ginsh/ginsh_parser.yy =================================================================== RCS file: /home/cvs/GiNaC/ginsh/ginsh_parser.yy,v retrieving revision 1.82 diff -r1.82 ginsh_parser.yy 7c7 < * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany --- > * GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany 319a320,329 > static ex f_abs(const exprseq &e) {return abs(e[0]);} > static ex f_acos(const exprseq &e) {return acos(e[0]);} > static ex f_acosh(const exprseq &e) {return acosh(e[0]);} > static ex f_asin(const exprseq &e) {return asin(e[0]);} > static ex f_asinh(const exprseq &e) {return asinh(e[0]);} > static ex f_atan(const exprseq &e) {return atan(e[0]);} > static ex f_atan2(const exprseq &e) {return atan2(e[0], e[1]);} > static ex f_atanh(const exprseq &e) {return atanh(e[0]);} > static ex f_beta(const exprseq &e) {return beta(e[0], e[1]);} > static ex f_binomial(const exprseq &e) {return binomial(e[0], e[1]);} 323a334,335 > static ex f_cos(const exprseq &e) {return cos(e[0]);} > static ex f_cosh(const exprseq &e) {return cosh(e[0]);} 325a338,339 > static ex f_function_derivative(const exprseq &e) {return function_derivative(e[0], e[1]);} > static ex f_eta(const exprseq &e) {return eta(e[0], e[1]);} 329a344 > static ex f_exp(const exprseq &e) {return exp(e[0]);} 330a346,348 > static ex f_factorial(const exprseq &e) {return factorial(e[0]);} > static ex f_G2(const exprseq &e) {return G(e[0], e[1]);} > static ex f_G3(const exprseq &e) {return G(e[0], e[1], e[2]);} 331a350 > static ex f_H(const exprseq &e) {return H(e[0], e[1]);} 335a355,359 > static ex f_lgamma(const exprseq &e) {return lgamma(e[0]);} > static ex f_Li2(const exprseq &e) {return Li(2, e[0]);} > static ex f_Li3(const exprseq &e) {return Li(3, e[0]);} > static ex f_Li(const exprseq &e) {return Li(e[0], e[1]);} > static ex f_log(const exprseq &e) {return log(e[0]);} 340a365 > static ex f_Order(const exprseq &e) {return Order(e[0]);} 341a367,371 > static ex f_psi1(const exprseq &e) {return psi(e[0]);} > static ex f_psi2(const exprseq &e) {return psi(e[0], e[1]);} > static ex f_S(const exprseq &e) {return S(e[0], e[1], e[2]);} > static ex f_sin(const exprseq &e) {return sin(e[0]);} > static ex f_sinh(const exprseq &e) {return sinh(e[0]);} 344a375,376 > static ex f_tan(const exprseq &e) {return tan(e[0]);} > static ex f_tanh(const exprseq &e) {return tanh(e[0]);} 345a378,380 > static ex f_tgamma(const exprseq &e) {return tgamma(e[0]);} > static ex f_zeta1(const exprseq &e) {return zeta(e[0]);} > static ex f_zeta2(const exprseq &e) {return zeta(e[0], e[1]);} 588a624,633 > {"abs", f_abs, 1}, > {"acos", f_acos, 1}, > {"acosh", f_acosh, 1}, > {"asin", f_asin, 1}, > {"asinh", f_asinh, 1}, > {"atan", f_atan, 1}, > {"atan2", f_atan2, 2}, > {"atanh", f_atanh, 1}, > {"beta", f_beta, 2}, > {"binomial", f_binomial, 2}, 595a641,642 > {"cos", f_cos, 1}, > {"cosh", f_cosh, 1}, 598a646 > {"Derivative", f_function_derivative, 2}, 603a652 > {"eta", f_eta, 2}, 609a659 > {"exp", f_exp, 1}, 610a661 > {"factorial", f_factorial, 1}, 612a664,665 > {"G", f_G2, 2}, > {"G", f_G3, 3}, 613a667 > {"H", f_H, 2}, 622a677,681 > {"lgamma", f_lgamma, 1}, > {"Li2", f_Li2, 1}, > {"Li3", f_Li3, 1}, > {"Li", f_Li, 2}, > {"log", f_log, 1}, 631a691 > {"Order", f_Order, 1}, 637a698,699 > {"psi", f_psi1, 1}, > {"psi", f_psi2, 2}, 641a704 > {"S", f_S, 3}, 642a706,707 > {"sin", f_sin, 1}, > {"sinh", f_sinh, 1}, 648a714,715 > {"tan", f_tan, 1}, > {"tanh", f_tanh, 1}, 649a717 > {"tgamma", f_tgamma, 1}, 654a723,724 > {"zeta", f_zeta1, 1}, > {"zeta", f_zeta2, 2}, 712,729d781 < static ex f_ginac_function(const exprseq &es, int serial) < { < return function(serial, es).eval(1); < } < < // All registered GiNaC functions < void GiNaC::ginsh_get_ginac_functions(void) < { < vector<function_options>::const_iterator i = function::registered_functions().begin(), end = function::registered_functions().end(); < unsigned serial = 0; < while (i != end) { < fcns.insert(make_pair(i->get_name(), fcn_desc(f_ginac_function, i->get_nparams(), serial))); < ++i; < serial++; < } < } < < 879c931 < cout << " __, _______ Copyright (C) 1999-2005 Johannes Gutenberg University Mainz,\n" --- > cout << " __, _______ Copyright (C) 1999-2006 Johannes Gutenberg University Mainz,\n" 899d950 < ginsh_get_ginac_functions();

4 9

GiNaC 1.3.4 released
by Jens Vollinga 13 Apr '06

13 Apr '06

Hi everybody, GiNaC 1.3.4 is out and available. The changes are: - More general index dimensions are now possible. - Improved algebraic substitutions. - Fixed wrong matching in .has(). - Fixed bug in differentiation of a power. - collect_common_factors also works with negative powers of common factors now. - Code clean-up and preparations for upcoming cln 1.2. As always, this release can be downloaded from http://www.ginac.de/Download.html Best wishes, Jens

1 0

Release 1.3.4
by Jens Vollinga 11 Apr '06

11 Apr '06

Hello, the time might be right to do a ginac 1.3.4 release. What do you think? Are there any reasons (aka bug->patch) to delay that release further? The NEWS file will contain the following items then: * More general index dimensions are now possible. * Improved algebraic substitutions. * Fixed bug in function bernoulli. * Fixed wrong matching in .has(). * Fixed bug in differentiation of a power. * collect_common_factors also works with negative powers of common factors now. * Code clean-up and preparations for upcoming cln 1.2. Hopefully nothing is missing. Regards, Jens

2 1

log(exp(Pi)) evaluates to Pi.
by Chris Dams 04 Apr '06

04 Apr '06

Dear all, I made log(exp(Pi)) evaluate to Pi. For this purpose I made the constructor of a constant take an optional flag to indicate the domain. I also cleaned up the code of that takes care of this simplification a bit. IMO it makes not so much sense to test here whether the x in log(exp(x)) is a symbol, constant or numeric. One can simply ask to x whether it is real. Best, Chris

1 0

Bug in differentiation of power.
by Chris Dams 04 Apr '06

04 Apr '06

Dear all, I found out that when differentiating a power it is tested whether the exponent is real. If this is the case, code is executed that assumes that this exponent is a numeric. Hower, because we also have real symbols nowadays this is not correct. Besides, why would differentiating x^I need a special case? Its derivative is simply I*x^(I-1). Therefore, I changed this test into a check whether the exponent is a numeric. I will also backport this change to 1.3. Best, Chris

1 0

Problems in compiling ginac on a iMac (Intel)
by Christian Fleck 02 Apr '06

02 Apr '06

Hello, I tried to compile ginac version 1.3.3 with gcc 4.0.0 (Apple version) on a Intel iMac running OS X 10.4.5. My version of libread is 5.1. I append the config.log. The first error is that the configure script can still not find the correct version of cln, although it should be installed correctly. I tried to play a bit, but was not able to fix it. Running 'configure --disable-clntest --enable-shared=no' and typing 'make' produces: ... ginsh_parser.yy:47:6: error: missing binary operator before token "wrapper" ginsh_parser.yy:53:6: error: missing binary operator before token "wrapper" ginsh_parser.yy:845:6: error: missing binary operator before token "wrapper" ginsh_parser.yy:855:6: error: missing binary operator before token "wrapper" I edited the file config.h changing the lines /* Major version of installed readline library. */ #define GINAC_RL_VERSION_MAJOR EditLine wrapper /* Minor version of installed readline library. */ #define GINAC_RL_VERSION_MINOR EditLine wrapper into /* Major version of installed readline library. */ #define GINAC_RL_VERSION_MAJOR 1 /* Minor version of installed readline library. */ #define GINAC_RL_VERSION_MINOR 3 But running 'make' then produced: ... Making all in ginsh if g++ -DHAVE_CONFIG_H -I. -I. -I.. -I./../ginac -I../ginac - DIN_GINAC -I/usr/local/include -g -O2 -MT ginsh_parser.o -MD -MP - MF ".deps/ginsh_parser.Tpo" -c -o ginsh_parser.o ginsh_parser.cc; \ then mv -f ".deps/ginsh_parser.Tpo" ".deps/ginsh_parser.Po"; else rm - f ".deps/ginsh_parser.Tpo"; exit 1; fi ginsh_parser.yy: In function 'char** fcn_completion(const char*, int, int)': ginsh_parser.yy:844: error: invalid conversion from 'const char*' to 'char*' ginsh_parser.yy:846: error: 'filename_completion_function' was not declared in this scope ginsh_parser.yy:846: error: 'completion_matches' was not declared in this scope ginsh_parser.yy:854: error: invalid conversion from 'const char*' to 'char*' ginsh_parser.yy:856: error: 'completion_matches' was not declared in this scope ginsh_parser.yy: In function 'int main(int, char**)': ginsh_parser.yy:908: error: invalid conversion from 'char** (*)(...)' to 'char** (*)(const char*, int, int)' ginsh_parser.yy:913: error: invalid conversion from 'const char*' to 'char*' make[2]: *** [ginsh_parser.o] Error 1 make[1]: *** [all-recursive] Error 1 make: *** [all] Error 2 Does anyone has an idea what is going wrong and how to fix it? Thanks, Christian

2 2

New .is_polynomial() method.
by Chris Dams 27 Mar '06

27 Mar '06

Dear all, I have added a method .is_polynomial() to GiNaC. Using .info(info_flags::polynomial) is not sufficient because in that case one cannot distinguish an expression that may be a polynomial in a variable but nog in another. Example: x*y*sin(y) is a polynomial in x but not in y. Hence (x*y*sin(y)).is_polynomial(x) returns true but (x*y*sin(y)).is_polynomial(lst(x,y)) returns false. I also patched the tutorial. Best, Chris

1 0