[PATCH] {add, power}::eval(): convert terms into primitive polynomials when possible.
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.
Almost the same patch for GiNaC 1.3, just in a case... The following transformation is done automagically now: (-x+y)*(x-y) -> -(x-y)^2 (3*x+y)/(6*x+2*y)^2 -> 1/4*(y+3*x)^(-1) --- 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 cbab7e7..8863949 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -440,7 +440,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 a63556b..0573dc6 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -463,6 +463,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.
Hi, Alexei Sheplyakov schrieb:
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.
thanks for the patch! I modified exam_series13() so that it doesn't use the new transformation. Regards, Jens
Hi Alexei! Alexei Sheplyakov wrote:
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.
This patch appears to come with a slight performance hit. On my machine, 14x14 Vandermonde matrix determinants degrades from 19s to 32s, Lewis-Wester test O1 from 7.5s to 11s, and Lewis-Wester test A from 0.055s to 0.22s. (I'm seriously surprised by the latter since it only involves integers, but it's reproducible.) Looking at the code it appears to me like mul::eval recurses now. Wouldn't it be faster to collect the coefficients in one sweep? Well, I recognize that most terms will be trivial in the following recursions, but still. Also, you write that "in some very unlucky event it can even loop forever". I don't quite understand. Do you have an example for this? Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi Alexei and Richy, Richard B. Kreckel schrieb:
Alexei Sheplyakov wrote:
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.
This patch appears to come with a slight performance hit. On my machine, 14x14 Vandermonde matrix determinants degrades from 19s to 32s, Lewis-Wester test O1 from 7.5s to 11s, and Lewis-Wester test A from 0.055s to 0.22s. (I'm seriously surprised by the latter since it only involves integers, but it's reproducible.)
Looking at the code it appears to me like mul::eval recurses now. Wouldn't it be faster to collect the coefficients in one sweep? Well, I recognize that most terms will be trivial in the following recursions, but still.
Also, you write that "in some very unlucky event it can even loop forever". I don't quite understand. Do you have an example for this?
something new here? I kind of feel like doing a release soon, 1.3.8 and also maybe then a 1.4, but this seems to be a sort of issue to be solved first. Regards, Jens
Hello! On Thu, Jul 12, 2007 at 11:23:02PM +0200, Richard B. Kreckel wrote:
Alexei Sheplyakov wrote:
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.
This patch appears to come with a slight performance hit.
Indeed, mul::eval() has more work to be done. advantages: - subsequent operations, in particular gcd, are faster. - no more ugly stuff like (-x+y)*(x-y). disadvantages: - subsequent operations, in particualr expand, may be slower. - there is some overhead in the trivial case, e.g. (x*y*z).eval().
On my machine, 14x14 Vandermonde matrix determinants degrades from 19s to 32s, Lewis-Wester test O1 from 7.5s to 11s,
On my box results are 14.847s versus 16.949s
and Lewis-Wester test A from 0.055s to 0.22s.
Ouch. I do observe some slowdown (from 0.292s to 0.3s) too, but it is not *that* bad.
(I'm seriously surprised by the latter since it only involves integers, but it's reproducible.)
a/b is actually mul(a, power(b, _ex_1)), and it has seq.size() >= 2. My patch adds empty loop over 2 elements, so there should be some slowdown.
Looking at the code it appears to me like mul::eval recurses now.
mul::eval calls evalchildren, so it is recursive even without my patch.
Wouldn't it be faster to collect the coefficients in one sweep?
I've tried this (see the patch [1]), it has almost no effect.
Well, I recognize that most terms will be trivial in the following recursions, but still.
My patch causes overhead even in the trivial case, e.g. when the polynomial is already unit normal. Consider e.g. e = (x+1)*(2*x^2+x+1)*...*(100*x^100+99*x^99+...+1); mul::eval will do 100 of useless checks (each of them is relatively expansive on its own). Any ideas how to avoid this? E.g. what is the best way to check if polynomial is already primitive (unit normal)?
Also, you write that "in some very unlucky event it can even loop forever". I don't quite understand.
The main variable(s) is (are) random (the first term of the add object). During re-evaluation expression ordering might change, so the recursion is not guaranteed to terminate.
Do you have an example for this?
No, this is purely theoretical. [1] [PATCH] {mul,power}::eval(): remove extra recursion (divide by integer in place). --- ginac/mul.cpp | 42 +++++++++++++++++++++++++++++++----------- ginac/power.cpp | 16 +++++++++++++--- ginac/compiler.h | 12 ++++++++++++ 3 files changed, 56 insertions(+), 14 deletions(-) diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 4ceeb79..f9f7be5 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 { @@ -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 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(); + 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((! i->coeff.is_equal(_ex1)) || (! is_a<add>(i->rest)))) { // 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) c = c.mul(*_num_1_p); - const ex primitive = (i->rest)/c; - s->push_back(expair(primitive, _ex1)); + + oc = oc.mul(c); + + 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)) + ex_to<numeric>(overall_coeff).mul_dyn(oc)) )->setflag(status_flags::dynallocated); } } diff --git a/ginac/power.cpp b/ginac/power.cpp index 238e12e..0103398 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -478,9 +478,19 @@ ex power::eval(int level) const 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 (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); + + return (new mul(power(*addp, *num_exponent), + icont.power(*num_exponent)))->setflag(status_flags::dynallocated); + } + if (canonicalizable && (! unit_normal)) { if (num_exponent->is_even()) 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 */ -- 1.4.4.4 Best regards, Alexei -- All science is either physics or stamp collecting.
Hi! Sheplyakov Alexei wrote:
This patch appears to come with a slight performance hit.
Indeed, mul::eval() has more work to be done.
advantages: - subsequent operations, in particular gcd, are faster. - no more ugly stuff like (-x+y)*(x-y).
disadvantages: - subsequent operations, in particualr expand, may be slower. - there is some overhead in the trivial case, e.g. (x*y*z).eval().
It is not some overhead that worries me. Rather, it is the algorithmic complexity of the new overhead. The new transformation appears to introduce quadratic behavior at the eval() level. That's bad. One of the points of choosing the algorithms that are done at the eval() stage to be only subquadratic algorithms was that otherwise one might be unable to merely construct the ex object.
and Lewis-Wester test A from 0.055s to 0.22s.
Ouch. I do observe some slowdown (from 0.292s to 0.3s) too, but it is not *that* bad.
(I'm seriously surprised by the latter since it only involves integers, but it's reproducible.)
a/b is actually mul(a, power(b, _ex_1)), and it has seq.size() >= 2. My patch adds empty loop over 2 elements, so there should be some slowdown.
Well, I still fail to see how that loop can make a measurable effect. But never mind this particular benchmark.
Looking at the code it appears to me like mul::eval recurses now.
mul::eval calls evalchildren, so it is recursive even without my patch.
Well, but that recurses into *children* only! At the end of the loop that you introduced into mul::eval(), a new mul object is returned. It is not evaluated yet and enters mul::eval() again. And it might contain most of the original terms. This is an entirely different quantity of recursion that is responsible for potentially quadratic behavior, IINM.
Wouldn't it be faster to collect the coefficients in one sweep?
I've tried this (see the patch [1]), it has almost no effect.
It still appears to be quadratic. Must it really be quadratic?
Well, I recognize that most terms will be trivial in the following recursions, but still.
My patch causes overhead even in the trivial case, e.g. when the polynomial is already unit normal. Consider e.g.
e = (x+1)*(2*x^2+x+1)*...*(100*x^100+99*x^99+...+1);
Indeed. This expression is well suited to study the asymptotic behavior: $ { echo -n "time((1+x)"; for i in $(seq 2 100); do echo -n "*(1+x"; for j in $(seq 2 $i); do echo -n "+${j}*x^${j}"; done; echo -n ")"; done; echo ");";} | ginsh-1.3.5 0.089994s $ { echo -n "time((1+x)"; for i in $(seq 2 100); do echo -n "*(1+x"; for j in $(seq 2 $i); do echo -n "+${j}*x^${j}"; done; echo -n ")"; done; echo ");";} | ginsh-head 0.169989s $ { echo -n "time((1+x)"; for i in $(seq 2 1000); do echo -n "*(1+x"; for j in $(seq 2 $i); do echo -n "+${j}*x^${j}"; done; echo -n ")"; done; echo ");";} | ginsh-1.3.5 11.4992s $ { echo -n "time((1+x)"; for i in $(seq 2 1000); do echo -n "*(1+x"; for j in $(seq 2 $i); do echo -n "+${j}*x^${j}"; done; echo -n ")"; done; echo ");";} | ginsh-head 68.2522s
mul::eval will do 100 of useless checks (each of them is relatively expansive on its own). Any ideas how to avoid this? E.g. what is the best way to check if polynomial is already primitive (unit normal)?
Wild idea: If the primitive predicate becomes really important, would it help to cache it like it's done for the status_flags::expanded predicate inside basic::flags? (Maybe using two bits to represent the tristate yes, no, unknown.) Feel free to flame me if this idea is bogus for some reason. Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
On Sat, Aug 11, 2007 at 11:11:57PM +0200, Richard B. Kreckel wrote:
It is not some overhead that worries me. Rather, it is the algorithmic complexity of the new overhead. The new transformation appears to introduce quadratic behavior at the eval() level.
[snipped]
Looking at the code it appears to me like mul::eval recurses now.
mul::eval calls evalchildren, so it is recursive even without my patch.
Well, but that recurses into *children* only! At the end of the loop that you introduced into mul::eval(), a new mul object is returned. It is not evaluated yet and enters mul::eval() again. And it might contain most of the original terms. This is an entirely different quantity of recursion that is responsible for potentially quadratic behavior, IINM.
mul::eval() *is quadratic* even without my patch: ex mul::eval(int level) const { std::auto_ptr<epvector> evaled_seqp = evalchildren(level); if (evaled_seqp.get()) { // do more evaluation later return (new mul(evaled_seqp, overall_coeff))-> setflag(status_flags::dynallocated); } If evaluation of children gives something non-trivial, a new mul object is returned. It is not evaluated and enters mul::eval() again.
That's bad.
My patch increases the coefficient. And that is indeed awful.
One of the points of choosing the algorithms that are done at the eval() stage to be only subquadratic algorithms was that otherwise one might be unable to merely construct the ex object.
I think this requirement is way too strict. I'd say anything of polynomial complexity (with "reasonable" coefficients) is OK.
Wouldn't it be faster to collect the coefficients in one sweep?
I've tried this (see the patch [1]), it has almost no effect.
It still appears to be quadratic. Must it really be quadratic?
Yes. After numeric coefficients are collected, some re-evaluation is necessary. For example, (3*x+3*y)*(x+y) -> (3*(x+y)*(x+y)).eval() -> 3*(x+y)^2 1/9*(9*x+9*y)*(x+y)^2 -> ((x+y)*(x+y)^2).eval() -> (x+y)^3 (N.B.: the result is not mul object)
Well, I recognize that most terms will be trivial in the following recursions, but still.
My patch causes overhead even in the trivial case, e.g. when the polynomial is already unit normal. Consider e.g.
e = (x+1)*(2*x^2+x+1)*...*(100*x^100+99*x^99+...+1);
Indeed. This expression is well suited to study the asymptotic behavior:
$ { echo -n "time((1+x)"; for i in $(seq 2 100); do echo -n "*(1+x"; for j in $(seq 2 $i); do echo -n "+${j}*x^${j}"; done; echo -n ")"; done; echo ");";} | ginsh-1.3.5 0.089994s
To make the test a little bit more interesting, let's take more than 2 points: { echo "plot '-' using 1:2 with lines title 'mul::eval() asymptotic behaviour'" echo for N in `seq 100 100 2000`; do TIME=$({ echo -n "time((1+x)" for i in $(seq 2 $N); do echo -n "*(1+x" for j in $(seq 2 $i); do echo -n "+${j}*x^${j}" done echo -n ")" done echo ");" ; } | ginsh | sed -e 's/s//g') echo "$N $TIME" done echo "e" ; } | gnuplot -persist '-' This gives a beautiful parabola (both with my patch or without it). Coefficients are different, but both versions of mul::eval() are quadratic. Best regards, Alexei -- All science is either physics or stamp collecting.
Hi! Sheplyakov Alexei wrote:
On Sat, Aug 11, 2007 at 11:11:57PM +0200, Richard B. Kreckel wrote:
One of the points of choosing the algorithms that are done at the eval() stage to be only subquadratic algorithms was that otherwise one might be unable to merely construct the ex object.
I think this requirement is way too strict. I'd say anything of polynomial complexity (with "reasonable" coefficients) is OK.
Not for add objects! And we really had evaluation of sums in mind when we said "no quadratic algorithms", well knowing that matrix::eval may be much worse without anybody realistically running into trouble (but that depends on how you count, in a sense). Your data shows that mul::eval() has never been subquadratic and nobody has noticed so far. This suggests that different criteria and benchmarks apply for varying kinds of objects.
My patch causes overhead even in the trivial case, e.g. when the polynomial is already unit normal. Consider e.g.
e = (x+1)*(2*x^2+x+1)*...*(100*x^100+99*x^99+...+1); [...] This gives a beautiful parabola (both with my patch or without it). Coefficients are different, but both versions of mul::eval() are quadratic.
This is interesting. Taking logarithms of both axis it appears like the behavior is indeed O(n^2.5) in both cases. I don't know what others think. For my own part, you've dissolved my worries. Jens wants to roll a release. I guess he would be interested whether your last patch cuts it: is it worth applying or not? Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hello! On Sun, Aug 12, 2007 at 11:30:57PM +0200, Richard B. Kreckel wrote:
Not for add objects! And we really had evaluation of sums in mind when we said "no quadratic algorithms",
Unfortunately, add::eval() is not subquadratic either. For example, 1 + \sum_{i=1}^N (i x_i + (i+1) x_i^i x_{i+1}^{i+1} + (i+2) x_i^i x_{i+2}^{i+2}) #!/bin/sh { echo "plot '-' using 1:2 with lines title 'add::eval() asymptotic behaviour'" echo for N in `seq 100 100 3000`; do TIME=$({ echo -n "time(1" for i in $(seq 1 $N); do j=$(expr $i + 1) k=$(expr $i + 2) printf "+${i}*x${i}+${j}*x${i}^${i}*x${j}^${j}+${k}*x${i}^${i}*x${k}^${k}" done echo ");" ; } | ginsh | sed -e 's/s//g') echo "$N $TIME" done echo "e" ; } | gnuplot -persist '-' Best regards, Alexei -- All science is either physics or stamp collecting.
Hi, Sheplyakov Alexei wrote:
On Sun, Aug 12, 2007 at 11:30:57PM +0200, Richard B. Kreckel wrote:
Not for add objects! And we really had evaluation of sums in mind when we said "no quadratic algorithms",
Unfortunately, add::eval() is not subquadratic either. For example,
1 + \sum_{i=1}^N (i x_i + (i+1) x_i^i x_{i+1}^{i+1} + (i+2) x_i^i x_{i+2}^{i+2})
#!/bin/sh { echo "plot '-' using 1:2 with lines title 'add::eval() asymptotic behaviour'" echo for N in `seq 100 100 3000`; do TIME=$({ echo -n "time(1" for i in $(seq 1 $N); do j=$(expr $i + 1) k=$(expr $i + 2) printf "+${i}*x${i}+${j}*x${i}^${i}*x${j}^${j}+${k}*x${i}^${i}*x${k}^${k}" done echo ");" ; } | ginsh | sed -e 's/s//g') echo "$N $TIME" done echo "e" ; } | gnuplot -persist '-'
Experimentally, this is O(N^2.2) on my machine. And I suppose it is clear by now that this is due to the fact that the problem is really quadratic in input size -- as was the original problem where the top-level object was a mul instead of an add. So, I dare say that O(N^2.x) is actually pretty good for a problem that has input size O(N^2), right? Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi! On Mon, Aug 13, 2007 at 10:44:59PM +0200, Richard B. Kreckel wrote:
And I suppose it is clear by now that this is due to the fact that the problem is really quadratic in input size -- as was the original problem where the top-level object was a mul instead of an add.
The sum has 3N terms, and the product has N terms. Why the problems are quadratic? Best regards, Alexei -- All science is either physics or stamp collecting.
Hi Alexei! Sheplyakov Alexei wrote:
On Mon, Aug 13, 2007 at 10:44:59PM +0200, Richard B. Kreckel wrote:
And I suppose it is clear by now that this is due to the fact that the problem is really quadratic in input size -- as was the original problem where the top-level object was a mul instead of an add.
The sum has 3N terms, and the product has N terms. Why the problems are quadratic?
Because of the children terms: { echo "plot '-' using 1:2 with lines title 'characters in term to be evaluated'" echo for N in `seq 100 100 2000`; do SIZE=$({ echo -n "(0001+0.0001*x)" for i in $(seq 2 $N); do echo -n "*(0001+0001*x" for j in $(seq 2 $i); do printf "+%04d*x^%04d" ${j} ${i} done echo -n ")" done echo ";" ; } | wc -c) echo "$N $SIZE" done echo "e" ; } | gnuplot -persist '-' :) Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
participants (4)
-
Alexei Sheplyakov
-
Jens Vollinga
-
Richard B. Kreckel
-
varg＠theor.jinr.ru