August 2007 - GiNaC-devel

Volunteer to continue the work on NTL factorisation bindings
by Remco Bloemen 06 Sep '07

06 Sep '07

I have created a simple factorisation for polynomials in Q[x] with real roots (source attached, based on [2]) but it's really slow. Creating a reasonable implementation would require much work. An NTL based implementation would be very powerfull and feasible, as said in [1]. I would like to volunteer to continue this work. I have the following questions: 1) Where can I obtain the source of this binding? (or should I start from the code in gTybalt, or from scratch?) 2) Is there any chance such a patch will be included? (this would imply dependecy on NTL) 2a) And if it is compile-time optional? 2b) And if it also included integration of rational functions? 2c) And if it also included integration of trancedental functions? 3) Do you have any related advice, hints or pointers? [1] http://www.ginac.de/pipermail/ginac-devel/2004-August/000695.html [2] http://yacas.sourceforge.net/Algochapter3.html

3 3

GiNaC release 1.4.0
by Jens Vollinga 04 Sep '07

04 Sep '07

Hello everybody, GiNaC 1.4.0 is out and available. The changes are: - New tinfo mechanism. - Removed rtt from class symbol. - Faster archiving by adding a map from strings to indices in the atoms vector. - Improved unarchiving: algorithms order N instead of order N^2. - Added compile function, excompiler class. - Added exset type. - Added step function to GiNaCs built-in functions. - Added is_polynomial() method. - Added real_part() and imag_part() methods. - Added matrix::is_zero_matrix() method. - Added evalm() method to pseries class. - Improved dummy index renaming. - Many improvements for class clifford. - New flag for positive symbols/numbers, added possymbol class. - Added programming examples in doc/examples *Note*: the tinfo mechanism has changed from fixed magic numbers to dynamically assigned tinfo-ids. If have written your own GiNaC classes, you need to adjust your code! As always, this release can be downloaded from http://www.ginac.de/Download.html Best wishes, Jens

2 5

GiNaC release 1.3.8
by Jens Vollinga 28 Aug '07

28 Aug '07

Hello everybody, GiNaC 1.3.8 is out and available. The changes are: - Drop support of ginac-config and ginac.m4. Please use pkg-config instead. - atan2(y,x) branch cut correction. - Fixed bug in series expansion. - Additional transformations for mul and power. - Clifford units now honor representation labels. - Correct csrc output for idx and fderivative. - Improved handling of convergence transformationis for Li/G. - Fix compilation issues with prereleases of GCC 4.3. - info_flags::expanded added. - input_parser also accepts pow(). - texinfo.tex removed from package. As always, this release can be downloaded from http://www.ginac.de/Download.html Best wishes, Jens

1 0

Release?
by Vladimir Kisil 28 Aug '07

28 Aug '07

Jens Vollinga <vollinga at physik.uni-wuppertal.de> Mon Aug 27 10:59:40 CEST 2007: > I'd like to do a 1.3.8 release today and (probably) a 1.4 release > tomorrow. If somebody has some patches in the making or wants certain > issues/problems to be resolved first, I will delay that plan, of course. Is file ginac/compiler.h from the Alexei's patch in the tree? Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv(a)maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/

2 1

Release?
by Jens Vollinga 27 Aug '07

27 Aug '07

Hi, I'd like to do a 1.3.8 release today and (probably) a 1.4 release tomorrow. If somebody has some patches in the making or wants certain issues/problems to be resolved first, I will delay that plan, of course. Regards, Jens

1 0

[PATCH] Change input parser to accept pow(x, y) and power(x, y).
by Alexei Sheplyakov 27 Aug '07

27 Aug '07

