October 2006 - GiNaC-devel

Custom power evaluation
by Vladimir Kisil 31 Oct '06

31 Oct '06

Dear All, In my GiNaC code I oftenly use some objects which have special properties with respect to power. For example, if values of a symbol s are assumed to be only -1, 0, or 1, then s^3=s. Another example is provided by "imaginary units" for double and dual numbers. It will be good for such objects to use such simplifying identities on the fly in power_eval() method rather then by user calls to subs() afterwards. I propose a patch to GiNaC which allows to supply such simplification rules from a user-defined classes. It requires only two new virtual methods to basic class and should not involve much overheads to the general GiNaC performance. I give an example below (which should work with the patch) for a symbol which is assumed to take only values -1 or 1, thus satisfying to s^2=1 rule. If general feeling toward this proposition will be positive I would send a more complete patch including documentation. Best wishes, Vladimir PS Somehow I was not able to use macros DEFAULT_UNARCHIVE() and DEFAULT_COMPARE() for my class, since compiler produces the error: error: expected constructor, destructor, or type conversion before 'int' What may be the reason for this? -- Vladimir V. Kisil email: kisilv(a)maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/ #include <iostream> #include <stdexcept> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; class myrealsymb : public symbol { GINAC_DECLARE_REGISTERED_CLASS(myrealsymb, symbol) public: myrealsymb(const std::string & initname); virtual ex eval_power_basis(const ex & exp) const; virtual bool info(unsigned inf) const; }; GINAC_IMPLEMENT_REGISTERED_CLASS(myrealsymb,symbol) myrealsymb::myrealsymb(): symbol(){}; myrealsymb::myrealsymb(const std::string & initname) : symbol(initname) {;} ex myrealsymb::eval_power_basis(const ex & exp) const { if (is_a<numeric>(exp)) { if (ex_to<numeric>(exp).is_even()) return numeric(1); if (ex_to<numeric>(exp).is_odd()) return *this; } return power(*this,exp).hold(); } bool myrealsymb::info(unsigned inf) const { return (inf == info_flags::is_power_basis) || symbol::info(inf); } myrealsymb::myrealsymb(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) { } void myrealsymb::archive(archive_node &n) const { symbol::archive(n); } ex myrealsymb::unarchive(const archive_node &n, lst &sym_lst) { return (new myrealsymb(n, sym_lst))->symbol::setflag(status_flags::dynallocated); } int myrealsymb::compare_same_type(const basic & other) const { return symbol::compare_same_type(other); } int main(){ try { ex e=myrealsymb("s"); cout << pow(e,2) << endl; // -> 1 cout << pow(1+e,2).expand() << endl; // -> 2+2*s } catch (exception &p) { cerr << "Got problem: " << p.what() << endl; } return 0; }

4 7

Custom power evaluation
by Vladimir Kisil 27 Oct '06

27 Oct '06

