Lazy Laurent series

[klee]

Introduce lazy Laurent series to Sage.

A lazy Laurent series computes coefficients only when demanded or needed. In a sense, lazy Laurent series are Laurent series of infinite precision.

A generating function example from the code:

Generating functions are laurent series over the integer ring::

    sage: from sage.rings.lazy_laurent_series_ring import LazyLaurentSeriesRing
    sage: L = LazyLaurentSeriesRing(ZZ, 'z')

This defines the generating function of Fibonacci sequence::

    sage: def coeff(s, i):
    ....:     if i in [0, 1]:
    ....:         return 1
    ....:     else:
    ....:         return s.coefficient(i - 1) + s.coefficient(i - 2)
    ....:
    sage: f = L.series(coeff, valuation=0); f
    1 + z + 2*z^2 + 3*z^3 + 5*z^4 + 8*z^5 + 13*z^6 + ...

The 100th element of Fibonacci sequence can be obtained from the generating
function::

    sage: f.coefficient(100)
    573147844013817084101 

Lazy Laurent series is of course useful for other things. This will be used to implement infinite precision Laurent series expansion of algebraic functions in function fields, as a sequel to #27418.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​5d24a34

Introduce lazy laurent series

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​ace4282

pyflakes fixes

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​389798e

Add series method

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​43ba785

Small fixes on docstrings

[embray]

Ticket retargeted after milestone closed (if you don't believe this ticket is appropriate for the Sage 8.8 release please retarget manually)

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​e7ff999

Add lazy laurent series

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​15e1b1c

Add lazy laurent series

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f70fca8

Add __bool__

[tscrim]

Some comments:

Is there a reason why you do not use the LazyPowerSeriesRing for the main element (with a non-zero constant term until the whole series is zero) and then just store the valuation? This is how things are done for the Laurent polynomials. That way there is less duplication of code and improvements to one improve both.

Why do you override __eq__ and not use the one provided by Element and the coercion framework? It is more flexible and allows you do implement comparisons (if you want) all in one function of _richcmp_.

Do not have bare except: statements (see coefficient).

Does the elements pickle? I am a little worried that the addition of two power series will not pickling because of the little helper add function. You should add a loads(dumps(foo)) test if it does work.

Typo laurent -> Laurent (it is proper noun).

if len(s) == 0: -> if s:

[klee]

Replying to tscrim:

Is there a reason why you do not use the LazyPowerSeriesRing for the main element (with a non-zero constant term until the whole series is zero) and then just store the valuation? This is how things are done for the Laurent polynomials. That way there is less duplication of code and improvements to one improve both.

I looked into the LazyPowerSeriesRing many years ago:

​https://groups.google.com/forum/#!searchin/sage-support/lazy$20power$20series%7Csort:date/sage-support/JZwDWKq3DtM/tyfxxdxPyXYJ

I think I didn't like it in several ways :-) Still I don't understand this behavior:

sage: L.<t>=LazyPowerSeriesRing(QQ)
sage: s=L([0,0,1,2])
sage: s.coefficient(0)
0
sage: s.coefficient(1)
0
sage: s.coefficient(2)
1
sage: s.coefficient(3)
2
sage: s
t^2 + 2*t^3 + O(x^4)
sage: s.get_order()
1
sage: s.get_aorder()
1

The reason was that it was easy to implement lazy laurent series with an interface and behaviors that I prefer, while I could not understand well enough LazyPowerSeriesRing...

I would reconsider basing lazy laurent series on LazyPowerSeriesRing, if you help me understand the strange behaviour above.

[tscrim]

So here is what refine_aorder does. It first checks if self.order has been set, if it has, do not do anything. Then it goes through all of the n computed entries if self.aorder is 0 and we have already computed at least one entry. Then find the first nonzero entry. If we have found one, set the order.

The problem is the bold part. After the first iteration when n = 1, the self.aorder gets set to 1, then on the second run through, that block does not get run, which then means self.aorder < n and we set the order to be 1.

sage: s.coefficient(1)
0
sage: s.aorder
1
sage: s.order
Unknown series order
sage: s.coefficient(2)
1
sage: s.aorder
1
sage: s.order
1

IMO, this is a bug.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​2f2c59c	Capitalize laurent as it is a proper noun
​146205a	No bare except clause
​fb4d65d	Use _richcmp_ instead of primitive __eq__

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​41f3a42

Lazy Laurent series is not picklable

[klee]

Why do you override __eq__ and not use the one provided by Element and the coercion framework? It is more flexible and allows you do implement comparisons (if you want) all in one function of _richcmp_.

Ok. Fixed.

Do not have bare except: statements (see coefficient).

Fixed.

Does the elements pickle? I am a little worried that the addition of two power series will not pickling because of the little helper add function. You should add a loads(dumps(foo)) test if it does work.

No. Your worry is real; lazy Laurent series is not picklable in general. This seems an inherent nature of my implementation...

I added a test to show this.

Typo laurent -> Laurent (it is proper noun).

Right. Fixed.

if len(s) == 0: -> if not s:

Done.

---

New commits:

​41f3a42

Lazy Laurent series is not picklable

[klee]

Replying to tscrim:

Some comments:

Is there a reason why you do not use the LazyPowerSeriesRing for the main element (with a non-zero constant term until the whole series is zero) and then just store the valuation? This is how things are done for the Laurent polynomials. That way there is less duplication of code and improvements to one improve both.

I thought about this. Now I think:

Lazy Laurent series is not just about Laurent series but also about how coefficients are computed lazily. A LazyPowerSeriesRing element gets its coefficients from a stream object, which essentially yields coefficients as required. On the other hand, lazy Laurent series s gets its coefficients from a python function that outputs n-th coefficient for input s and n. This allows coefficients to be computed recursively. For example, it is very easy to define the Fibonacci series.

So it is impossible to base my lazy Laurent series code on LazyPowerSeriesRing without abandoning the essential feature of the implementation.

Frankly, I would rather provide new lazy power series as a subclass of my lazy Laurent series. But this is of course highly biased opinion :-)

