I haven't studied all the notes prior to this, but it may be helpful to look at Macsyma/ Maxima. Series can be extended to several variables in different ways, e.g. series in x to order xn with coefs as series in y to order yn.... etc
or to total order that is degree in(x) + degree in (y) + .... <=N Also the index set could be other than the integers 0,1,2, ... including negative and fractional powers. Also, a key procedure in some algorithms is to find the lowest non-zero coefficient, a task which can be easy, hard, or non-computable. So you might have to think about this. Others have resolved it by heuristics. Apologies if these comments have already been taken into account or if they are irrelevant in this context. RJF On Wednesday, February 11, 2015 at 5:43:17 PM UTC-8, Ondřej Čertík wrote: > > Hi Mario, > > On Tue, Nov 26, 2013 at 5:38 AM, mario <[email protected] <javascript:>> > wrote: > > > > In PR 609 and PR 2630 power series are implemented with dictionaries, > that > > is in a sparse representation. > > The sparse representation is efficient when there are many variables. PR > > 2630 uses ``ring()``. > > > > For power series in one variable the dense representation is usually > faster. > > In `compose_add` in PR 2630 series in one variable are used, so using > the > > dense representation > > would be faster, but in a negligible way, since `compose_add` takes > roughly > > 5% of the time in `minpoly`. > > > > In PR 609 there is lpoly.py, which deals with power series, and > ltaylor.py. > > lpoly.py can now be implemented using ``ring()``; this is partially done > in > > ring_series.py in PR 2630, > > where there are only the functions used in the computed sum. It seems to > me > > that ring_series.py fits > > well in SymPy; if one wants to compute efficiently power series, > > one can use ``ring`` with ring_series.py functions (once all the series > > functions in lpoly.py are ported). > > > > To do a thorough treatment of power series, one should also implement > > the dense representation; consider the case of Laurent series; > > add functions to manipulate series at a higher level, possibly using a > > Series class with an Order function. > > > > Mateusz Paprocki might want to comment on this. > > > > PowerSeriesRing in Sage could be an example of higher level > representation. > > > > ltaylor.py uses lpoly.py to compute the taylor series of SymPy > functions; > > if the function is not smooth enough, the computation is passed to > > series.series. > > So if the function is smooth enough, the taylor series is computed fast > > (more or less the speed of `taylor` in Sage); otherwise it can be > > slightly slower than series.series. > > I do not know if it is a good approach; anyway, since it uses power > series, > > for the moment I will leave it aside. > > I finally got to this. I checked out your branch: > > https://github.com/sympy/sympy/pull/609 > > Well, actually the Mateusz' improved branch, and I updated it further > to run with the latest gmpy2: > https://github.com/certik/sympy/tree/lpoly2-improved > > And tested an example on my computer with gmpy: > > In [1]: %time s = series(sin(x)*cos(x), x, 0, 1000) > CPU times: user 4.52 s, sys: 14 ms, total: 4.54 s > Wall time: 4.5 s > > In [2]: %time s = series(log(1+x)*exp(x), x, 0, 500) > CPU times: user 4.15 s, sys: 14 ms, total: 4.16 s > Wall time: 4.13 s > > > Then I ran the same benchmark using the latest sympy master, using > your code that got merged at https://github.com/sympy/sympy/pull/2630, > I put the scripts here: > https://github.com/certik/sympy/compare/formal_power_series, it's just > the a.py and b.py. Your merged code doesn't seem to support sin/cos > benchmark [1] above, so I just ran the benchmark [2]: > > > https://github.com/certik/sympy/blob/59492520b443ea5f0ef31fc018e9bc700b93b818/b.py > > > $ python b.py > import > done > initialize series > done > start > stop > 1.85430693626 > > Seems to be a little bit faster. Finally, I then wrote my own > implementation: > > > https://github.com/certik/sympy/blob/59492520b443ea5f0ef31fc018e9bc700b93b818/a.py > > > $ python a.py > import > done > initialize series > done > start > 1.17188405991 > stop > 1.17299985886 > > > It's using the same algorithm, the speed improvement is just caused by > sorting the dictionary first before multiplying things out in rs_mul: > > items1 = list(p1._data.items()) > items1.sort(key=lambda e: e[0]) > items2 = list(p2._data.items()) > items2.sort(key=lambda e: e[0]) > > > You only sort the items2 in sympy, while I noticed a speedup if I sort > items1 as well. But this is minor. > > > All this is using sparse representation, as a dictionary. I also > implemented dense representation, e.g. you can try the dense > representation on the example [2] by applying the following patch: > > diff --git a/a.py b/a.py > index 5fabdd9..f7caf4a 100644 > --- a/a.py > +++ b/a.py > @@ -179,20 +179,20 @@ def div(a, b): > for i in range(1, n): > t = t/i > data.append(t) > -exp = FormalPowerSeriesSparse.from_dense(data) > +exp = FormalPowerSeries(data) > # log(1+x) > data = [QQ.dtype(0)] > t = QQ.dtype(1) > for i in range(1, n): > data.append(t/i) > t = -t > -log = FormalPowerSeriesSparse.from_dense(data) > +log = FormalPowerSeries(data) > #x = FormalPowerSeries([0, 1, 0, 0, 0, 0, 0, 0, 0]) > #onemx = FormalPowerSeries([1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) > print "done" > print "start" > t1 = clock() > -s = mul_sparse(log, exp, n) > +s = mul(log, exp) > #s = mul(sin, cos) > t2 = clock() > print "stop" > > > Then you get: > > > $ python a.py > import > done > initialize series > done > start > 1.11554312706 > stop > 1.11559510231 > > So it's a little bit faster, but the difference is small. > > Conclusion --- a.py now contains simple implementation of the taylor > series that is as fast as Mario's fast code. > > > Anyway, I think this is the right approach, we should support more > functions (sin, cos, ...) in sympy.polys.ring_series and then create a > class like the FormalPowerSeriesSparse() in a.py that would hold the > series and would know how to print it, and also had methods to > manipulate it (or those can be external functions like in a.py). > Finally, this should then be hooked into our series expansion in > sympy. > > Now, this works great for a Taylor expansion. We then need to extend > this to support Laurent series and eventually Puiseux series (that's > what SymPy's series is doing). An example of a Puiseux series is: > > In [1]: %time print (1/cos(x/log(x))).series(x, 0, 10) > 1 + x**2/(2*log(x)**2) + 5*x**4/(24*log(x)**4) + > 61*x**6/(720*log(x)**6) + 277*x**8/(8064*log(x)**8) + O(x**10) > CPU times: user 3 s, sys: 14 ms, total: 3.01 s > Wall time: 3.03 s > > As you can see, currently it takes 3s, and it should be immediate in > my opinion. You can read the references that I quoted in the docstring > of FormalPowerSeries() in a.py. Essentially, I think these are all > series of the form a*X^n, where n = 0, 1, 2, ... for Taylor series, n > = n0, n0+1, n0+2, ... for negative n0 for Laurent series, and X is a > symbol. For Puiseux series, the X is then something more complicated, > either a fractional power like x^(1/2), or something like x/log(x) in > the last example. But the contents of X is stored separately, and the > rest is just like polynomials. I am currently not sure if things like > multiplication of two Puiseux series needs to be changed somehow. > > Sparse/Dense: I would just go with Sparse, it seems the speed is very > good both for really sparse series, as well as dense series. > > Overall, I think this is the way to go. Let me know what you think. > > Ondrej > -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/c063f642-635f-4fc8-aee5-a119e6312382%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
