Hi, in order to avoid doubled effort I'd like to make an announcement: I am currently working on a better series expansion mechanism using recurrence relations. Attached is a patch that shows the principal idea (cf. TAOCP or papers from A.C.Norman). This code is very veeeery poorly written and it won't go into CVS. It's even very buggy ... Don't mess with the argument of sqrt! So, it is just a show case (PIV, 3GHz): doing series( sqrt( (1-x)/(1+x) ), x, N) takes N=50: old: 54s new: 0.2s N=100: old: 32m3s (and 360MB of RAM ...!) new: 1.4s In the end, the recurrence relations should be determined automatically by GiNaC for any (nice) function or operation together with an on-demand calculation of arbitrary expansion coefficients (cf. Norman 1975). BTW, has anybody thought about ways to do Puiseux series? If one doesn't care about the strange and obscure functions but concentrates on a practical solution for "simple" functions like exp(sqrt(x)) for example, it seems to me that the only thing needed would be an additional number kept with every pseries that represents the denominator of non-integer exponents (like 2 in the case of exp(sqrt(x))). Of course, this reduces the performance for the small(?) benefit of rarely(?) used Puiseux expansions. Regards, Jens Index: pseries.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/pseries.cpp,v retrieving revision 1.88 diff -r1.88 pseries.cpp 1039a1040,1092
ex res; // series expansion by recurrence relations -- work in progress if (r.rhs().is_zero() && (exponent == _ex1_2)) { const ex& xpt = r.lhs(); std::vector<ex> v(order+1); ex basisexp = basis.series(r, order, options); for (int i=0; i<order; ++i) { v[i] = basisexp.coeff(xpt, i); } std::vector<ex> w(order+1); w[0] = basis.subs(r);
// prep for recurrence relation for (int i=1; i<=order; ++i) { w[i] = symbol(); }
symbol x("x"); ex V; ex W; for (int i=0; i<=order; ++i) { V += v[i] * pow(x, i); W += w[i] * pow(x, i); } ex Vp; ex Wp; for (int i=1; i<=order; ++i) { Vp += i * v[i] * pow(x, i-1); Wp += i * w[i] * pow(x, i-1); }
// W = sqrt(V) ex dgl = 2 * W * Wp - Vp; dgl = dgl.expand();
epvector ser; ser.push_back(expair(w[0], 0)); for (int i=0; i<order-1; ++i) { ex co = dgl.coeff(x, i); ex co2 = -co / co.coeff(w[i+1]); ex co3 = co2.subs(w[i+1]==0); co3 = co3.normal();
dgl= dgl.subs(w[i+1] == co3); w[i+1] = co3; ser.push_back(expair(co3, i+1)); } ser.push_back(expair(Order(xpt), order)); res = pseries(r, ser);
return res; }
1061d1113 <
Dear Jens, On Wed, 26 Apr 2006, Jens Vollinga wrote:
in order to avoid doubled effort I'd like to make an announcement: I am currently working on a better series expansion mechanism using recurrence relations.
I think this is a great idea. I have since long thought that the current series expansion is bad.
BTW, has anybody thought about ways to do Puiseux series? If one doesn't care about the strange and obscure functions but concentrates on a practical solution for "simple" functions like exp(sqrt(x)) for example, it seems to me that the only thing needed would be an additional number kept with every pseries that represents the denominator of non-integer exponents (like 2 in the case of exp(sqrt(x))). Of course, this reduces the performance for the small(?) benefit of rarely(?) used Puiseux expansions.
Maybe the current representation of a series could be turned into something that uses an exvector instead of an expairseq. You would store the order at which the series starts and the stepsize by which the exponent increases for every term. These could both be rational numbers. This may waste some space in case that, say, half the coefficients are zero, but now space is often wasted by storing the exponents in cases where they are subsequent natural numbers. Best, Chris
participants (2)
-
Chris Dams
-
Jens Vollinga