- GiNaC-devel - ginac.de

[PATCH] {mul,power}::eval(): remove extra recursions
by Alexei Sheplyakov 19 Aug '07

19 Aug '07

Hello! The patch below makes code converting terms into primitive polynomials (i.e. (-x+y)*(x-y) -> -(x-y)^2, (3*x+y)/(6*x+2*y)^2 -> 1/4*(y+3*x)^(-1)) a little bit more efficient: numeric coefficient is calculated in one sweep, sums divided by numbers in place. While at it, fix minor signed vs unsiged issue. --- ginac/compiler.h | 12 ++++++++++++ ginac/mul.cpp | 48 ++++++++++++++++++++++++++++++++++-------------- ginac/power.cpp | 33 +++++++++++++++++++-------------- 3 files changed, 65 insertions(+), 28 deletions(-) diff --git a/ginac/compiler.h b/ginac/compiler.h new file mode 100644 index 0000000..66f6e92 --- /dev/null +++ b/ginac/compiler.h @@ -0,0 +1,12 @@ +#ifndef GINAC_COMPILER_DEP_HH +#define GINAC_COMPILER_DEP_HH + +#ifdef __GNUC__ +#define unlikely(cond) __builtin_expect((cond), 0) +#define likely(cond) __builtin_expect((cond), 1) +#else +#define unlikely(cond) (cond) +#define likely(cond) (cond) +#endif + +#endif /* GINAC_COMPILER_DEP_HH */ diff --git a/ginac/mul.cpp b/ginac/mul.cpp index b4ef897..9ed27c6 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -34,6 +34,7 @@ #include "lst.h" #include "archive.h" #include "utils.h" +#include "compiler.h" namespace GiNaC { @@ -422,7 +423,7 @@ ex mul::eval(int level) const return *this; } - int seq_size = seq.size(); + size_t seq_size = seq.size(); if (overall_coeff.is_zero()) { // *(...,x;0) -> 0 return _ex0; @@ -448,14 +449,18 @@ ex mul::eval(int level) const ex_to<numeric>(addref.overall_coeff). mul_dyn(ex_to<numeric>(overall_coeff))) )->setflag(status_flags::dynallocated | status_flags::evaluated); - } else if (seq_size >= 2) { + } else if ((seq_size >= 2) && (! (flags & status_flags::expanded))) { // Strip the content and the unit part from each term. Thus // things like (-x+a)*(3*x-3*a) automagically turn into - 3*(x-a)2 epvector::const_iterator last = seq.end(); epvector::const_iterator i = seq.begin(); + epvector::const_iterator j = seq.begin(); + std::auto_ptr<epvector> s(new epvector); + numeric oc = *_num1_p; + bool something_changed = false; while (i!=last) { - if (! (is_a<add>(i->rest) && i->coeff.is_equal(_ex1))) { + if (likely(! (is_a<add>(i->rest) && i->coeff.is_equal(_ex1)))) { // power::eval has such a rule, no need to handle powers here ++i; continue; @@ -472,32 +477,47 @@ ex mul::eval(int level) const // very unlucky event it can even loop forever). Hopefully the main // variable will be the same for all terms in *this const bool unit_normal = lead_coeff.is_pos_integer(); - if ((c == *_num1_p) && ((! canonicalizable) || unit_normal)) { + if (likely((c == *_num1_p) && ((! canonicalizable) || unit_normal))) { ++i; continue; } - std::auto_ptr<epvector> s(new epvector); - s->reserve(seq.size()); + if (! something_changed) { + s->reserve(seq_size); + something_changed = true; + } - epvector::const_iterator j=seq.begin(); - while (j!=i) { + while ((j!=i) && (j!=last)) { s->push_back(*j); ++j; } - if (! unit_normal) { + if (! unit_normal) c = c.mul(*_num_1_p); - } - const ex primitive = (i->rest)/c; - s->push_back(expair(primitive, _ex1)); - ++j; + oc = oc.mul(c); + + // divide add by the number in place to save at least 2 .eval() calls + const add& addref = ex_to<add>(i->rest); + add* primitive = new add(addref); + primitive->setflag(status_flags::dynallocated); + primitive->clearflag(status_flags::hash_calculated); + primitive->overall_coeff = ex_to<numeric>(primitive->overall_coeff).div_dyn(c); + for (epvector::iterator ai = primitive->seq.begin(); + ai != primitive->seq.end(); ++ai) + ai->coeff = ex_to<numeric>(ai->coeff).div_dyn(c); + + s->push_back(expair(*primitive, _ex1)); + + ++i; + ++j; + } + if (something_changed) { while (j!=last) { s->push_back(*j); ++j; } - return (new mul(s, ex_to<numeric>(overall_coeff).mul_dyn(c)) + return (new mul(s, ex_to<numeric>(overall_coeff).mul_dyn(oc)) )->setflag(status_flags::dynallocated); } } diff --git a/ginac/power.cpp b/ginac/power.cpp index 8829bba..354519b 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -41,6 +41,7 @@ #include "archive.h" #include "utils.h" #include "relational.h" +#include "compiler.h" namespace GiNaC { @@ -480,25 +481,29 @@ ex power::eval(int level) const // (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4) if (num_exponent->is_integer() && is_exactly_a<add>(ebasis)) { - const numeric icont = ebasis.integer_content(); + numeric icont = ebasis.integer_content(); const numeric& lead_coeff = ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div_dyn(icont); const bool canonicalizable = lead_coeff.is_integer(); const bool unit_normal = lead_coeff.is_pos_integer(); - - if (icont != *_num1_p) { - return (new mul(power(ebasis/icont, *num_exponent), power(icont, *num_exponent)) - )->setflag(status_flags::dynallocated); - } - - if (canonicalizable && (! unit_normal)) { - if (num_exponent->is_even()) { - return power(-ebasis, *num_exponent); - } else { - return (new mul(power(-ebasis, *num_exponent), *_num_1_p) - )->setflag(status_flags::dynallocated); - } + if (canonicalizable && (! unit_normal)) + icont = icont.mul(*_num_1_p); + + if (canonicalizable && (icont != *_num1_p)) { + const add& addref = ex_to<add>(ebasis); + add* addp = new add(addref); + addp->setflag(status_flags::dynallocated); + addp->clearflag(status_flags::hash_calculated); + addp->overall_coeff = ex_to<numeric>(addp->overall_coeff).div_dyn(icont); + for (epvector::iterator i = addp->seq.begin(); i != addp->seq.end(); ++i) + i->coeff = ex_to<numeric>(i->coeff).div_dyn(icont); + + const numeric c = icont.power(*num_exponent); + if (likely(c != *_num1_p)) + return (new mul(power(*addp, *num_exponent), c))->setflag(status_flags::dynallocated); + else + return power(*addp, *num_exponent); } } -- 1.4.4.4 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH] fast check for expanded expressions: introduce info_flags::expanded.
by Alexei Sheplyakov 17 Aug '07