Dear All, Following my previous message on the subject (http://www.ginac.de/pipermail/ginac-devel/2006-October/001071.html) I am including the full patch. Its functionality is now described in ginac.info. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv(a)maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/

2 3

Re: Bug(?) in reposition_dummy_indices: test case
by Chris Dams 24 Oct '06

24 Oct '06

Dear Alexei, > First of all, I think it is perfectly sensible for eval() to return zero > (e.g., what about traceless tensors?). Yes. However, if we have a traceless tensor both T~mu.mu and T.mu~mu should evaluated themselves into 0. It does not make sense to only have one of them evaluate into zero. The point I was trying to make is that tensors with some index lowered/raised always (?) have the same properties so that an evaluate step that aplies to the version with an up-index should also aply to the one with a down-index. > Secondly, your patch does not catch this: > [I admit that example is a little bit weird] Your example is perfectly fine. I should have deleted all indices that are dummy indices of the indexed object e that is given as a parameter to reposition_dummy_indices in any case, not just when a cyclic symmetry is involved. Because you also mentioned the possibility that users would want to implement things like the Riemann curvature tensor, it makes sense that the evaluation of an indexed object could do any weird kind of symmetrization it would want, so I also removed the condition that a cyclic symmetry should be involved when trying all combinations of repositioning local dummy indices. > > On the other hand, as a matter of principle you are right that it is > > best to assume as little about the rest of the code as possbile. For > > this reason I would accept a patch that would do the substitutions as > > they are currently but does not cause eval to be called. > > The patch I've posted does that. Is there something wrong/stupid/etc > with it? Oh well, the thing seemed rather complicated to me and didn't do all simplifications possible, as noted in http://www.ginac.de/pipermail/ginac-devel/2006-August/001055.html. I now more or less applied it, but with some modifications. I hope that I haven't broken anything ;-). (1) Now that the raising and lowering of indices is done on an exvector, we do not need to raise/lower both indices at the same time. Another reason why this isn't necessary anymore is that indices that occur in dummy pairs within the indexed object e should already have been removed from the vector variant_dummy_indices. I think this means that the for_each that you had can be omitted. (2) It is not correct to return the indexed object by doing indexed ei_ = ex_to<indexed>(e); ei_.seq = seq; e = ei_; This should be e = ex_to<indexed>(e).thiscontainer(seq); because we may be dealing with something that inherits from indexed rather than an indexed object. I also simplified the clifford exam by removing a huge macro. A patch is in CVS. Best wishes, Chris

2 1

[BUG] clifford::op() versus tensor::replace_contr_index
by varg＠theor.jinr.ru 18 Oct '06

18 Oct '06

Dear all, I've tried making ex_to<T> functions [more] type-safe, like this: diff --git a/ginac/ex.h b/ginac/ex.h index 3070b7a..79f872b 100644 --- a/ginac/ex.h +++ b/ginac/ex.h @@ -949,7 +949,8 @@ inline bool is_exactly_a(const ex &obj) template <class T> inline const T &ex_to(const ex &e) { - GINAC_ASSERT(is_a<T>(e)); + if (__builtin_expect(! is_a<T>(e), 0)) + throw std::bad_cast(); return static_cast<const T &>(*e.bp); } and got failure in exam_clifford (to be more specific, last expression in clifford_check4 triggered bad_cast). I've found a simple[r] example which fails in the same way: #include <iostream> #include <stdexcept> #include "ginac.h" using namespace GiNaC; int main(int argc, char** argv) { symbol d("D"); varidx mu(symbol("mu"), d), nu(symbol("nu"), d); ex e = lorentz_g(mu, nu)*dirac_gamma(nu.toggle_variance()); e = e.simplify_indexed(); return 0; } The reason of this error seems to be inconsistency between clifford::op and tensor::replace_contr_index. tensor::replace_contr_index expects any subexpression of `indexed' referenced via op() to be `idx' (see tensor.cpp:398). On the other hand, `clifford' objects make their representation label accessible via op(), so this assumption is incorrect. So, the questions are: 1) Any ideas how to fix this? 2) In general, is it OK for class derived from indexed return something which is *not* idx (or derived from it) via op()? Best regards, Alexei -- All science is either physics or stamp collecting.

2 3

Bug(?) in reposition_dummy_indices: test case and patch
by varg＠theor.jinr.ru 10 Oct '06

10 Oct '06

Hello! This simple program #include <iostream> #include <stdexcept> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main(int argc, char** argv) { varidx mu1(symbol("mu1"), 4); varidx mu2(symbol("mu2"), 4); varidx mu3(symbol("mu3"), 4); varidx mu4(symbol("mu4"), 4); varidx mu5(symbol("mu5"), 4); varidx mu6(symbol("mu6"), 4); exvector ev2; ev2.push_back(mu3.toggle_variance()); ev2.push_back(mu6); ev2.push_back(mu5.toggle_variance()); ev2.push_back(mu6.toggle_variance()); ev2.push_back(mu5); ev2.push_back(mu3); ex test_cycl = indexed(symbol("A"), sy_cycl(), ev2); test_cycl = test_cycl.simplify_indexed(); if (test_cycl.get_free_indices().size() != 0) { std::cerr << e << std::endl; throw std::logic_error("Inconsistent indices in expression"); } return 0; } fails (both with GiNaC 1.4 CVS and 1.3.x). I believe that the reason is bug in indexed.cpp:reposition_dummy_indices() and propose attached patch to fix it. Note that patch *seems* to be correct, but IMHO it is somewhat ugly and probably inefficient. Could anyone suggest a better solution? Best regards, Alexei. -- All science is either physics or stamp collecting.

2 7

[PATCH 3/3] power::to_polinomial: improve handling of negative powers
by Alexei Sheplyakov 06 Oct '06

06 Oct '06

Now collect_common_factors((x*y*a+x*y*b)^(-3) + (x*z + x*y)^(-2)) will evaluate to x^(-2)*(x^(-1)*(a+b)^(-3)*y^(-3)+(z+y)^(-2)) which is actually what I expect. --- ginac/normal.cpp | 11 ++++++++++- 1 files changed, 10 insertions(+), 1 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index fc89710..b2a4e8c 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -2455,7 +2455,16 @@ ex power::to_polynomial(exmap & repl) co if (exponent.info(info_flags::posint)) return power(basis.to_rational(repl), exponent); else if (exponent.info(info_flags::negint)) - return power(replace_with_symbol(power(basis, _ex_1), repl), -exponent); + { + ex basis_pref = collect_common_factors(basis); + if (is_exactly_a<mul>(basis_pref) || is_exactly_a<power>(basis_pref)) { + // (A*B)^n will be automagically transformed to A^n*B^n + ex t = power(basis_pref, exponent); + return t.to_polynomial(repl); + } + else + return power(replace_with_symbol(power(basis, _ex_1), repl), -exponent); + } else return replace_with_symbol(*this, repl); } -- All science is either physics or stamp collecting.

1 0

[PATCH 2/3] collect_common_factors: improve handling of powers.
by Alexei Sheplyakov 06 Oct '06

06 Oct '06

Now collect_common_factors((a*x+a*y)^(-2)) evaluates to (x+y)^(-2)*a^(-2) as one might expect. --- ginac/normal.cpp | 4 ++-- 1 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 49f0bd7..fc89710 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -2513,7 +2513,7 @@ static ex find_common_factor(const ex & for (size_t i=0; i<num; i++) { ex x = e.op(i).to_polynomial(repl); - if (is_exactly_a<add>(x) || is_exactly_a<mul>(x)) { + if (is_exactly_a<add>(x) || is_exactly_a<mul>(x) || is_a<power>(x)) { ex f = 1; x = find_common_factor(x, f, repl); x *= f; @@ -2573,7 +2573,7 @@ term_done: ; } else if (is_exactly_a<power>(e)) { const ex e_exp(e.op(1)); - if (e_exp.info(info_flags::posint)) { + if (e_exp.info(info_flags::integer)) { ex eb = e.op(0).to_polynomial(repl); ex factor_local(_ex1); ex pre_res = find_common_factor(eb, factor_local, repl); -- All science is either physics or stamp collecting.

1 0

[PATCH 1/3] divide: try to avoid expanding partially factored expressions.
by Alexei Sheplyakov 06 Oct '06

06 Oct '06

Dear all, this ginsh code snippet: collect_common_factors((x*y*a+x*y*b)^(3) + (x*z + x*y)^(2)) gives x^2*(x*(3*a^2*y^3*b+3*a*y^3*b^2+y^3*b^3+a^3*y^3)+2*z*y+z^2+y^2) Thus instead of collecting common factors from the terms of sums (as manual says) collect_common_factors [sometimes] expands term[s]! --- ginac/normal.cpp | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 files changed, 66 insertions(+), 0 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 7b90030..49f0bd7 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -614,6 +614,72 @@ #endif if (!get_first_symbol(a, x) && !get_first_symbol(b, x)) throw(std::invalid_argument("invalid expression in divide()")); + // Try to avoid expanding partially factored expressions. + if (is_exactly_a<mul>(b)) { + // Divide sequentially by each term + ex rem_new, rem_old = a; + for (size_t i=0; i < b.nops(); i++) { + if (! divide(rem_old, b.op(i), rem_new, false)) + return false; + rem_old = rem_new; + } + q = rem_new; + return true; + } else if (is_exactly_a<power>(b)) { + const ex& bb(b.op(0)); + int exp_b = ex_to<numeric>(b.op(1)).to_int(); + ex rem_new, rem_old = a; + for (int i=exp_b; i>0; i--) { + if (! divide(rem_old, bb, rem_new, false)) + return false; + rem_old = rem_new; + } + q = rem_new; + return true; + } + + if (is_exactly_a<mul>(a)) { + // Divide sequentially each term. If some term in a is divisible + // by b we are done... and if not, we can't really say anything. + size_t i; + ex rem_i; + bool divisible_p = false; + for (i=0; i < a.nops(); ++i) { + if (divide(a.op(i), b, rem_i, false)) { + divisible_p = true; + break; + } + } + if (divisible_p) { + exvector resv; + resv.reserve(a.nops()); + for (size_t j=0; j < a.nops(); j++) { + if (j==i) + resv.push_back(rem_i); + else + resv.push_back(a.op(j)); + } + q = (new mul(resv))->setflag(status_flags::dynallocated); + return true; + } + } else if (is_exactly_a<power>(a)) { + // The base itself might be divisible by b, in that case we don't + // need to expand a + const ex& ab(a.op(0)); + int a_exp = ex_to<numeric>(a.op(1)).to_int(); + ex rem_i; + if (divide(ab, b, rem_i, false)) { + q = rem_i*power(ab, a_exp - 1); + return true; + } + for (int i=2; i < a_exp; i++) { + if (divide(power(ab, i), b, rem_i, false)) { + q = rem_i*power(ab, a_exp - i); + return true; + } + } // ... so we *really* need to expand expression. + } + // Polynomial long division (recursive) ex r = a.expand(); if (r.is_zero()) { -- All science is either physics or stamp collecting.

1 0