lazy series/pushout experiment

[mantepse]

This pull request provides a method define_implicitly for lazy rings that allows to produce lazy power series solutions for systems of functional equations, in relatively great generality, superseding #37033.

We use a special functor to make sure that coercion between CoefficientRing and its base ring works as desired. There is an alternative idea, using an action, which I did not explore. However, it should be a local change, if the current approach eventually shows its shortcomings.

Dependencies: #38729

[mantepse]

There is a time-sink I don't completely understand:

sage: h = SymmetricFunctions(QQ).h()
sage: L.<t, u> = LazyPowerSeriesRing(h.fraction_field())
sage: D = L.undefined()
sage: s1 = L.sum(lambda n: h[n]*t^(n+1)*u^(n-1), 1)
sage: L.define_implicitly([D], [u*D - u - u*s1*D - t*(D - D(t, 0))])
sage: D[0]
h[]
sage: D[1]
0
sage: D[2]
h[1]*t^2
sage: %prun D[3]
         3357653 function calls (3326279 primitive calls) in 3.081 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
[...]
  777/111    0.038    0.000    2.599    0.023 {method 'gcd' of 'sage.rings.polynomial.multi_polynomial.MPolynomial' objects}
[...]
  793/110    0.012    0.000    2.226    0.020 unique_factorization_domains.py:134(_gcd_univariate_polynomial)
[...]
sage: D[3]
0

[mantepse]

That time sink is also present with L.<t, u> = LazyPowerSeriesRing(QQbar["x"].fraction_field()), so I guess it cannot be easily fixed.

[mantepse]

Notes to myself, concerning the time sink. If I understand correctly, it is dominated by the code in Stream_uninitialized.__getitem__ that creates the new undetermined element of the stream.

I tried a bit to optimize this in the special case where isinstance(self._coefficient_ring, MPolynomialRing_polydict), mainly by bypassing MPolynomialRing_polydict.__call__ and the like. Curiously, it turned out that creating the string representation of the monomials takes a lot of time (it invokes multiplication and gcd!). However, even with all these optimizations, the performance was still quite bad. Essentially, we are using a lot of time in gcd computations, and it is not clear which of these can be avoided.

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
...
  1886         2       1883.0    941.5      0.0          for n0 in range(len(self._cache) + self._approximate_order, n+1):
  1887                                                       # WARNING: coercing the new variable to self._PF slows
  1888                                                       # down the multiplication enormously
  1889         1      12524.0  12524.0      0.0              if self._coefficient_ring == self._base_ring:
  1890                                                           x = (self._pool.new_variable(self._name + "[%s]" % n0)
  1891                                                                * self._terms_of_degree(n0, self._P)[0])
  1892                                                       else:
  1893         2   26085266.0    1e+07     11.6                  x = sum(self._pool.new_variable(self._name + "[%s]" % m) * m
  1894         1    4885896.0    5e+06      2.2                          for m in self._terms_of_degree(n0, self._P))
  1895         1  194465377.0    2e+08     86.2              x = self._U(x)
  1896         1       1333.0   1333.0      0.0              self._cache.append(x)
  1897                                           
  1898         1       2886.0   2886.0      0.0          return x

[mantepse]

Hm, somebody is calling Rings.ParentMethods.__getitem__ an awful lot.

[mantepse]

Oh wow. I extracted the bit. Creating FESDUMMY_13/h[]*t^4 + FESDUMMY_12/h[]*t^3*u + FESDUMMY_11/h[]*t^2*u^2 + FESDUMMY_10/h[]*t*u^3 + FESDUMMY_9/h[]*u^4 takes more than a second on my computer!

def testit(n0, _P, _U, _terms_of_degree):
    """
    sage: from sage.data_structures.stream import CoefficientRing, VariablePool
    sage: _base_ring = SymmetricFunctions(QQ).h().fraction_field()
    sage: _coefficient_ring = _base_ring["t", "u"]
    sage: _PF = CoefficientRing(_base_ring)
    sage: _U = _coefficient_ring.change_ring(_PF); _U
    sage: _P = _PF.base(); _P
    sage: _terms_of_degree = lambda n, R: [m.change_ring(R) for m in _coefficient_ring.monomials_of_degree(n)]
    sage: %time testit(4, _P, _U, _terms_of_degree)
    CPU times: user 1.3 s, sys: 4.02 ms, total: 1.3 s
    Wall time: 1.3 s
    FESDUMMY_13/h[]*t^4 + FESDUMMY_12/h[]*t^3*u + FESDUMMY_11/h[]*t^2*u^2 + FESDUMMY_10/h[]*t*u^3 + FESDUMMY_9/h[]*u^4
    """
    _pool = VariablePool(_P)
    x = sum(_pool.new_variable("[%s]" % m) * m for m in _terms_of_degree(n0, _P))
    return _U(x)

[mantepse]

Apparently, there is a hidden cost in the dense implementation of the infinite polynomial ring, I'll ask a question on sage-devel.

[mantepse]

I now have an experimental version of the gcd implementation which reduces the time to compute D above from 30 seconds to 4 seconds!