You may be reluctant to accept my implementation of lazy Laurent series into sage.rings as it then gets kind of standard status among different possible implementations of lazy Laurent series in Sage. Then I am willing to relocate it into sage.rings.function_field.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b09cd84

Fix pickling doctest

[tscrim]

Replying to klee:

Replying to tscrim:

Some comments:

Is there a reason why you do not use the LazyPowerSeriesRing for the main element (with a non-zero constant term until the whole series is zero) and then just store the valuation? This is how things are done for the Laurent polynomials. That way there is less duplication of code and improvements to one improve both.

I thought about this. Now I think:

Lazy Laurent series is not just about Laurent series but also about how coefficients are computed lazily. A LazyPowerSeriesRing element gets its coefficients from a stream object, which essentially yields coefficients as required. On the other hand, lazy Laurent series s gets its coefficients from a python function that outputs n-th coefficient for input s and n. This allows coefficients to be computed recursively. For example, it is very easy to define the Fibonacci series.

Well, I am fairly certain you could do this with an RecursivelyEnumeratedSet and a bit of know-how, but I agree that it is more complicated and less intuitive. I also agree that having general functions as input would probably be a good thing for LazyPowerSeriesRing to have, but I believe that was designed for more specific use for the combinat/species code.

One benefit I do see is that you do not need to explicitly know the valuation at creation-time. With my suggested implementation, it does need to be computed at element creation.

Frankly, I would rather provide new lazy power series as a subclass of my lazy Laurent series. But this is of course highly biased opinion :-)

I don't have an opinion on this matter, but there are likely some considerations required to replace or extend the current implementation.

You may be reluctant to accept my implementation of lazy Laurent series into sage.rings as it then gets kind of standard status among different possible implementations of lazy Laurent series in Sage. Then I am willing to relocate it into sage.rings.function_field.

This I don't really care about. I just am a little unhappy with the large difference between the semantics (and syntax) between the two implementations. I guess that is a bit unavoidable here because this definitely is serving a purpose.

