Implement infinite sums and products for lazy series

[tscrim]

We allow lazy series to compute infinite sums \sum_{i \in I} p_i and products \prod_{i \in I} (1 + p_i) over arbitrary (enumerated) index sets I subject to a somewhat mild technical constraint that the order of each input p_i is weakly increasing wrt the iteration order and the preimage of each order is finite.

📝 Checklist

• The title is concise, informative, and self-explanatory.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation accordingly.

⌛ Dependencies

• Improve approximate order on getitem calls #35265

[tscrim]

@mantepse Here is a sketch for infinite products using lazy series based upon our conversion about this from many months ago. There are likely many ways it could be improved, but I wanted you to take a look at it before I started to polish it and copy to infinite sums.

In principle, I could get rid of the _ring dependency (and being a LazySeries type object) by having it only work with streams, but this was very cumbersome to implement and means we had to separate out the non-Cauchy products. In many ways, this is just a Stream_function wrapping around a very fancy iterator with allowing random access. It seems like I should have an InfiniteProductList or similar such class within lazy_series.py that is then passed to a Stream_function, but this way cuts out that as a middleman.

Anyways, I am off to bed. If you get a chance today, let me know what you think.

[mantepse]

I will look properly tomorrow - maybe you could meanwhile merge develop?

One thing: essentially you only use self._ring.one(). Maybe it would be better (and more in the spirit of stream.py) to just pass this, i.e., one?

Also, it would be good to include a check that it interacts nicely with recursive definitions. Spoiler: looks good.

sage: f = L.undefined(valuation=0)
sage: f.define(L.prod(lambda n: t^n*f, PositiveIntegers()))
sage: f
1 + t + 2*t^2 + 5*t^3 + 12*t^4 + 30*t^5 + 77*t^6 + O(t^7)
sage: f[:10]
[1, 1, 2, 5, 12, 30, 77, 201, 532, 1427]
sage: oeis(f[:10])
0: A145267: G.f. satisfies A(x) = Product_{k>0} (1+x^k*A(x)).

In the docstring of prod, index=True should be index_set=True.

Did you check whether padics has similar functionality?

One thing that I do not really like is the dependency of the semantics on the input:

sage: L.prod([t^k for k in range(5)])
t^10

Wouldn't it be better to separate the two functionalities? I don't know whether there is an established name for the infinite product of this specific form, though.

[mantepse]

Possibly euler_product?

Besides, here is possibly a bug:

sage: D = LazyDirichletSeriesRing(QQ, "s")
sage: D.prod(lambda p: (1+D(1, valuation=p)).inverse()-1, index_set=Primes())

does not return.

[tscrim]

I thought a lot about how to name it, and I couldn't come up with something natural beyond prod() (since it should also work for Dirichlet series). Ideally, it would be able to take in an infinite iterator and know to pass it off to the infinite product or just evaluate it out explicitly. However, that means we have to know whether an iterator is infinite or not, which is an impossible problem. Consequently I think of this prod as extending the current functionality, but if you feel it is better to separate the two, euler_product could work.

Indeed, that seems like a bug. It does "return" but it seems to have trouble getting coefficients. It probably has to do with not taking in the minimal valuation. This likely should be passed to the stream (which would be the default for _cur_order).

[tscrim]

Also, I did not check if padics has similar functionality. I agree with you about passing one instead of the ring. I thought I would use the ring more.

[mantepse]

Another (different) approach: wouldn't it be equally good to require that the factors passed are of the form $1+ f(x)$ with the valuation of $f$ increasing? In the Dirichlet series example above, I have to manually subtract $1$ to get it into the required form.

I guess that your approach is slightly faster, but then we could have two methods, prod which just does the product and euler_prod which first adds one.

[tscrim]

Another (different) approach: wouldn't it be equally good to require that the factors passed are of the form $1+f(x)$ with the valuation of $f$ increasing? In the Dirichlet series example above, I have to manually subtract 1 to get it into the required form.

I thought about this, but I quickly came to the conclusion it is not a good idea. We would have to either manually implement something to check the next nonzero coefficient (and essentially duplicate our valuation() code) or subtract one and do a computation on something we are not actually going to use (which could also run into trouble with exact versus looks-like-1, say, with coefficients being lazy series). Indeed, it is not the "natural" input (and why you need to subtract 1), but the implementation aspects suggest something better.