17 Aug '07

Dear all, Sometimes it would be nice to check if the expression is already expanded. This information is cached by GiNaC, but the user code had no access to it. --- ginac/expairseq.cpp | 2 ++ ginac/flags.h | 5 ++++- ginac/numeric.cpp | 1 + ginac/power.cpp | 2 ++ ginac/symbol.cpp | 3 ++- 5 files changed, 11 insertions(+), 2 deletions(-) diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 4dffafa..757f3c3 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -267,6 +267,8 @@ void expairseq::do_print_tree(const print_tree & c, unsigned level) const bool expairseq::info(unsigned inf) const { + if (inf == info_flags::expanded) + return (flags & status_flags::expanded); return inherited::info(inf); } diff --git a/ginac/flags.h b/ginac/flags.h index a416958..86d8238 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -250,7 +250,10 @@ public: has_indices, // object has at least one index // answered by class idx - idx + idx, + + // answered by classes numeric, symbol, add, mul, power + expanded }; }; diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 1858cbd..ef61ce7 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -576,6 +576,7 @@ bool numeric::info(unsigned inf) const case info_flags::numeric: case info_flags::polynomial: case info_flags::rational_function: + case info_flags::expanded: return true; case info_flags::real: return is_real(); diff --git a/ginac/power.cpp b/ginac/power.cpp index 8829bba..faaef18 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -231,6 +231,8 @@ bool power::info(unsigned inf) const case info_flags::algebraic: return !exponent.info(info_flags::integer) || basis.info(inf); + case info_flags::expanded: + return (flags & status_flags::expanded); } return inherited::info(inf); } diff --git a/ginac/symbol.cpp b/ginac/symbol.cpp index a84ae4a..9f3b1c6 100644 --- a/ginac/symbol.cpp +++ b/ginac/symbol.cpp @@ -214,7 +214,8 @@ bool symbol::info(unsigned inf) const inf == info_flags::cinteger_polynomial || inf == info_flags::rational_polynomial || inf == info_flags::crational_polynomial || - inf == info_flags::rational_function) + inf == info_flags::rational_function || + inf == info_flags::expanded) return true; if (inf == info_flags::real) return domain==domain::real || domain==domain::positive; -- 1.4.4.4 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH] {add, power}::eval(): convert terms into primitive polynomials when possible.
by Alexei Sheplyakov 14 Aug '07

14 Aug '07