However, I do think we should at least try to address the pickling issues. I believe this means you cannot use this in parallel implementations (IIRC this uses pickling to communicate between the processes). It also means you cannot (easily) store the data you compute. I think it is sufficient to separate out the add and similar operations into functions, but maybe they need to be small little helper class, such as

class LaurentSeriesOperator(object):
    def __init__(self, lps, op):
        self.lps = lps
        self.op = op
    def __call__(self, s, n):
        return self.op(self.lps[n], s[n])
    def __reduce__(self):
        return (type(self), (self.lps, self.op), {})
    def __eq__(self, other):
        return (isinstance(other, LaurentSeriesOperator)
                and self.lps == other.lps and self.op == other.op)

where op is, e.g., operator.add. This way you might be able to do something with comparisons in some semi-reasonable capacity too.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​89992b8

Reimplement arithmetic operations for series to be picklable

[klee]

Replying to tscrim:

However, I do think we should at least try to address the pickling issues. I believe this means you cannot use this in parallel implementations (IIRC this uses pickling to communicate between the processes). It also means you cannot (easily) store the data you compute. I think it is sufficient to separate out the add and similar operations into functions, but maybe they need to be small little helper class, such as

class LaurentSeriesOperator(object):
    def __init__(self, lps, op):
        self.lps = lps
        self.op = op
    def __call__(self, s, n):
        return self.op(self.lps[n], s[n])
    def __reduce__(self):
        return (type(self), (self.lps, self.op), {})
    def __eq__(self, other):
        return (isinstance(other, LaurentSeriesOperator)
                and self.lps == other.lps and self.op == other.op)

where op is, e.g., operator.add. This way you might be able to do something with comparisons in some semi-reasonable capacity too.

Using module-level classes to define operators is a good idea. Hinted by your template class, I could reimplement lazy Laurent series to be picklable. Great!

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​989a3ac

Comparison now works in semi-reasonable capacity

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​ce4db37	Now pass _test_pickling
​e12fd5b	Minor fixes in string formatting

[git]

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

​afeae09	Merge branch 'function-fields-infinity-precision' into function-fields
​bf81509	Merge branch 'lazy-power-series-trac27347' into function-fields
​e552770	fix
​6c5fc73	Merge branch 'function-fields-infinity-precision' into function-fields
​f616d20	Merge branch 'lazy-power-series-trac27347' into function-fields-infinity-precision
​8d77ecc	Merge branch 'function-fields-infinity-precision' into function-fields
​bbeebef	Merge branch 'function-fields' into curves-irreducible
​64c9c6b	Add delta-invariant
​6142750	Add ag codes
​5214a3a	Add list operator and allow list input

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​eac7047

Add list operator and allow list input

[klee]

Oops. Pushed a wrong branch!

[tscrim]

This is looking good.

Past experience tells me you should also implement a __ne__ at the base class. I also think you should include more ABCs to avoid the extra duplication with the __eq__ methods. Similarly, why not use the op version for +-*/ to avoid class duplication? You may also want a __hash__ and incorporate that into the hash of the lazy Laurent series. With how many operator classes you have, it may also be worthwhile pulling that out to a separate file for better separations-of-concerns between files.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f5b5792

Add hash and fix other issues from reviewer comments

[klee]

Replying to tscrim:

This is looking good.

Past experience tells me you should also implement a __ne__ at the base class.

Done.

I also think you should include more ABCs to avoid the extra duplication with the __eq__ methods.

Ok.

Similarly, why not use the op version for +-*/ to avoid class duplication?

I don't get it. They cannot be treated uniformly. Look at the __call__ method.

You may also want a __hash__ and incorporate that into the hash of the lazy Laurent series.

Done.

With how many operator classes you have, it may also be worthwhile pulling that out to a separate file for better separations-of-concerns between files.

Done.

[tscrim]

Replying to klee:

Replying to tscrim:

Similarly, why not use the op version for +-*/ to avoid class duplication?

I don't get it. They cannot be treated uniformly. Look at the __call__ method.

You can do what I suggested in comment:20 for these classes:

• LazyLaurentSeriesOperator_add (by passing operator.add)
• LazyLaurentSeriesOperator_sub (by passing operator.sub)

