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 < < #include "fderivative.h" < #include "operators.h" < #include "archive.h" < #include "utils.h" < < namespace GiNaC { < < GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(fderivative, function, < print_func(&fderivative::do_print). < print_func(&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 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 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 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(other)); < const fderivative & o = static_cast(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(other)); < const fderivative & o = static_cast(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(other)); < const fderivative & o = static_cast(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 < < #include "function.h" < < namespace GiNaC { < < < typedef std::multiset 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 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 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 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=< < #include < < // CINT needs to work properly with < #include < < #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 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_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 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 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 & 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 < inline bool is_the_function(const ex & x) < { < return is_exactly_a(x) < && ex_to(x).get_serial() == T::serial; < } < < // Check whether OBJ is the specified symbolic function. < #define is_ex_the_function(OBJ, FUNCNAME) (GiNaC::is_the_function(OBJ)) < < } // namespace GiNaC < < #endif // ndef __GINAC_FUNCTION_H__ < < END_OF_INTERFACE < < $implementation=< < #include < #include < #include < < #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 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::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 &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(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(c).delta_indent; < for (size_t i=0; i(c)) { < < // Print function name in lowercase < std::string lname = opt.name; < size_t num = lname.size(); < for (size_t i=0; i(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 1 && !(opt.symtree.is_zero())) { < exvector v = seq; < GINAC_ASSERT(is_a(opt.symtree)); < int sig = canonicalize(v.begin(), ex_to(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(serialevalf(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; iop(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 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(other)); < const function & o = static_cast(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(other)); < const function & o = static_cast(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(other)); < const function & o = static_cast(other); < < return serial == o.serial; < } < < unsigned function::return_type() const < { < GINAC_ASSERT(serialreturn_type(); < } < } < < tinfo_t function::return_type_tinfo() const < { < GINAC_ASSERT(serialreturn_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(serialsetflag(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::registered_functions() < { < static std::vector * rf = new std::vector; < 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=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::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(serialfunction.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 < #include < 26a24 > 39a38,40 > #include > #include > 43c44 < // complex conjugate --- > // absolute value 46,57c47,49 < static ex conjugate_evalf(const ex & arg) < { < if (is_exactly_a(arg)) { < return ex_to(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(&abs_function::do_print_csrc_float). > print_func(&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(conjugate_print_latex). < conjugate_func(conjugate_conjugate). < set_name("conjugate","conjugate")); < < ////////// < // absolute value < ////////// < < static ex abs_evalf(const ex & arg) < { < if (is_exactly_a(arg)) < return abs(ex_to(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(arg)) { > return abs(ex_to(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(exp) && ex_to(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(exp) && ex_to(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(abs_print_latex). < print_func(abs_print_csrc_float). < print_func(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(&conjugate_function::do_print_latex)) > > ex conjugate_function::conjugate() const 132,135c103 < if (is_exactly_a(arg)) < return step(ex_to(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(arg)) < return step(ex_to(arg)); < < else if (is_exactly_a(arg) && < is_exactly_a(arg.op(arg.nops()-1))) { < numeric oc = ex_to(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(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(arg)) { > return ex_to(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(arg)) < return csgn(ex_to(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(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(arg)) > return csgn(ex_to(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(x); < const numeric ny = ex_to(y); < const numeric nxy = ex_to(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(x); < const numeric ny = ex_to(y); < const numeric nxy = ex_to(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(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(x)) < return Li2(ex_to(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(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 const ex& arg = seq[0]; > > if (is_exactly_a(arg)) > return step(ex_to(arg)); > > else if (is_exactly_a(arg) && > is_exactly_a(arg.op(arg.nops()-1))) { > numeric oc = ex_to(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(x_pt).is_real() && ex_to(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(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 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(arg)) > return step(ex_to(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(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(&factorial_function::do_print_dflt_latex). > print_func(&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(x) || < is_exactly_a(x) || < is_exactly_a(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(x) || is_exactly_a(x) || is_exactly_a(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(factorial_print_dflt_latex). < print_func(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(y)); --- > return sym(x, ex_to(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(&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(x); --- > const mul& m = ex_to(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(r.lhs()); --- > const symbol& s = ex_to(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(&function_derivative_function::do_print_dflt). > print_func(&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(seq[0])); > GINAC_ASSERT(is_a(seq[1])); > > const ex& func = seq[1]; > ex result; > for (size_t i=0; i ex arg_diff = func.op(i).diff(s); > if (!arg_diff.is_zero()) { > lst params = ex_to(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(seq[1]).pderivative((ex_to(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(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(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].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 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 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 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 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 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 < 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 < 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(const ex& x) < { < return is_the_function(x) || is_the_function(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 < 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(const ex& x) < { < return is_the_function(x) || is_the_function(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 < 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 < 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(const ex & x) < { < return is_the_function(x) || is_the_function(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 < #include < < #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(x)) { < try { < return lgamma(ex_to(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(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(x)) { < try { < return tgamma(ex_to(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(x); < if (two_x.is_even()) { < // tgamma(n) -> (n-1)! for postitive n < if (two_x.is_positive()) { < return factorial(ex_to(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(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(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(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(x) && is_exactly_a(y)) { < try { < return tgamma(ex_to(x))*tgamma(ex_to(y))/tgamma(ex_to(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(x); < const numeric &ny = ex_to(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(rel.lhs())); < const symbol &s = ex_to(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(x)) { < try { < return psi(ex_to(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(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(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(n) && is_exactly_a(x)) { < try { < return psi(ex_to(n),ex_to(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(n); < const numeric &nx = ex_to(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 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(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 < #include < #include < #include < < #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 > 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 buf(xninitsize); < std::vector::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 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 buf(xninitsize); < std::vector::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 buf(xninitsize/2); < std::vector::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 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::const_iterator it = Xn[0].begin(); < std::vector::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::const_iterator it = Xn[n-2].begin(); < std::vector::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=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(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(cln::realpart(value))); < else if (!x.imag().is_rational()) < prec = cln::float_format(cln::the(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& s, const std::vector& x) < { < const int j = s.size(); < < std::vector 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 m_int; < std::vector x_cln; < for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) { < m_int.push_back(ex_to(*itm).to_int()); < x_cln.push_back(ex_to(*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 gsyms; < < < // type used by the transformation functions for G < typedef std::vector 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 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(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 sortmap; < int size = 0; < for (int i=0; i(abs(x[i]), i)); < ++size; < } < } < // include upper limit (scale) < sortmap.insert(std::pair(abs(y), x.nops())); < < // generate missing dummy-symbols < int i = 1; < gsyms.clear(); < gsyms.push_back(symbol("GSYMS_ERROR")); < ex lastentry; < for (std::multimap::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::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(x_) ? ex_to(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(x_) ? ex_to(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(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(x_) ? ex_to(x_) : lst(x_); < lst s = is_a(s_) ? ex_to(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(x_) ? ex_to(x_) : lst(x_); < lst s = is_a(s_) ? ex_to(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(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(m_).to_int(), ex_to(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(m_).to_int(), ex_to(x_val)); < } < } < } < // multiple polylogs < if (is_a(m_) && is_a(x_)) { < < const lst& m = ex_to(m_); < const lst& x = ex_to(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(m_)) { < if (is_a(x_)) { < // multiple polylogs < const lst& m = ex_to(m_); < const lst& x = ex_to(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(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(m_).to_int(), ex_to(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(m) || is_a(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=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(m_)) { < m = m_.op(0); < } else { < m = m_; < } < ex x; < if (is_a(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(m_)) { < m = ex_to(m_); < } else { < m = lst(m_); < } < lst x; < if (is_a(x_)) { < x = ex_to(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(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 > 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 buf(initsize); < std::vector::iterator it = buf.begin(); < std::vector::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 buf(initsize); < std::vector::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::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::iterator it = Yn[n].begin(); < std::vector::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 ynsize+1) { < for (int i=ynsize; i= 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(cln::realpart(value))); < else if (!x.imag().is_rational()) < prec = cln::float_format(cln::the(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 1) { < < cln::cl_N result; < < for (int s=0; s(x)) { < return S_num(ex_to(n).to_int(), ex_to(p).to_int(), ex_to(x)); < } else { < ex x_val = x.evalf(); < if (is_a(x_val)) { < return S_num(ex_to(n).to_int(), ex_to(p).to_int(), ex_to(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(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(n).to_int(), ex_to(p).to_int(), ex_to(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 presubsum, subsum; < subsum.push_back(0); < for (int i=1; i=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(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(e) || is_a(e)) { < return e.map(*this); < } < if (is_a(e)) { < std::string name = ex_to(e).get_name(); < if (name == "H") { < lst parameter; < if (is_a(e.op(0))) { < parameter = ex_to(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(e) || is_a(e)) { < return e.map(*this); < } < if (is_a(e)) { < std::string name = ex_to(e).get_name(); < if (name == "H") { < lst parameter; < if (is_a(e.op(0))) { < parameter = ex_to(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(e) || is_a(e)) { < return e.map(*this); < } < if (is_a(e)) { < std::string name = ex_to(e).get_name(); < if (name == "H") { < lst parameter; < if (is_a(e.op(0))) { < parameter = ex_to(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 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(h1.op(0)); < } else { < hshort = h1.op(0).op(0); < if (h2nops > 1) { < hlong = ex_to(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(e)) { < return e.map(*this); < } < < if (is_a(e)) { < < ex result = 1; < ex firstH; < lst Hlst; < for (int pos=0; pos(e.op(pos)) && is_a(e.op(pos).op(0))) { < std::string name = ex_to(e.op(pos).op(0)).get_name(); < if (name == "H") { < for (ex i=0; i(e.op(pos))) { < std::string name = ex_to(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(e)) { < name = ex_to(e).get_name(); < } < if (name == "H") { < h = e; < } else { < for (int i=0; i(e.op(i))) { < std::string name = ex_to(e.op(i)).get_name(); < if (name == "H") { < h = e.op(i); < } < } < } < } < if (h != 0) { < lst newparameter = ex_to(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(e)) { < name = ex_to(e).get_name(); < } < if (name == "H") { < h = e; < } else { < for (int i=0; i(e.op(i))) { < std::string name = ex_to(e.op(i)).get_name(); < if (name == "H") { < h = e.op(i); < } < } < } < } < if (h != 0) { < lst newparameter = ex_to(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(e)) { < name = ex_to(e).get_name(); < } < if (name == "H") { < h = e; < } else { < for (int i=0; i(e.op(i))) { < std::string name = ex_to(e.op(i)).get_name(); < if (name == "H") { < h = e.op(i); < } < } < } < } < if (h != 0) { < lst newparameter = ex_to(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(e)) { < name = ex_to(e).get_name(); < } < if (name == "H") { < h = e; < } else { < for (int i=0; i(e.op(i))) { < std::string name = ex_to(e.op(i)).get_name(); < if (name == "H") { < h = e.op(i); < } < } < } < } < if (h != 0) { < lst newparameter = ex_to(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(e)) { < name = ex_to(e).get_name(); < } < if (name == "H") { < h = e; < } else { < for (int i=0; i(e.op(i))) { < std::string name = ex_to(e.op(i)).get_name(); < if (name == "H") { < h = e.op(i); < } < } < } < } < if (h != 0) { < lst newparameter = ex_to(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(e) || is_a(e)) { < return e.map(*this); < } < < if (is_a(e)) { < std::string name = ex_to(e).get_name(); < if (name == "H") { < < lst parameter = ex_to(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; i0; 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; i0; 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(buffer)) { < for (int i=0; i 1/x transformation < struct map_trafo_H_1overx : public map_function < { < ex operator()(const ex& e) < { < if (is_a(e) || is_a(e)) { < return e.map(*this); < } < < if (is_a(e)) { < std::string name = ex_to(e).get_name(); < if (name == "H") { < < lst parameter = ex_to(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(buffer)) { < for (int i=0; i(buffer)) { < for (int i=0; i (1-x)/(1+x) transformation < struct map_trafo_H_1mxt1px : public map_function < { < ex operator()(const ex& e) < { < if (is_a(e) || is_a(e)) { < return e.map(*this); < } < < if (is_a(e)) { < std::string name = ex_to(e).get_name(); < if (name == "H") { < < lst parameter = ex_to(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(buffer)) { < for (int i=0; i(buffer)) { < for (int i=0; i& m, const cln::cl_N& x) < { < const int j = m.size(); < < std::vector 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(x1)) { < < cln::cl_N x; < if (is_a(x2)) { < x = ex_to(x2).to_cl_N(); < } else { < ex x2_val = x2.evalf(); < if (is_a(x2_val)) { < x = ex_to(x2_val).to_cl_N(); < } < } < < for (int i=0; i(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 m_int; < std::vector 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(*it_int).to_int()); < x_cln.push_back(ex_to(*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 m_int; < for (lst::const_iterator it = m_lst.begin(); it != m_lst.end(); it++) { < m_int.push_back(ex_to(*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 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(m_)) { < m = ex_to(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(m_)) { < m = ex_to(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(m_)) { < m = ex_to(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(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(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 > f_kj; < std::vector crB; < std::vector > crG; < std::vector crX; < < < void halfcyclic_convolute(const std::vector& a, const std::vector& b, std::vector& c) < { < const int size = a.size(); < for (int n=0; n& 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 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 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 >::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& 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 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& s) < { < std::vector 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 maxr) { < maxr = r[i]; < } < } < < calc_f(maxr); < < const cln::cl_N r0factorial = cln::factorial(r[0]-1); < < std::vector 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& r) < { < const int j = r.size(); < < // buffer for subsums < std::vector 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& m_, const std::vector& s_) < { < // prepare parameters < // holds Li arguments in [BBB] notation < std::vector s = s_; < std::vector m_p = m_; < std::vector m_q; < // holds Li arguments in nested sums notation < std::vector 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_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(x) && (x.nops()>1)) { < < // multiple zeta value < const int count = x.nops(); < const lst& xlst = ex_to(x); < std::vector r(count); < < // check parameters and convert them < lst::const_iterator it1 = xlst.begin(); < std::vector::iterator it2 = r.begin(); < do { < if (!(*it1).info(info_flags::posint)) { < return zeta(x).hold(); < } < *it2 = ex_to(*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(x) && (x != 1)) { < try { < return zeta(ex_to(x)); < } catch (const dunno &e) { } < } < < return zeta(x).hold(); < } < < < static ex zeta1_eval(const ex& m) < { < if (is_exactly_a(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(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(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(m_)) { < const lst& m = ex_to(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(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(x)) { < < // alternating Euler sum < const int count = x.nops(); < const lst& xlst = ex_to(x); < const lst& slst = ex_to(s); < std::vector xi(count); < std::vector si(count); < < // check parameters and convert them < lst::const_iterator it_xread = xlst.begin(); < lst::const_iterator it_sread = slst.begin(); < std::vector::iterator it_xwrite = xi.begin(); < std::vector::iterator it_swrite = si.begin(); < do { < if (!(*it_xread).info(info_flags::posint)) { < return zeta(x, s).hold(); < } < *it_xwrite = ex_to(*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(s_)) { < const lst& s = ex_to(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(m)) { < return _ex0; < } else { < if ((is_exactly_a(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(m_)) { < m = ex_to(m_); < } else { < m = lst(m_); < } < lst s; < if (is_a(s_)) { < s = ex_to(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(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 < #include < < #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(x)) < return exp(ex_to(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(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(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(x)) < return log(ex_to(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(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(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(arg.series(rel, order+extra_ord, options)); < ++extra_ord; < } while (!argser.is_terminating() && argser.nops()==1); < < const symbol &s = ex_to(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(acc)); < acc = (ex_to(rest)).mul_series(ex_to(acc)); < } < return acc; < } < const ex newarg = ex_to((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true); < return pseries(rel, seq).add_series(ex_to(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(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(x)) < return sin(ex_to(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(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(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(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(x)) < return cos(ex_to(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(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(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(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(x)) < return tan(ex_to(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(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(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(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(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(x)) < return asin(ex_to(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(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(x)) < return acos(ex_to(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(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(x)) < return atan(ex_to(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(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(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(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(y) && is_exactly_a(x)) < return atan(ex_to(y), ex_to(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(y), ex_to(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(x)) < return sinh(ex_to(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(x)); < < // sinh() is odd < if (x.info(info_flags::negative)) < return -sinh(-x); < } < < if ((x/Pi).info(info_flags::numeric) && < ex_to(x/Pi).real().is_zero()) // sinh(I*x) -> I*sin(x) < return I*sin(x/I); < < if (is_exactly_a(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(x)) < return cosh(ex_to(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(x)); < < // cosh() is even < if (x.info(info_flags::negative)) < return cosh(-x); < } < < if ((x/Pi).info(info_flags::numeric) && < ex_to(x/Pi).real().is_zero()) // cosh(I*x) -> cos(x) < return cos(x/I); < < if (is_exactly_a(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(x)) < return tanh(ex_to(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(x)); < < // tanh() is odd < if (x.info(info_flags::negative)) < return -tanh(-x); < } < < if ((x/Pi).info(info_flags::numeric) && < ex_to(x/Pi).real().is_zero()) // tanh(I*x) -> I*tan(x); < return I*tan(x/I); < < if (is_exactly_a(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(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(x)) < return asinh(ex_to(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(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(x)) < return acosh(ex_to(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(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(x)) < return atanh(ex_to(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(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(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(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($3)).eval(1); --- > factory_p p = find_func_factory(n+"_function"); > // exprseq -> exvector > const exprseq& in = ex_to($3); > exvector out; > for (exprseq::const_iterator it = in.begin(); it != in.end(); ++it) { > out.push_back(*it); > } > $$ = (*p)(out); > // $$ = n(ex_to($3)); > //TODO > // unsigned i = function::find_function(n, $3.nops()); > // $$ = function(i, ex_to($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(ebasis)) < return ex_to(ebasis).power(eexponent); --- > if (is_a(ebasis)) > return ex_to(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(i->rest)) { 270c270 < if (is_order_function(seq[i].rest)) --- > if (is_exactly_a(seq[i].rest)) 502c502 < if (is_order_function(it->rest)) { --- > if (is_exactly_a(it->rest)) { 515c515 < if (is_order_function(it->rest)) { --- > if (is_exactly_a(it->rest)) { 535c535 < if (is_order_function(it->rest)) { --- > if (is_exactly_a(it->rest)) { 547c547 < return seq.empty() || !is_order_function((seq.end()-1)->rest); --- > return seq.empty() || !is_exactly_a((seq.end()-1)->rest); 682c682 < if (is_order_function((*a).rest)) --- > if (is_exactly_a((*a).rest)) 688c688 < if (is_order_function((*b).rest)) --- > if (is_exactly_a((*b).rest)) 693c693 < if (is_order_function((*a).rest) || is_order_function((*b).rest)) { --- > if (is_exactly_a((*a).rest) || is_exactly_a((*b).rest)) { 750c750 < if (!is_order_function(it->rest)) --- > if (!is_exactly_a(it->rest)) 791c791 < if (is_order_function(coeff(var, a_max))) --- > if (is_exactly_a(coeff(var, a_max))) 793c793 < if (is_order_function(other.coeff(var, b_max))) --- > if (is_exactly_a(other.coeff(var, b_max))) 805c805 < if (!is_order_function(a_coeff) && !is_order_function(b_coeff)) --- > if (!is_exactly_a(a_coeff) && !is_exactly_a(b_coeff)) 940c940 < if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative()) --- > if (seq.size() == 1 && is_exactly_a(seq[0].rest) && p.real().is_negative()) 951c951 < if (is_order_function(c)) { --- > if (is_exactly_a(c)) { 966c966 < if (is_order_function(co[i])) { --- > if (is_exactly_a(co[i])) { 1117c1117 < if (is_order_function(currcoeff)) --- > if (is_exactly_a(currcoeff)) 1133c1133 < if (is_order_function(currcoeff)) --- > if (is_exactly_a(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& 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 < < #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=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()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 & remember_table::remember_tables() < { < static std::vector * rt = new std::vector; < 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 < #include < #include < < 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 { < 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 { < 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_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::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();