Hello! The following transformations are done automagically now: (-x+y)*(x-y) -> -(x-y)^2 (3*x+y)/(6*x+2*y)^2 -> 1/4*(y+3*x)^(-1) Please note: exam_series13() from check/exam_pseries.cpp fails due to this patch, because the expression 1/x*(-(1+x)/(1-x)) + (1+x)/x/(1-x) gets evaluated to zero. --- ginac/mul.cpp | 55 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- ginac/power.cpp | 22 ++++++++++++++++++++++ 2 files changed, 76 insertions(+), 1 deletions(-) diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 9bbcd61..4902cf1 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -448,7 +448,60 @@ ex mul::eval(int level) const ex_to<numeric>(addref.overall_coeff). mul_dyn(ex_to<numeric>(overall_coeff)))) ->setflag(status_flags::dynallocated | status_flags::evaluated); - } + } else if (seq_size >= 2) { + // Strip the content and the unit part from the each term. Thus + // things like (-x+a)*(3*x-3*a) automagically turn into - 3*(x-a)^2 + + epvector::const_iterator last = seq.end(); + epvector::const_iterator i = seq.begin(); + while (i!=last) { + if (! (is_a<add>(i->rest) && i->coeff.is_equal(_ex1))) { + // power::eval has such a rule, no need to handle powers here + ++i; + continue; + } + + // XXX: What is the best way to check if the polynomial is a primitive? + numeric c = i->rest.integer_content(); + const numeric& lead_coeff = + ex_to<numeric>(ex_to<add>(i->rest).seq.begin()->coeff).div_dyn(c); + const bool canonicalizable = lead_coeff.is_integer(); + + // XXX: The main variable is chosen in a random way, so this code + // does NOT transform the term into the canonical form (thus, in some + // very unlucky event it can even loop forever). Hopefully the main + // variable will be the same for all terms in *this + const bool unit_normal = lead_coeff.is_pos_integer(); + if ((c == *_num1_p) && ((! canonicalizable) || unit_normal)) { + ++i; + continue; + } + + std::auto_ptr<epvector> s(new epvector); + s->reserve(seq.size()); + + epvector::const_iterator j=seq.begin(); + while (j!=i) { + s->push_back(*j); + ++j; + } + + if (! unit_normal) + c = c.mul(*_num_1_p); + const ex primitive = (i->rest)/c; + s->push_back(expair(primitive, _ex1)); + ++j; + + while (j!=last) { + s->push_back(*j); + ++j; + } + return (new mul(s, + ex_to<numeric>(overall_coeff).mul_dyn(c)) + )->setflag(status_flags::dynallocated); + } + } + return this->hold(); } diff --git a/ginac/power.cpp b/ginac/power.cpp index fb5e351..2b8bfc5 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -477,6 +477,28 @@ ex power::eval(int level) const return expand_mul(ex_to<mul>(ebasis), *num_exponent, 0); } + // (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4) + if (num_exponent->is_integer() && is_exactly_a<add>(ebasis)) { + const numeric icont = ebasis.integer_content(); + const numeric& lead_coeff = + ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div_dyn(icont); + + const bool canonicalizable = lead_coeff.is_integer(); + const bool unit_normal = lead_coeff.is_pos_integer(); + + if (icont != *_num1_p) + return (new mul(power(ebasis/icont, *num_exponent), + power(icont, *num_exponent)))->setflag(status_flags::dynallocated); + + if (canonicalizable && (! unit_normal)) { + if (num_exponent->is_even()) + return power(-ebasis, *num_exponent); + else + return (new mul(power(-ebasis, *num_exponent), + *_num_1_p))->setflag(status_flags::dynallocated); + } + } + // ^(*(...,x;c1),c2) -> *(^(*(...,x;1),c2),c1^c2) (c1, c2 numeric(), c1>0) // ^(*(...,x;c1),c2) -> *(^(*(...,x;-1),c2),(-c1)^c2) (c1, c2 numeric(), c1<0) if (is_exactly_a<mul>(ebasis)) { -- 1.4.4.4 Best regards, Alexei -- All science is either physics or stamp collecting.

4 12

small clifford.cpp patch
by Vladimir Kisil 01 Aug '07

01 Aug '07

Dear all, Here is a few lines patch to clifford.cpp which makes clifford_moebius_map() function a little bit more intelligent. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv(a)maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/

2 1

Clifford patch from May 25th
by Jens Vollinga 13 Jul '07

13 Jul '07

Hi, I put Vladimir's patch from May 25th into CVS. Now Clifford units honor the representation labels. ... and a belated thanks for the patch! Regards, Jens

1 0

patch: valid C code for idx class with print_csrc
by Markus Grabner 11 Jul '07

11 Jul '07

Hi! Printing an indexed expression using an instance of the print_csrc_double class produces invalid C code (the same syntax that is used in the tutorial http://www.ginac.de/tutorial/Indexed-objects.html#Indexed-objects) in ginac-1.3.7. I understand that the index class in ginac is more powerful than the indexing facilities of the C language, but still I expected print_csrc_double to produce valid C code simply due to its name. So I added a few lines which handle the case of integer indices (numeric or symbolic) stored in the idx class (patch attached). Kind regards, Markus -- Markus Grabner - Computer Graphics and Vision Graz University of Technology, Inffeldgasse 16a/II, 8010 Graz, Austria Phone: +43/316/873-5041, Fax: +43/316/873-5050 WWW: http://www.icg.tugraz.at/Members/grabner

3 4

Re: Clifford base names WAS: Re: [GiNaC-list] How can I get a GiNaC::ex's name
by Vladimir Kisil 25 May '07

25 May '07

Dear Javier, >>>>> "JRG" == Javier Ros Ganuza <jros(a)unavarra.es> writes: JRG> It woul be interesting to have this posibility, tough elements JRG> of diferent basis are right now plotted in the same way. Indeed it may be worth to see a difference in the output for Clifford units with different representation labels. I include a patch which give this in the way which is shown by such an example: #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main() { try { varidx mu(symbol("nu", "\\nu"), 3); ex basis1 = clifford_unit(mu, diag_matrix(lst(1, 1, 1)),0), basis2 = clifford_unit(mu, diag_matrix(lst(1, 1, 1)),1); ex e00 = basis1.subs(mu == 0), e10 = basis2.subs(mu == 0); cout << e00 << endl; // <-- e~0 cout << e10 << endl; // <-- e[1]~0 } catch (exception &p) { cerr << "Got problem: " << p.what() << endl; } } Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv(a)maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/

1 0

Symbolic Integration
by Remco Bloemen 16 May '07

16 May '07

Dear developers, I have been using the Ginac library for some of my hobby programs now and I believe a generic C++ CAS would be very valuable. I have grown interest in symbolic integration and tried implementing polynomial factorization, it seems quite hard to do efficiently, but doable. I have even bought the book mentioned in the todo list (the second edition). Altough I have just began to read it the book mentions that some recent algorithms for trancedental integration do not require polynomial factorization. So I would like to try and implement one of these algo's. My question is, what are your thoughts on how to implement integration in a nice fashion? I tought about implementing some kind of antiderivative operator wich can be used in expressions and gives the indefinite integral. The reason I preffer operators to member function is that it allows unsolvable integrals to remain in the equation and still be numerically evaluable. In this fashion I would like to create a differential operator so ODE's and PDE's can be expressed in Ginac and numerically solved. The hobby program I am trying to make should become a generic numerical integrator for arbitrary PDE's by dynamicaly generating code. I already have a class that converts an expression to its Chebyshev-Pade approximation, compiles it to double arithmetic and dynamicly links it in. So gentleman, what are your thoughts? Sincerely, R. Bloemen

2 1

Problem with power series
by Stefan Weinzierl 16 May '07

16 May '07

Gentlemen, I came across a small problem, when calculating series expansions. The following code hangs, when attempting to compute the expansion: #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main() { symbol x("x"); ex e1 = pow(2,x)*( 1/x*(-(1+x)/(1-x)) + (1+x)/x/(1-x)); cout << e1 << endl; cout << e1.expand() << endl; cout << e1.series(x,1) << endl; } I tried it with GiNaC 1.3.7 as well as 1.3.1 and 1.3.0. e1 is actually a complicated zero, and the correct answer should we O(x). Best wishes, Stefan

3 3

Optimized C-style output
by Stefan Weinzierl 06 Apr '07

06 Apr '07

Hello, recently I was generating C-style output from GiNaC expressions and ran into the following problems: - the compiled code was inefficient, as common subexpressions are calculated over and over again, - the compiler didn't really appreciate single lines with about one million or more characters. On the mailing list I found that the issue of long lines was already discussed in July 2003, but it seems that this discussion did not lead to additional functionality of GiNaC. To overcome these problems, I wrote a routine which generates optimized C code. This routine - substitutes every subexpression, which occurs more than once with a temporary variable. - splits long expressions into several lines. As this could be of use to other people as well, it could be included into GiNaC. The source files are attached to this mail. An example for this routine would be the following code fragment symbol x("x"), y("y"); ex f = pow(sin(x+y),2)+sin(x+y)+log(sin(x+y)); optimized_C_code(std::cout, std::string("double"), std::string("tmp"), f); which prints double tmp1 = sin(y+x); double tmp0 = tmp1+(tmp1*tmp1)+log(tmp1); to standard output. Best wishes, Stefan

5 14