so symbol x("x"); lst syms(x); ex("pow(x, 2)", syms); does NOT trigger parser error any more. --- ginac/input_parser.yy | 4 ++++ 1 files changed, 4 insertions(+), 0 deletions(-) diff --git a/ginac/input_parser.yy b/ginac/input_parser.yy index 8601524..74beb4c 100644 --- a/ginac/input_parser.yy +++ b/ginac/input_parser.yy @@ -105,6 +105,10 @@ exp : T_NUMBER {$$ = $1;} if ($3.nops() != 1) throw (std::runtime_error("too many arguments to sqrt()")); $$ = sqrt($3.op(0)); + } else if (n == "pow" || n == "power") { + if ($3.nops() != 2) + throw std::invalid_argument("wrong number of arguments to pow()"); + $$ = power($3.op(0), $3.op(0)); } else { unsigned i = function::find_function(n, $3.nops()); $$ = function(i, ex_to<exprseq>($3)).eval(1); -- 1.4.4.4 Best regards, Alexei -- All science is either physics or stamp collecting.

2 1

lst_to_clifford() and clifford_to_lst() patch
by Vladimir Kisil 27 Aug '07

27 Aug '07

Dear All, I wish to propose a small patch in clifford.cpp which will enable functions lst_to_clifford() and clifford_to_lst() to use pseudo-vector representation along with the vector one (which was implemented before). The new behaviour is reflected in the tutorial and exam_clifford.cpp. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv(a)maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/ Index: doc/tutorial/ginac.texi =================================================================== RCS file: /home/cvs/GiNaC/doc/tutorial/ginac.texi,v retrieving revision 1.197 diff -u -r1.197 ginac.texi --- doc/tutorial/ginac.texi 2 Aug 2007 14:56:44 -0000 1.197 +++ doc/tutorial/ginac.texi 13 Aug 2007 13:07:05 -0000 @@ -3401,8 +3401,20 @@ with @samp{e.k} directly supplied in the second form of the procedure. In the first form the Clifford unit @samp{e.k} is generated by the call of -@code{clifford_unit(mu, metr, rl)}. The previous code may be rewritten -with the help of @code{lst_to_clifford()} as follows +@code{clifford_unit(mu, metr, rl)}. +@cindex pseudo-vector +If the number of components supplied +by @code{v} exceeds the dimensionality of the Clifford unit @code{e} by +1 then function @code{lst_to_clifford()} uses the following +pseudo-vector representation: +@tex +$v^0 {\bf 1} + v^1 e_0 + v^2 e_1 + ... + v^{n+1} e_n$ +@end tex +@ifnottex +@samp{v~0 ONE + v~1 e.0 + v~2 e.1 + ... + v~[n+1] e.n} +@end ifnottex + +The previous code may be rewritten with the help of @code{lst_to_clifford()} as follows @example @{ @@ -3432,17 +3444,25 @@ @ifnottex @samp{v = (v~0, v~1, ..., v~n)} @end ifnottex -such that +such that the expression is either vector @tex $e = v^0 c_0 + v^1 c_1 + ... + v^n c_n$ @end tex @ifnottex @samp{e = v~0 c.0 + v~1 c.1 + ... + v~n c.n} @end ifnottex -with respect to the given Clifford units @code{c} and with none of the -@samp{v~k} containing Clifford units @code{c} (of course, this +or pseudo-vector +@tex +$v^0 {\bf 1} + v^1 e_0 + v^2 e_1 + ... + v^{n+1} e_n$ +@end tex +@ifnottex +@samp{v~0 ONE + v~1 e.0 + v~2 e.1 + ... + v~[n+1] e.n} +@end ifnottex +with respect to the given Clifford units @code{c}. Here none of the +@samp{v~k} should contain Clifford units @code{c} (of course, this may be impossible). This function can use an @code{algebraic} method -(default) or a symbolic one. With the @code{algebraic} method the @samp{v~k} are calculated as +(default) or a symbolic one. With the @code{algebraic} method the +@samp{v~k} are calculated as @tex $(e c_k + c_k e)/c_k^2$. If $c_k^2$ @end tex Index: ginac/clifford.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/clifford.cpp,v retrieving revision 1.102 diff -u -r1.102 clifford.cpp --- ginac/clifford.cpp 1 Aug 2007 14:27:25 -0000 1.102 +++ ginac/clifford.cpp 13 Aug 2007 13:07:05 -0000 @@ -1183,13 +1183,20 @@ if (min == 1) { if (dim == max) return indexed(v, mu_toggle) * e; - else + else if (max - dim == 1) { + if (ex_to<matrix>(v).cols() > ex_to<matrix>(v).rows()) + return v.op(0) * dirac_ONE(ex_to<clifford>(e).get_representation_label()) + indexed(sub_matrix(ex_to<matrix>(v), 0, 1, 1, dim), mu_toggle) * e; + else + return v.op(0) * dirac_ONE(ex_to<clifford>(e).get_representation_label()) + indexed(sub_matrix(ex_to<matrix>(v), 1, dim, 0, 1), mu_toggle) * e; + } else throw(std::invalid_argument("lst_to_clifford(): dimensions of vector and clifford unit mismatch")); } else throw(std::invalid_argument("lst_to_clifford(): first argument should be a vector (nx1 or 1xn matrix)")); } else if (v.info(info_flags::list)) { if (dim == ex_to<lst>(v).nops()) return indexed(matrix(dim, 1, ex_to<lst>(v)), mu_toggle) * e; + else if (ex_to<lst>(v).nops() - dim == 1) + return v.op(0) * dirac_ONE(ex_to<clifford>(e).get_representation_label()) + indexed(sub_matrix(matrix(dim+1, 1, ex_to<lst>(v)), 1, dim, 0, 1), mu_toggle) * e; else throw(std::invalid_argument("lst_to_clifford(): list length and dimension of clifford unit mismatch")); } else @@ -1274,19 +1281,28 @@ or (not is_a<numeric>(pow(c.subs(mu == i, subs_options::no_pattern), 2)))) algebraic = false; lst V; + ex v0 = remove_dirac_ONE(canonicalize_clifford(e+clifford_prime(e)).normal())/2; + if (not v0.is_zero()) + V.append(v0); + ex e1 = canonicalize_clifford(e - v0 * dirac_ONE(ex_to<clifford>(c).get_representation_label())); if (algebraic) { for (unsigned int i = 0; i < D; i++) V.append(remove_dirac_ONE( - simplify_indexed(canonicalize_clifford(e * c.subs(mu == i, subs_options::no_pattern) + c.subs(mu == i, subs_options::no_pattern) * e)) + simplify_indexed(canonicalize_clifford(e1 * c.subs(mu == i, subs_options::no_pattern) + c.subs(mu == i, subs_options::no_pattern) * e1)) / (2*pow(c.subs(mu == i, subs_options::no_pattern), 2)))); } else { - ex e1 = canonicalize_clifford(e); try { for (unsigned int i = 0; i < D; i++) V.append(get_clifford_comp(e1, c.subs(c.op(1) == i, subs_options::no_pattern))); } catch (std::exception &p) { /* Try to expand dummy summations to simplify the expression*/ e1 = canonicalize_clifford(expand_dummy_sum(e1, true)); + V.remove_all(); + v0 = remove_dirac_ONE(canonicalize_clifford(e1+clifford_prime(e1)).normal())/2; + if (not v0.is_zero()) { + V.append(v0); + e1 = canonicalize_clifford(e1 - v0 * dirac_ONE(ex_to<clifford>(c).get_representation_label())); + } for (unsigned int i = 0; i < D; i++) V.append(get_clifford_comp(e1, c.subs(c.op(1) == i, subs_options::no_pattern))); } Index: check/exam_clifford.cpp =================================================================== RCS file: /home/cvs/GiNaC/check/exam_clifford.cpp,v retrieving revision 1.36 diff -u -r1.36 exam_clifford.cpp --- check/exam_clifford.cpp 31 Jan 2007 22:29:19 -0000 1.36 +++ check/exam_clifford.cpp 13 Aug 2007 13:07:05 -0000 @@ -410,13 +410,23 @@ result += check_equal(canonicalize_clifford(e), 0); /* lst_to_clifford() and clifford_inverse() check*/ - realsymbol x("x"), y("y"), t("t"), z("z"); + realsymbol s("s"), t("t"), x("x"), y("y"), z("z"); ex c = clifford_unit(nu, A, 1); e = lst_to_clifford(lst(t, x, y, z), mu, A, 1) * lst_to_clifford(lst(1, 2, 3, 4), c); e1 = clifford_inverse(e); result += check_equal_simplify_term2((e*e1).simplify_indexed(), dirac_ONE(1)); +/* lst_to_clifford() and clifford_to_lst() check for vectors*/ + e = lst(t, x, y, z); + result += check_equal_lst(clifford_to_lst(lst_to_clifford(e, c), c, false), e); + result += check_equal_lst(clifford_to_lst(lst_to_clifford(e, c), c, true), e); + +/* lst_to_clifford() and clifford_to_lst() check for pseudovectors*/ + e = lst(s, t, x, y, z); + result += check_equal_lst(clifford_to_lst(lst_to_clifford(e, c), c, false), e); + result += check_equal_lst(clifford_to_lst(lst_to_clifford(e, c), c, true), e); + /* Moebius map (both forms) checks for symmetric metrics only */ matrix M1(2, 2), M2(2, 2); c = clifford_unit(nu, A);

2 2

Please remove texinfo.tex from CVS.
by varg＠theor.jinr.ru 22 Aug '07

22 Aug '07

Hello! GiNaC CVS holds an obsolete copy of texinfo.tex (doc/tutorial/texinfo.tex). Due to this file tutorial and reference fail to build with tetex >= 3.0 (and probably any other more-or-less up to date distribution of TeX). So, could anyone with write access purge that file, please? Storing such files in CVS (or any other SCM system) is just asking for troubles. Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

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

19 Aug '07

The same patch for GiNaC 1.3. Make 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: the numeric coefficient is calculated in one sweep, sums divieded by numbers in place. --- ginac/compiler.h | 12 ++++++++++++ ginac/mul.cpp | 48 ++++++++++++++++++++++++++++++++++-------------- ginac/power.cpp | 34 ++++++++++++++++++++-------------- 3 files changed, 66 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 230717a..3c87d5a 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 { @@ -414,7 +415,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; @@ -440,14 +441,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; @@ -464,32 +469,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 0dfbc8f..a94cc47 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -40,6 +40,8 @@ #include "lst.h" #include "archive.h" #include "utils.h" +#include "relational.h" +#include "compiler.h" namespace GiNaC { @@ -465,25 +467,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] {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