a lazy solver for functional equations

[mantepse]

We enable sage to solve systems of functional equations of lazy power series, generalizing the existing possibility of defining series recursively. For example:

sage: L.<z> = LazyPowerSeriesRing(QQ)
sage: A = L.undefined()
sage: B = L.undefined()
sage: FA = A^2 + B^2 - 2 - z*B
sage: FB = B^3 + 2*A^3 - 3 - z*(A + B)
sage: L.define_implicitly([(A, [1]), (B, [1])], [FA, FB])
sage: A
1 + 1/6*z - 7/72*z^2 - 125/1296*z^3 - 3289/31104*z^4 - 23591/186624*z^5 - 1048607/6718464*z^6 + O(z^7)
sage: B
1 + 1/3*z + 7/36*z^2 + 47/324*z^3 + 1903/15552*z^4 + 2959/23328*z^5 + 500663/3359232*z^6 + O(z^7)

Dependencies: #37451, #37687, #37736, #37766, #37761, #37916

[tscrim]

I don't understand this change. Sorry if we've already discussed it and I don't remember.

[mantepse]

❤️ I am grateful that you keep working on this ❤️!

No, I don't think we discussed this specifically. The reason for the hack (which I hate) is as follows:

• recall that g is a tuple of (at least 2) elements of lazy series.
• if c is an element of the base ring B, c(g) cannot work correctly (and thus just return c) - e.g., integers are not callable, polynomials would give nonsense.
• if c is a proper element of the _laurent_poly_ring $R$, we are in the easy case
• otherwise, c should be an element of the modified _laurent_poly_ring. E.g., If $R=\mathbb Q[x,y]$, we have $c\in \mathbb Q[\vec a][x,y]$ (where $\vec a$ is FESDUMMY). In this case, to make the substitution work, I make sure that the elements of g are lazy series over the modified base ring, i.e. $\mathbb Q[\vec a]$.

I just saw that the error changed after #37916. We now have the following:

sage: U.<a> = InfinitePolynomialRing(QQ)
sage: R.<x,y> = U[]
sage: L.<z> = LazyPowerSeriesRing(QQ)
sage: c = x
sage: c(z, 0)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[7], line 1
----> 1 c(z, Integer(0))

File ~/sage/src/sage/rings/polynomial/multi_polynomial_element.py:182, in MPolynomial_element.__call__(self, *x, **kwds)
    180 y = K(0)
    181 for m, c in self.element().dict().items():
--> 182     y += c * prod(v ** e for v, e in zip(x, m) if e)
    183 return y
...
RuntimeError: There is a bug in the coercion code in Sage.
Both x (=1) and y (=z) are supposed to have identical parents but they don't.
In fact, x has parent 'Integer Ring'
whereas y has parent 'Infinite polynomial ring in a over Lazy Taylor Series Ring in z over Rational Field'
Original elements 1 (parent Infinite polynomial ring in a over Rational Field) and z (parent Lazy Taylor Series Ring in z over Rational Field) and maps
<class 'sage.structure.coerce_maps.DefaultConvertMap_unique'> (map internal to coercion system -- copy before use)
Coercion map:
  From: Infinite polynomial ring in a over Rational Field
  To:   Infinite polynomial ring in a over Lazy Taylor Series Ring in z over Rational Field
<class 'sage.categories.morphism.SetMorphism'> (map internal to coercion system -- copy before use)
Generic morphism:
  From: Lazy Taylor Series Ring in z over Rational Field
  To:   Infinite polynomial ring in a over Lazy Taylor Series Ring in z over Rational Field

[tscrim]

This code is so hacky and liable to pick up far more than it actually needs and likely could be quite slow because it needs to check potentially a lot of objects. I think we should consider caching the result. I would still return a list for consistency, but with a warning to the user not to mutate it (since it isn't part of the public facing stuff, only a more serious developer will get into this). Although perhaps that is slightly dangerous even within our code...? Perhaps these that need to do some computations be cached and every such method returns a tuple?

[mantepse]

This code is only executed once for each operation when creating the power series. For example,

sage: L.<z> = LazyPowerSeriesRing(SR)
sage: G = L.undefined(0)
sage: L.define_implicitly([(G, [ln(2)])], [diff(G) - exp(-G(-z))])
sage: %prun G[:20]
         532675 function calls (516627 primitive calls) in 0.815 seconds

   Ordered by: internal time
   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
...
        2    0.000    0.000    0.000    0.000 stream.py:1027(input_streams)

[mantepse]

I have one other question, concerning equality of parents: should I be using "==" or "is"? For example, in stream.py, _subs_in_caches:

            for i0, i in enumerate(indices):
                c = s._cache[i]
                if self._coefficient_ring == self._base_ring:
                    if c.parent() == self._PF:
                        c = retract(subs(c, var, val))
                        if c.parent() is not self._base_ring:
                            good = m - i0 - 1
                elif c.parent() == self._U:
                    c = c.map_coefficients(lambda e: subs(e, var, val))
                    try:
                        c = c.map_coefficients(lambda e: retract(e), self._base_ring)
                    except TypeError:
                        good = m - i0 - 1
                s._cache[i] = c

I actually think that I did it backwards here. I understand that computations could in principle return an element with parent that is not identical to the one I'm expecting. Do we have a nice parent where == is frequently (or at least sometimes) different from is?

[tscrim]

Permutation (sub)groups are good examples, but I think their group algebras might not fall under this. Hom spaces between equal-but-not-identical parents also have the same behavior IIRC. I'm not sure how many rings fall under this though. I would probably just be consistent with using == and != for all of them in _subs_in_cache just to be safe, even if it is marginally slower. Running a check with an example to verify this doesn't noticeably change the computation time would be good to do.

[mantepse]

OK, I can do that.

Meanwhile, I discovered another oddity. There is a bug in LazyModuleElement.change_ring:

    def change_ring(self, ring):
        r"""
        Return ``self`` with coefficients converted to elements of ``ring``.
...
        """
        P = self.parent()
        if P._names is None:
            Q = type(P)(ring, sparse=P._sparse)
        else:
            Q = type(P)(ring, names=P.variable_names(), sparse=P._sparse)
        return Q.element_class(Q, self._coeff_stream)

This doesn't follow the specification. However, if I change this, I seem to get coercion errors. I have to double check, though.

Update: sorry, I think it actually does, because __getitem__ converts to the correct ring in the end. The problem with #37975 must be elsewhere.

[mantepse]

I'm afraid I found the problem. Consider, as an example, the lazy power series ring $B[[q, t]]$ over a base ring $B$. Let $B(\vec f)$ denote the fraction field of rational functions in variables $f_1, f_2, \dots$. Then it is probably unavoidable that we will sometimes multiply elements from $B(\vec f)$ and the polynomial ring $B[q, t]$, and, unfortunately, the result will live in $B[q, t](\vec f)$ instead of $B(\vec f)[q, t]$.

I admit I do not fully understand why this happens only rarely, and if so, it is unconsequential most of the time. Mostly, it seems to happen with exact series, where we do not make sure that their data (initial_coefficients and constant) live in the correct parent.

However, at least sometimes it is sheer luck that everything works out. For example, in LazyCauchyProductSeries.__pow__ we happen to coerce the initial coefficients of an exact series to the proper parent even if n=1.

I am not sure how to proceed. Do you think it makes sense to have another meeting?

[mantepse]

Superseeded by #38108