You're right that multiplication must be treated separately. Although you could also do this for a scalar multiplication version of this (which will be faster than coercion and multiplication). For 2 cases, there is a less compelling reason to factor the common code out like this.

There is not really a benefit for doing this for unitary operations (as there is only 1 that could work this way).

In some ways it feels like we are reinventing the wheel here because we are effectively building an evaluation tree, and we already have something like this within SR (or sympy/pynac/etc.). However, I don't see a way around this by using SR (and we might have to work a little bit to avoid doing computations unnecessarily inside of there). So I think we have to keep doing things this way. Just an observation.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3bf2101

Add scalar multiplication operator

[klee]

Replying to tscrim:

Although you could also do this for a scalar multiplication version of this (which will be faster than coercion and multiplication).

Added scalar multiplication.

In some ways it feels like we are reinventing the wheel here because we are effectively building an evaluation tree, and we already have something like this within SR (or sympy/pynac/etc.). However, I don't see a way around this by using SR (and we might have to work a little bit to avoid doing computations unnecessarily inside of there). So I think we have to keep doing things this way. Just an observation.

True. It is just inevitable that sometimes we make wheels for different needs :-)

[mantepse]

I find this quite wonderful! It would be great if we could replace LazyPowerSeries with this. One feature of LazyPowerSeries was that it allows recursive definitions "guessing" the correct initialisations, as for example in

            sage: L = LazyPowerSeriesRing(QQ)
            sage: one = L(1)
            sage: monom = L.gen()
            sage: s = L()
            sage: s._name = 's'
            sage: s.define(one+monom*s*s)
            sage: [s.coefficient(i) for i in range(6)]
            [1, 1, 2, 5, 14, 42]

Although this is based on the original code by Ralf Hemmecke and myself, which in turn is based on code of Nicolas Thiery (I think), I am not sure anymore whether it is so important. Of course, if we could have something similar, even better.

[mantepse]

One question with respect to equality. It seems to me that your implementation of __eq__ is quite brittle. It is clear that one cannot have proper equality, so it should be clearly stated what it can do and what it can't.

Apart from that, in the case of truncation, one could make it exact easily.

[tscrim]

Some other comments on the code:

For __bool__, I think erroring out immediately when self._constant is None is not the way to go. I think you should check the computed coefficients in that case to see if any of them are non-zero. If they are all 0, then you should error out. Also, a fast first check when self._constant is not None is to see if that constant is non-zero 0.

In _repr_, I do not see how self.valuation() would normally raise a ValueError. That method does not raise one, nor does self.coefficient() (the only other method it calls).

Instead of self.parent()(1), which calls the coercion framework, you should use self.parent().one().

Are you sure you want _div_ to be two separate operations and not one combined one? I feel like the latter would be faster (at least since it is implemented at the Python level).

Actually, you may want to make the operator file into a Cython file (since is forms the key computational part and it is in a separate file).

I think the truncate should return an honest Laurent polynomial. Furthermore, I don't think there are coercions (or at least conversions) from the corresponding Laurent polynomial ring. This feels like a very natural to have. I don't want to go too far with ticket expansion, but this feels like something people will try soon after having this.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​4bf0f2f	Add coercion from Laurent polynomials
​86ccdc4	Add polynomial method and coercions from Laurent polynomial ring

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​e028ed6

A little change in truncate method

[klee]

Replying to tscrim:

For __bool__, I think erroring out immediately when self._constant is None is not the way to go. I think you should check the computed coefficients in that case to see if any of them are non-zero. If they are all 0, then you should error out.

Good idea. Done.

Also, a fast first check when self._constant is not None is to see if that constant is non-zero 0.

This check is already there.

In _repr_, I do not see how self.valuation() would normally raise a ValueError. That method does not raise one, nor does self.coefficient() (the only other method it calls).

Right. Fixed.

Instead of self.parent()(1), which calls the coercion framework, you should use self.parent().one().

Ok. Fixed.

Are you sure you want _div_ to be two separate operations and not one combined one? I feel like the latter would be faster (at least since it is implemented at the Python level).

I don't feel so. I don't see a way to implement series division such that 1 + 1 < 2 with a significant margin.

