- GiNaC-devel - ginac.de

[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

Unarchiving algorithms of order N.
by Chris Dams 13 Sep '06

13 Sep '06

Dear all, I noticed that the unarchiving alogithms of container-like classes are of order N^2. This makes it impossible to unarchive large expressions within reasonable time. A patch to make this of order N is in CVS. Best wishes, Chris

1 0

Re: [GiNaC-devel] Vladimirs patch has been applied.
by Vladimir Kisil 08 Sep '06

08 Sep '06

Dear Chris, Many thanks for making the patch. From my side I revised once more clifford.cpp and change some varidx to idx. I also amend exam_clifford.cpp in away which checks metrics given by matrices both with idx and varidx. A strange thing occurred however. If I put both calls to macro CHECK6 one after another within clifford_check6 routine then compilation of exam_clifford.cpp on my computer never ends. Thus I separate them (rather ugly) to the second clone clifford_check6a. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv(a)maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/

2 1

AsyForGiNaC - an output extension producing Asymptote files
by Daniel Seidel 06 Sep '06

06 Sep '06

Hi! For the last months I worked, inspired by Vladimir Kisil, on an output extension for GiNaC to plot functions with Asymptote. Now the first version is published as a stand-alone library. My intention is to include the functions to the GiNaC library. The design is already focused on a later implementation as part of the GiNaC library. Please have a look. Tarball and Documentation can be found at <http://www.maths.leeds.ac.uk/~kisilv/asyforginac> http://www.maths.leeds.ac.uk/~kisilv/asyforginac. I would be glad to get some feedback. Greetings, Daniel.

7 28

Re: [GiNaC-devel] Vladimirs patch has been applied.
by Chris Dams 21 Aug '06

21 Aug '06

Dear all, I enabled automatic dummy index renaming also for ncmuls. This was done by moving the class make_flat_inserter from expairseq.cpp to expairseq.h and using that same class also in ncmul.cpp. After this Vladimirs test program gives as output: [[1],[2],[3],[4]].symbol8*((ONE+[[t],[x],[y],[z]].i*e.i)*e.symbol8) Best wishes, Chris

1 0

Proposal: matrices may no longer be given varidxes.
by Chris Dams 18 Aug '06

18 Aug '06

Dear all, As discussed recently on this list, if we have a matrix M, both M~0~0 and M.0.0 evaluate to the same matrix element. This idicates that the variance of a varidx given to a matrix really doesn't mean anything. Therefore I propose to add code to the method matrix::eval_indexed that checks whether the indices are idxes (in the sense of is_exactly_a<idx>) and to throw an exception if they are not. What do other developers think? Obviously this could break some code out there. Best wishes, Chris

3 3

Vladimirs patch has been applied.
by Chris Dams 18 Aug '06

18 Aug '06

Dear all, I applied Vladimirs patch from the "Improved dummy index renaming -> clifford exam fails"-thread. I hope I did everything correctly because I had to remove some casting to varidxes. This is undefined behaviour if the object being casted is not a varidx. Not to mention the operation of asking such an object to toggle its variance. I have another question for Vladimir: at some places in clifford.cpp locally within some functions a varidx is created. But now that clifford objects can also carry indices without variance, is this still correct? Maybe it should depend on the kind of index that wat supplied to such a function? Best wishes, Chris

2 2

Release of GiNaC 1.3.5
by Jens Vollinga 18 Aug '06

18 Aug '06

Hi everybody, GiNaC 1.3.5 is out and available. The changes are: - Re-built bison related files that caused bugs with gcc 4. - Fixed bugs in fderivative::match_same_type(), expairseq::match(), expand_dummy_sum(), and expairseq::map(). - Fixed template specialization for class container to be gcc 4.2 compliant. - Output in csrc mode is now -x instead of -1.0*x. - Fixed tutorial and reference generation in html/pdf/ps format. - Modified autotool configuration to allow for compilation on MinGW. - Debian packaging files updated. As always, this release can be downloaded from http://www.ginac.de/Download.html Best wishes, Jens

1 0

New: exset
by Chris Dams 17 Aug '06

17 Aug '06

Dear all, We have exvector and exmap in GiNaC, so why not have an exset too? I added it and also provided support for printing them. These containers are now be printed as: exvec: [a,b,c] exset: <a,b,c> exmap: {a==b,c==d} Best wishes, Chris

1 0

Improved dummy index renaming -> clifford exam fails
by Chris Dams 15 Aug '06

15 Aug '06

Dear Vladimir, > CD> Here one > CD> clearly sees that neither contractiong g~mu~mu nor g.mu.mu makes > CD> any sense. > > The mentioned expression indexed(squared_metric, alpha, alpha) did > not stand for a contractions, it was simply an abbreviation for > pow(e.alpha,4)=pow(g.alpha.alpha,2). However I will rewrite > clifford.cpp without this sort of "optimisation". Yes, the this interpretation was never in agreement with other code that has been in GiNaC for ages. The question has to be asked how you can guarantee that this alpha is not to be summed over? What if somewhere else in the expression there is a v~alpha? Then we would have an Einstein summation convention involving three indices, however, GiNaC always had the convention that dummy indices were to be contracted over in pairs. Also not that B~mu~mu + 1 has yielded an exception as long as I can remember when, for instance, asked what its free indices are. If you decide to go with the idea of interpreting a matrix that is given as metric to clifford objects as really being the metric in the sense as is common in physics, I have an even better optimisation for you. And it doesn't even depend on the metric being anti-commuting or not. Realizing that in that case the metric with up-indices is the inverse of the metric with down-indices, we arive at the identity e~mu e~alpha e.mu = (2 - dim) * e~alpha, where dim is the number of row/columns of the metric. > CD> This is why I was wondering in an > CD> earlier email whether we should even allow matrices to carry > CD> indices with a variance. If g is the metric tensor g with up > CD> indices should be the inverse of g with down indices. > > Why not use the agreement from the second sentence to make sense for > the first one. That certainly is a possibility. However, If I create M = [[2, 0], [0, 2]] then one of M.1.1 and M~1~1 should evaluate into 2 and the other into 1/2. Which of the two? Seems a bit arbitrary and confusing to me. Besides that, what I just wrote is not even true. It is only true if M is considered to be the metric. Not all matrices are the metric... And who says what matrix is the metric and what matrix is not... IMO it would be much less problematic to simply disallow varidxes as indices of matrices. Or maybe GiNaC should have the possibility to globally set what the metric is considered to be. Not that I would be very much in favour of that. It makes one think "global variable... blech...". > CD> If g is the metric tensor, it automatically obeys g.i~j = > CD> delta.i~j. Actually, this simplification is done automatically > CD> in GiNaC. It is in the function ex > CD> tensmetric::eval_indexed(const basic & i) const. > > So I only need take care on it if metric is given by a matrix? Well, if you are really set on giving clifford units with a matrix as metric varidxes, I suppose the answer is yes. Best wishes, Chris

2 3