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
<