Hi, here is a very odd performance behaviour of ginac/series: The operation ex e = 1/(a+sqrt(b*b+x*x)); e = e.series(x==1, 21); is more than a factor of 10^4(!) faster then outputting the result using cout << e << endl; I have measured these times using clocks(). You find the source of the program at the end of this email. What does this mean? This means, that series is very efficient (both time- and space-wise), because it reuses intermediate results heavily. On the other hand the result is nearly useless, because it is actually _very_ complicated, which shows up in any subsequent usage. In this specific example, a hand-written taylor-series using derivatives was 2 orders of magnitude slower (for 60th order) to compute the series, but the series (20th order) still was printed 1000 times faster than the standard series. In the light of this finding, something needs to be done. I suggest to essentially revert the logic of series, namely to compute series by differentiation and to fall back to an explicit composition if it does not work. Mathematica does it that way according to its documentation. At least ginac should provide this mode as an option, either by an additional taylor() method for ex or by means of series_options. Do You agree with this analysis and my conclusions? Are there any other ways to solve this problem? Regards Frieder Kalisch #include <ginac/ginac.h> #include <iostream> #include <ctime> using namespace GiNaC; using namespace std; int main(int argc, char **argv) { time_t t0, t1, t2; t0 = clock(); symbol x("x"), a("a"), b("b"); ex e = 1/(a+sqrt(b*b+x*x)); e = e.series(x==1, 21); t1 = clock(); cout << e << endl; t2 = clock(); double res = CLOCKS_PER_SEC; cerr << (t1-t0)/res << endl << (t2-t1)/res << endl; return 0; } // int main(int argc, char **argv) -- Frieder Kalisch Institut für theoretische Physik kalisch@tphys.uni-heidelberg.de Philosophenweg 19 +49-6221-549-433 D-69120 Heidelberg
Dear Frieder, On Wed, 19 Jan 2005, Frieder Kalisch wrote:
here is a very odd performance behaviour of ginac/series.
Okay, I had a quick look at this. There are two questions that I am asking myself now: (1) How should .series() perform its calculation? (2) What should the result be like? Ad (1): I think series expansion should most of the time not be done by repeated differentiation. What if some large addition appears as a subexpression? The chain rule together with expansion of the derivatives (in the sense of .expand()) will cause expression explosion and happyness is not going to occur. Ad (2): I think that the expansion coefficients should be in expandend form (again in the sense of .expand()). If the main developpers agree on this, then you have just located a bug. The coefficients of the result of your example look more than a bit unexpanded.... Best wishes, Chris
Hi Chris, Am Mittwoch 19 Januar 2005 18:31 schrieb Chris Dams:
Ad (1): I think series expansion should most of the time not be done by repeated differentiation. What if some large addition appears as a subexpression? The chain rule together with expansion of the derivatives (in the sense of .expand()) will cause expression explosion and happyness is not going to occur.
I admit that series expansion by repeated differentiation is not always the right solution. However, given the patch at the end of this mail, which cures the unexpandedness of power_const, series is still slower by a factor of 10 than series expansion by differentiation. Maybe there just isn't one single best method for series expansion... Therefore I think that ginac should provide a way to compute series by differientiation. I suggest adding a method to ex, which is essentially a call to basic::series. Given Your reasoning, I don't understand the implementation of power::series anymore. E.g. pow(sin(something very complicated without singularity), 5) is expanded by repeated differentiation. basic::series is called in the if statement just below "// Is the expression of type something^(-int)?".
Ad (2): I think that the expansion coefficients should be in expandend form (again in the sense of .expand()). If the main developpers agree on this, then you have just located a bug. The coefficients of the result of your example look more than a bit unexpanded....
The following patch makes pseries::power_const and pseries::mul_series return expanded series. Regards Frieder Kalisch *** GiNaC-1.3.0/ginac/pseries.cpp.save 2005-01-18 10:13:46.000000000 +0100 --- GiNaC-1.3.0/ginac/pseries.cpp 2005-01-20 12:34:28.000000000 +0100 *************** *** 756,762 **** ex a_coeff = coeff(var, i); ex b_coeff = other.coeff(var, cdeg-i); if (!is_order_function(a_coeff) && !is_order_function(b_coeff)) ! co += a_coeff * b_coeff; } if (!co.is_zero()) new_seq.push_back(expair(co, numeric(cdeg))); --- 756,762 ---- ex a_coeff = coeff(var, i); ex b_coeff = other.coeff(var, cdeg-i); if (!is_order_function(a_coeff) && !is_order_function(b_coeff)) ! co += (a_coeff * b_coeff).expand(); } if (!co.is_zero()) new_seq.push_back(expair(co, numeric(cdeg))); *************** *** 893,906 **** for (int i=1; i<deg; ++i) { ex sum = _ex0; for (int j=1; j<=i; ++j) { ! ex c = coeff(var, j + ldeg); if (is_order_function(c)) { co.push_back(Order(_ex1)); break; } else ! sum += (p * j - (i - j)) * co[i - j] * c; } ! co.push_back(sum / coeff(var, ldeg) / i); } // Construct new series (of non-zero coefficients) --- 893,906 ---- for (int i=1; i<deg; ++i) { ex sum = _ex0; for (int j=1; j<=i; ++j) { ! ex c = (coeff(var, j + ldeg) / coeff(var, ldeg) / i).expand(); if (is_order_function(c)) { co.push_back(Order(_ex1)); break; } else ! sum += ((p * j - (i - j)) * co[i - j] * c).expand(); } ! co.push_back(sum); } // Construct new series (of non-zero coefficients) -- Frieder Kalisch Institut für theoretische Physik kalisch@tphys.uni-heidelberg.de Philosophenweg 19 +49-6221-549-433 D-69120 Heidelberg
Hi! On Wed, Jan 19, 2005 at 06:31:41PM +0100, Chris Dams wrote:
Ad (2): I think that the expansion coefficients should be in expandend form (again in the sense of .expand()).
In ginsh: series((a+x)^(100), x==1, 3); ((1+a)^100)+(100*(1+a)^99)*(-1+x)+(4950*(1+a)^98)*(-1+x)^2+Order((-1+x)^3) I don't think anyone wants these coefficient to be .expand()'-ed... Best regards, Alexei.
Hello, Am Freitag 21 Januar 2005 13:10 schrieb Sheplyakov Alexei:
On Wed, Jan 19, 2005 at 06:31:41PM +0100, Chris Dams wrote:
Ad (2): I think that the expansion coefficients should be in expandend form (again in the sense of .expand()).
In ginsh:
series((a+x)^(100), x==1, 3); ((1+a)^100)+(100*(1+a)^99)*(-1+x)+(4950*(1+a)^98)*(-1+x)^2+Order((-1+x)^3)
I don't think anyone wants these coefficient to be .expand()'-ed...
Perhaps the coefficients should preserve the "expandedness" of the original expression. The patch below implements this. In each step of the composition of a series in pseries::mul_series and pseries::power_const the coefficients are multiplied term by term. In this way deeply nested expressions are avoided. Another way to implement the same feature would be to add an option to expand to not descend into children . Any comments about these two suggestions and the patch? Please note that I am not aware of any coding standards within ginac. Regards Frieder Kalisch *** ginac/pseries.cpp.save 2005-01-18 10:13:46.000000000 +0100 --- ginac/pseries.cpp 2005-02-03 13:41:45.000000000 +0100 *************** *** 693,698 **** --- 693,738 ---- } + /** Helper struct for multiply_term_by_term. + * map_function, that multiplies its argument by a constant. If the + * argument or the constant is a sum, they are multiplied term by term. */ + struct multiply_by : public map_function { + ex c; + + multiply_by(const ex &e) : c(e) {} ; + + ex operator()(const ex &e) + { + if (is_a<add>(c)) { + multiply_by mult_func(e); + return c.map(mult_func); + } + + return e*c; + } + + }; + + + /** Helper function, that uses struct multiply_by. + * Multiply two expressions. If either is a sum, they are multiplied + * term by term. */ + ex multiply_term_by_term(const ex &e1, const ex &e2) + { + if (is_a<add>(e1)) { + multiply_by mult_func(e2); + return e1.map(mult_func); + } + + if (is_a<add>(e2)) { + multiply_by mult_func(e1); + return e2.map(mult_func); + } + + return e1*e2; + } + + /** Multiply a pseries object with a numeric constant, producing a pseries * object that represents the product. * *************** *** 756,762 **** ex a_coeff = coeff(var, i); ex b_coeff = other.coeff(var, cdeg-i); if (!is_order_function(a_coeff) && !is_order_function(b_coeff)) ! co += a_coeff * b_coeff; } if (!co.is_zero()) new_seq.push_back(expair(co, numeric(cdeg))); --- 796,802 ---- ex a_coeff = coeff(var, i); ex b_coeff = other.coeff(var, cdeg-i); if (!is_order_function(a_coeff) && !is_order_function(b_coeff)) ! co += multiply_term_by_term(a_coeff, b_coeff); } if (!co.is_zero()) new_seq.push_back(expair(co, numeric(cdeg))); *************** *** 893,906 **** for (int i=1; i<deg; ++i) { ex sum = _ex0; for (int j=1; j<=i; ++j) { ! ex c = coeff(var, j + ldeg); if (is_order_function(c)) { co.push_back(Order(_ex1)); break; } else ! sum += (p * j - (i - j)) * co[i - j] * c; } ! co.push_back(sum / coeff(var, ldeg) / i); } // Construct new series (of non-zero coefficients) --- 933,946 ---- for (int i=1; i<deg; ++i) { ex sum = _ex0; for (int j=1; j<=i; ++j) { ! ex c = multiply_term_by_term(coeff(var, j + ldeg), _ex1/i/coeff(var, ldeg)); if (is_order_function(c)) { co.push_back(Order(_ex1)); break; } else ! sum += multiply_term_by_term((p * j - (i - j)) * co[i - j], c); } ! co.push_back(sum); } // Construct new series (of non-zero coefficients) -- Frieder Kalisch Institut für theoretische Physik kalisch@tphys.uni-heidelberg.de Philosophenweg 19 +49-6221-549-433 D-69120 Heidelberg
participants (3)
-
Chris Dams
-
Frieder Kalisch
-
varg＠thsun1.jinr.ru