- GiNaC-devel - ginac.de

new function class
by Jens Vollinga 15 Feb '06

15 Feb '06

Hi, here is a short example code to demonstrate the changes I plan to make to GiNaC's function system (like discussed in November on this list). - every function is now a C++ class derived from class function (here still called newfunction to avoid name clashes). - the 14 functions already present in the standard C library have a special naming convention (here _ginac is appended to the class name) and supporting C functions for creation. - the definition of the functions is rather long, but will be made more comfortable with the help of new macros like GINAC_DECLARE_FUNCTION ... (not included here yet). Comments? Regards, Jens #include <iostream> #include <fstream> #include <stdexcept> #include <cmath> #include <ginac/ginac.h> namespace GiNaC { //////////////////////////////////////////////////////////////////////////////// class newfunction : public exprseq { GINAC_DECLARE_REGISTERED_CLASS_NO_CTORS(newfunction, exprseq) virtual newfunction* duplicate() const { return new newfunction(*this); } virtual void accept(GiNaC::visitor& v) const { if (visitor* p = dynamic_cast<visitor*>(&v)) { p->visit(*this); } else { inherited::accept(v); } } protected: virtual int compare_same_type(const basic& other) const { return exprseq::compare_same_type(other); } virtual bool is_equal_same_type(const basic& other) const { return exprseq::compare_same_type(other); } public: virtual void print(const print_context& c, unsigned level = 0) const; public: newfunction(tinfo_t ti) { tinfo_key = ti; } newfunction(tinfo_t ti, const ex& x1) : inherited(x1) { tinfo_key = ti; } newfunction(tinfo_t ti, const ex& x1, const ex& x2) : inherited(x1, x2) { tinfo_key = ti; } newfunction(tinfo_t ti, const ex& x1, const ex& x2, const ex& x3) : inherited(x1, x2, x3) { tinfo_key = ti; } // ... }; GINAC_IMPLEMENT_REGISTERED_CLASS(newfunction, exprseq) newfunction::newfunction(const GiNaC::archive_node& n, GiNaC::lst& sym_lst) : inherited(n, sym_lst) { } void newfunction::archive(GiNaC::archive_node& n) const { inherited::archive(n); } GiNaC::ex newfunction::unarchive(const GiNaC::archive_node& n, GiNaC::lst& sym_lst) { return (new newfunction(n, sym_lst))->setflag(status_flags::dynallocated); } void newfunction::print(const print_context& c, unsigned level) const { const std::vector<print_functor>& pdt = get_class_info().options.get_print_dispatch_table(); unsigned id = c.get_class_info().options.get_id(); if (id >= pdt.size() || !(pdt[id].is_valid())) { if (is_a<print_tree>(c)) { c.s << std::string(level, ' ') << class_name() << " @" << this << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec << ", nops=" << nops() << std::endl; unsigned delta_indent = static_cast<const print_tree&>(c).delta_indent; for (size_t i=0; i<seq.size(); ++i) { seq[i].print(c, level + delta_indent); } c.s << std::string(level + delta_indent, ' ') << "=====" << std::endl; } else if (is_a<print_latex>(c)) { c.s << "\\mbox{" << class_name() << "}"; inherited::do_print(c,level); } else { c.s << class_name(); inherited::do_print(c,level); } } else { pdt[id](*this, c, level); } } //////////////////////////////////////////////////////////////////////////////// class polylog : public newfunction { // a new GINAC_DECLARE_FUNCTION macro will replace this GINAC_DECLARE_REGISTERED_CLASS_NO_CTORS(polylog, newfunction) // will be part of GINAC_DECLARE_FUNCTION macro virtual polylog* duplicate() const { return new polylog(*this); } // will be part of GINAC_DECLARE_FUNCTION macro virtual void accept(GiNaC::visitor& v) const { if (visitor* p = dynamic_cast<visitor*>(&v)) { p->visit(*this); } else { inherited::accept(v); } } public: polylog(const ex& x1, const ex& x2) : inherited(&polylog::tinfo_static, x1, x2) { } protected: void do_print_latex(const print_context& c, unsigned level) const; public: virtual ex eval(int level = 0) const; }; GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(polylog, exprseq, print_func<print_latex>(&polylog::do_print_latex)) // will be part of GINAC_IMPLEMENT_FUNCTION macro polylog::polylog(const GiNaC::archive_node& n, GiNaC::lst& sym_lst) : inherited(n, sym_lst) { } // will be part of GINAC_IMPLEMENT_FUNCTION macro void polylog::archive(GiNaC::archive_node& n) const { inherited::archive(n); } // will be part of GINAC_IMPLEMENT_FUNCTION macro GiNaC::ex polylog::unarchive(const GiNaC::archive_node& n, GiNaC::lst& sym_lst) { return (new polylog(n, sym_lst))->setflag(status_flags::dynallocated); } void polylog::do_print_latex(const print_context& c, unsigned level) const { c.s << "\\mbox{Li}"; inherited::do_print(c,level); } ex polylog::eval(int level) const { if (op(0).is_equal(1)) { return zeta(op(1)); } else { return this->hold(); } } //////////////////////////////////////////////////////////////////////////////// class log10_ginac : public newfunction { GINAC_DECLARE_REGISTERED_CLASS_NO_CTORS(log10_ginac, newfunction) virtual log10_ginac* duplicate() const { return new log10_ginac(*this); } virtual void accept(GiNaC::visitor& v) const { if (visitor* p = dynamic_cast<visitor*>(&v)) { p->visit(*this); } else { inherited::accept(v); } } public: log10_ginac(const ex& x) : inherited(&log10_ginac::tinfo_static, x) { } protected: void do_print(const print_context& c, unsigned level) const; void do_print_latex(const print_context& c, unsigned level) const; public: virtual ex eval(int level = 0) const; }; GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(log10_ginac, newfunction, print_func<print_dflt>(&log10_ginac::do_print). print_func<print_latex>(&log10_ginac::do_print_latex)) log10_ginac::log10_ginac(const GiNaC::archive_node& n, GiNaC::lst& sym_lst) : inherited(n, sym_lst) { } void log10_ginac::archive(GiNaC::archive_node& n) const { inherited::archive(n); } GiNaC::ex log10_ginac::unarchive(const GiNaC::archive_node& n, GiNaC::lst& sym_lst) { return (new log10_ginac(n, sym_lst))->setflag(status_flags::dynallocated); } void log10_ginac::do_print(const print_context& c, unsigned level) const { std::string classname(class_name()); c.s << classname.erase(classname.find("_ginac",0),6); inherited::do_print(c,level); } void log10_ginac::do_print_latex(const print_context& c, unsigned level) const { c.s << "\\log"; inherited::do_print(c,level); } ex log10_ginac::eval(int level) const { if (op(0).is_equal(0)) { return 666; } else { return this->hold(); } } template<typename T1> const GiNaC::log10_ginac log10(const T1& x) { return log10_ginac(x); } const GiNaC::log10_ginac log10(double x) { return log10_ginac(x); } const GiNaC::log10_ginac log10(float x) { return log10_ginac(x); } //////////////////////////////////////////////////////////////////////////////// } // namespace GiNaC //using namespace std; using namespace GiNaC; int main() { try { symbol x("x"); ex r1 = polylog(x,3) + polylog(1,3) + 6*GiNaC::log10(3.0f) - log10(x); std::cout << dflt << r1 << std::endl; std::cout << tree << r1 << std::endl; std::cout << latex << r1 << std::endl; } catch (const std::exception& e) { std::cout << e.what() << std::endl; } return 0; }

