[GiNaC-devel] series expansion

26 Apr 2006

      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
<

[GiNaC-devel] series expansion

Jens Vollinga