I guess that your approach is slightly faster, but then we could have two methods, prod which just does the product and euler_prod which first adds one.

I also believe it to be faster. It is also better in terms of memory as the caches don't store all of those intermediate 0's. Moreover, I do not think the behaviors of the two methods should be different. In the proposal, the user is explicitly passing a different signature, so it is far less of a "gotcha" IMO.

TL;DR We can make it more natural by having $1 + f(x)$ as input, but there are tricky details with such an implementation.

[tscrim]

Perhaps another option is we assume the user is giving good input (which we very clearly must warn them about). It would have an advantage of avoiding a trivial +1 - 1 computation in your example. However, we would still have to manually (or through some transient object) find the valuation of $g(x) = 1 - (1 + f(x))$.

[tscrim]

Perhaps it actually isn't so bad; I forgot that I did a custom check for the valuation on the components since I didn't need to compute their actual order.

[mantepse]

I could not find anything in padics.

I would still prefer to have a separate method for products of the form $\prod_i 1+f_i(x)$. I do not see the benefit of having the same name for these. It seems fitting to me to call this method euler_product, even though this has a much more specialised meaning in the context of Dirichlet series.

I agree that accepting the factors as input (i.e., including the 1) doesn't look very attractive, either.

[tscrim]

I strongly think that all infinite products should be under prod; it is by far the most natural place one would look. If we just explicitly assume the user has passed good input (which I already do by assuming that every $f_i(x) \neq 0$) of the form $\prod_i p_i(x)$ with $p_i(x) = 1 + f_i(x)$ with $f_i(x)$ having weakly increasing orders, then I think I can make it work with very minimal changes overall (at the cost of cache inflation for dense cases).

[mantepse]

If it must be, and you make prod accept factors, that's OK. I still think that having a second method would be a good idea, but I really don't want to argue about this. I am only slightly worried that at some point we could have symbolic infinite products and lazy infinite products mean (slightly) different things.

Currently sage doesn't seem to have any infinite product at all.

[tscrim]

I changed it to have the iterator output $p_i(x) = 1 + f_i(x)$ and fixed the Dirichlet series bug (which was not limited to those series actually). However, in order to have these factors, I had to remove a safeguard against products involving smaller order $f_i(x)$ as I did not think it was useful to force the computation of all of those coefficients. We are already having to assume a fair amount of "good" user input for this.

To get rid of the infinite loop trap, I can change it so the order is only increased one step at each time. Although for something more "sparse" in regards to the orders of $f_i(x)$, this will mean multiplying more factors sooner, which can get very computationally expensive I think.

Anyways, that's all for me tonight; I need some sleep.

[tscrim]

If we are mostly okay with the general format of the infinite product, I will also add the infinite sums now too.

[tscrim]

I added the doctests and fixed a few small things (it handles trivial zeros in the iterator). This is now ready for review.

[tscrim]

Closes #34404.

[mantepse]

I discovered a problem, which I don't understand yet:

L.<z> = LazyPowerSeriesRing(QQ)
T = L.undefined(1)
D = L.undefined(0)
T.define(z*exp(T)*D)
H = L(lambda n: sum(T(z^k)/k for k in range(2, n+1))[n], valuation=2)
D.define(exp(H))

works, but

L.<z> = LazyPowerSeriesRing(QQ)
T2 = L.undefined(1)
D2 = L.undefined(0)
T2.define(z*exp(T2)*D2)
H2 = L.sum(lambda k: T2(z^k)/k, 2)
D2.define(exp(H2))

raises a StopIteration. It might be a user error, of course, but I don't see it yet.

[tscrim]

That is a good example. It shows that my implementation was computing coefficients when it shouldn't have been. I have changed things around so that it computes fewer coefficients unless it absolutely needs to. The drawback is that we need to be a bit more careful with having the approximate order weakly increasing when using undefined series in infinite sums/products.

[github-actions]

Documentation preview for this PR (built with commit 2afd7f3; changes) is ready! 🎉

[mantepse]

Looks good!

[tscrim]

Thank you.