2 6

Fixed problem on 64-bit machines and introduced has_options::algebraic.
by Chris Dams 14 Feb '06

14 Feb '06

Dear all, I found out that the new tinfo() system does not compile on AMD 64-bit machines because integers are still 32 bit on them. This makes the cast from void* to unsigned illegal. I solved this problem. I also introduced an extra optional argument for .has() that can be used to pass a has_option. I defined has_options::algebraic. This is analogous to subs_options::algebraic. For example (x*y*z).has(x*y,has_options::algebraic) will return true. I also included some checks for correct behavior of this feature in the check directory and documented it in the tutorial. Best, Chris

1 0

Index dimensions.
by Chris Dams 11 Jan '06

11 Jan '06

Dear all, I have modified the function minimal_dim in idx.cpp so that non-numerical dimensions do not have to be symbols. For instance minimal_dim(4,4-2*delta) will return 4. Generally a numeric dimension will always be considered smaller than a non-numeric dimension. This makes sense in the context of dimensional regularization. Best wishes, Chris

2 2

Another bug in clifford.cpp
by Chris Dams 05 Jan '06

05 Jan '06

Dear all, I found out that canonicalizing clifford objects was not functioning properly. This was because a dirac_slash was automatically getting zero as its metric. When canonicalizing clifford objects the anti-commutator of two objects that were ordered in the wrong way is necesary. However, if the metric is automatically set to zero, this anti-commutator is zero too. I fixed this in CVS. I do have two questions about the implementation (1) Why are the indices in dirac_gamma (and with my patch in dirac_slash) static? Does the dimension of them not matter? It now seems that the dimension as used in the metric is always equal to the dimension of the very first dirac_gamma/dirac_slash that is constructed. (2) Why does the metric have indices? In view of the implementation of get_metric the primary purpose of these indices seems to be to be replaced by other indices. If that is the case why are they present in the first place? Best, Chris