I think the truncate should return an honest Laurent polynomial.

I don't agree. One might still want a truncated series with laziness. For your cases, I added polynomial method to convert to an honest polynomial.

Furthermore, I don't think there are coercions (or at least conversions) from the corresponding Laurent polynomial ring. This feels like a very natural to have. I don't want to go too far with ticket expansion, but this feels like something people will try soon after having this.

Right. I added coercions from (Laurent) polynomial rings.

Actually, you may want to make the operator file into a Cython file (since is forms the key computational part and it is in a separate file).

Hmm. I doubt if there would be significant speed boost, compensating increased build time (that I hate). Moreover I am not versed with Cython. So I don't want to do that, at least for this ticket. But if you do it, I would not object.

[klee]

Replying to mantepse:

One question with respect to equality. It seems to me that your implementation of __eq__ is quite brittle. It is clear that one cannot have proper equality, so it should be clearly stated what it can do and what it can't.

The implementation gives an answer if possible without computing coefficients indefinitely. Otherwise raise an exception.

I added a short explanation to the method.

Apart from that, in the case of truncation, one could make it exact easily.

Truncation still gives a lazy series. I added polynomial method to convert to an exact polynomial.

[klee]

Replying to mantepse:

I find this quite wonderful! It would be great if we could replace LazyPowerSeries with this. One feature of LazyPowerSeries was that it allows recursive definitions "guessing" the correct initialisations, as for example in

            sage: L = LazyPowerSeriesRing(QQ)
            sage: one = L(1)
            sage: monom = L.gen()
            sage: s = L()
            sage: s._name = 's'
            sage: s.define(one+monom*s*s)
            sage: [s.coefficient(i) for i in range(6)]
            [1, 1, 2, 5, 14, 42]

Although this is based on the original code by Ralf Hemmecke and myself, which in turn is based on code of Nicolas Thiery (I think), I am not sure anymore whether it is so important. Of course, if we could have something similar, even better.

Recursive definition is possible as a result of being lazy and being power series. Now you can do the same in the following way:

sage: from sage.rings.lazy_laurent_series_ring import LazyLaurentSeriesRing
sage: L = LazyLaurentSeriesRing(ZZ, 'z')
sage: 
sage: z = L.gen()
sage: def f(s,n):
....:     return (1 + z*s^2).coefficient(n)
....: 
sage: L.series(f, valuation=0)
1 + z + 2*z^2 + 5*z^3 + 14*z^4 + 42*z^5 + 132*z^6 + ...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​8f8adf2

Add __getitem__ special method to series

[mantepse]

Apart from that, in the case of truncation, one could make it exact easily.

Truncation still gives a lazy series. I added polynomial method to convert to an exact polynomial.

excellent point! I missed that!

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b3a6792

Add division operator

[klee]

Are you sure you want _div_ to be two separate operations and not one combined one? I feel like the latter would be faster (at least since it is implemented at the Python level).

I don't feel so. I don't see a way to implement series division such that 1 + 1 < 2 with a significant margin.

You are right. I was wrong!

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​575022c

Typos

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​c55925f

Added a recursive definition example

[mantepse]

Are composition and reversion planned in a later ticket?

I am really looking forward to this! I hope very much we can rebase the species code on it. Should LazyPowerSeriesRing be a subclass?

[klee]

Replying to mantepse:

Are composition and reversion planned in a later ticket?

Not by me. Anyone is invited to add them in a later ticket.

I hope very much we can rebase the species code on it. Should LazyPowerSeriesRing be a subclass?

Yes, in a separate file alongside with lazy Laurent series. It would be great. I don't have a plan to do it in any short time, simply because I don't need it now.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​410ca38

Add prec and approximate_series method

[tscrim]

Last little thing:

R(0) -> R.zero() (as it is cached and R is always a LazyLaurentSeriesRing).

Once done you can set a positive review on my behalf.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​989f009

Fix R(0)

[klee]

Thank you, Travis!

[slelievre]

Nice work! Follow-up ticket at #27694 for allowing syntactic sugar

sage: L.<x> = LazyLaurentSeriesRing(ZZ)