2 3

Bug in clifford::subs().
by Chris Dams 04 Jan '06

04 Jan '06

Dear all, I found out that if there is nothing to substitute for clifford::subs it does not return an expression that is trivially equal (in the sense of are_ex_trivially_equal) to the original one. Besides that there were a lot of overloaded clifford::subs methods of which the purpose is unclear to me. I removed these spurious methods and fixed the bug. The patch is in CVS. Best wishes, Chris

1 0

Improvements for algebraic substitutions.
by Chris Dams 03 Jan '06

03 Jan '06

Happy new year! I found out that if we do (notation simplified to be more readable but less compilable): (a^n*b^m*b^p).subs($0^$1*$0^$2==$0^($1+$2),subs_options::algebraic) nothing happens because a^n already matches a factor in the pattern and in that case no other matches are attempted. I fixed this. Also, I allowed algebraic substitution to match multiple times on different factors with the same pattern. E.g. the above .subs() can be used to simplify a^n*a^m*b^p*b^q -> a^(n+m)*b^(p+q) These changes are in CVS. Bye, Chris

1 0

Bugfixes in 1.3
by Chris Dams 16 Dec '05

16 Dec '05

Dear GiNaCers, I put the speed-up for simplify_indexed and the fix for the yestday-discovered bug for user-defined inner products in CVS. Best, Chris

1 0

Improvements for dummy index renaming.
by Chris Dams 16 Dec '05

16 Dec '05

Dear all, I found out that the function get_all_dummy_indices only really gets all dummy indices if the it acts on something that is expanded (i.e., it works fine when called form simplify_indexed). However, if it is called from the new index renaming code, it does not see dummy indices that occur in a factor that looks, for instance like (A.i*A.i+alpha). To make this work, I added a new function get_all_dummy_indices_safely which is more costly computatunally but is supposed to really get always all dummy indices. Other changes: (1) Removed power::get_free_indices(). It is a bad idea that and index inside a power could be contracted with something. Therefore it should not be considered a free index. Besides, now that all dummy indices are found, this would break algebraic substitutions when used for something like U.i.j*Uinv.j.k==delta.i.k. (2) As mentioned in the "Note:" of my previous mail, removed the try/catch block that was introduce there. Now we do not try to look up inner products of expressions with more than one index. (3) Please don't cast an expession to a lst with (lst)expression. This is undefined behaviour and at my place it does not work as intended. What was the purpose of ex_to<whatever> again? Best, Chris

1 0

Crash in simplify_indexed.
by Chris Dams 15 Dec '05

15 Dec '05

Dear all, I found out that if one tries to use simplify_indexed for a product that contains, say, T~a.mu*T~a~mu and the dimensions of the indices a and mu are incomparable, this throws an exception. This is because GiNaC tries to apply user defined scalar_products in order to contract factors. As the user may not even be using scalar_products, I think this is not a good idea. I solved this by wrapping the part of the code that does this in a try-catch block. A patch is in CVS. Note: The try/catch-approach is not very elegant IMO. An alternative solution might be to stop trying to use user-defined scalar_products in case of considering contraction of indexed objects with more than one index. Is anybody using such a feature anyway? Best, Chris

2 1

tinfo
by Jens Vollinga 30 Nov '05

30 Nov '05

Hi everybody, just an idea for a better(?) tinfo mechanism in ginac: if we added two static constants to every class like class SOMECLASS { public: const static char* tinfo_name; const static unsigned tinfo; // ... rest of class ... } const char* SOMECLASS::tinfo_name = "SOMECLASS"; const unsigned SOMECLASS::tinfo = (unsigned)SOMECLASS::idstring; it seems that we would get a tinfo mechanism for ginac without having to set the tinfo numbers by hand. But maybe I am mistaken. So I'd like to ask somebody who knows C++ better than me, if there might be some hidden pitfalls. Is the conversion a problem (what about 64bit machines)? Is it guaranteed that the address of the static string is unique to every class? Maybe some optimizer will break this? Does it only work with gcc on linux etc.? Regards, Jens

3 12