A computation straightforward in Mathematica and almost impossible in sage

videlec commented 1 year ago

During an admcycles workshop we wanted to compute the inverse (in terms of composition) of the power series in $\mathbb{Q}[p,q,\sqrt{p+q}][[z]]$ given by $$- \frac{1}{\sqrt{p + q}} \frac{p}{q} \left( 1 - \frac{1 - \exp \left( \frac{2 \cdot q}{\sqrt{p + q} } \cdot z \right)}{1 - \exp \left( - \ \frac{2 \cdot p}{\sqrt{p+q}} \cdot z \right)} \right)$$

Missing simplification when working over `QQ[p, q, u]`

sage: R.<p,q,u> = QQ[]
sage: P.<z> = PowerSeriesRing(R)
sage: S = (1 - (2 * q / u * z).exp()) / (1 - (- 2 * p / u * z).exp())
sage: S[1]
(-28269553036454149273332760011886696253239742350009903329945699220681916416*p*q - 28269553036454149273332760011886696253239742350009903329945699220681916416*q^2)/(28269553036454149273332760011886696253239742350009903329945699220681916416*p*u)

A "fix" is to work over ZZ[p, q, u] instead

sage: R.<p,q,u> = ZZ[]
sage: P.<z> = PowerSeriesRing(R)
sage: S = (1 - (2 * q / u * z).exp()) / (1 - (- 2 * p / u * z).exp())
sage: S[1]
(-p*q - q^2)/(p*u)

The fraction field of QQ[p, q, u] ought to be the one of ZZ.

Working in quotient `QQ[p, q, u] / (u^2 - p - q)` and `ZZ[p, q, u] / (u^2 - p - q)`

Working with ZZ[p, q, u][[z]] does not simplify enough

sage: S[3].numerator() / S[3].denominator().reduce(R.ideal(u^2 - p - q))
(-p*q^2 - q^3)/(3*p*u)
sage: S[3]
(-p^2*q^2 - 2*p*q^3 - q^4)/(3*p*u^3)

and it is hence tempting to work over quotients. Neither ZZ[p, q, u] / (u^2 - p - q) nor QQ[p, q, u] / (u^2 - p - q) is satisfying. The former is just broken with respect to arithmetic

sage: Rbar = R.quotient(u^2 - p - q, names='p,q,u')
sage: p, q, u = Rbar.gens()
sage: p / u
Traceback (most recent call last):
.../sage/structure/element.pyx in sage.structure.element.Element.__truediv__ (build/cythonized/sage/structure/element.c:13274)()
   1734         cdef int cl = classify_elements(left, right)
   1735         if HAVE_SAME_PARENT(cl):
-> 1736             return (<Element>left)._div_(right)
   1737         if BOTH_ARE_ELEMENT(cl):
   1738             return coercion_model.bin_op(left, right, truediv)

.../sage/structure/element.pyx in sage.structure.element.RingElement._div_ (build/cythonized/sage/structure/element.c:19235)()
   2728         return l
   2729 
-> 2730     cpdef _div_(self, other):
   2731         """
   2732         Default implementation of division using the fraction field.

.../sage/rings/quotient_ring_element.py in _div_(self, right)
    420         B  = I.groebner_basis()
    421         try:
--> 422             XY = L.lift((R,) + tuple(B))
    423         except ValueError:
    424              raise ArithmeticError("Division failed. The numerator is not "

.../sage/rings/polynomial/multi_polynomial_libsingular.pyx in sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular.lift (build/cythonized/sage/rings/polynomial/multi_polynomial_libsingular.cpp:35422)()
   4560         global errorreported
   4561         if not self._parent._base.is_field():
-> 4562             raise NotImplementedError("Lifting of multivariate polynomials over non-fields is not implemented.")
   4563 
   4564         cdef ideal *fI = idInit(1,1)

NotImplementedError: Lifting of multivariate polynomials over non-fields is not implemented.

While QQ[p,q,u] / (u^2 - p - q)[[z]] is still on strike with respect to simplification

age: R.<p,q,u> = QQ[]
sage: P.<z> = PowerSeriesRing(R)
sage: S = (1 - (2 * q / u * z).exp()) / (1 - (- 2 * p / u * z).exp())
sage: S[1]
(-28269553036454149273332760011886696253239742350009903329945699220681916416*p*q - 28269553036454149273332760011886696253239742350009903329945699220681916416*q^2)/(28269553036454149273332760011886696253239742350009903329945699220681916416*p*u)

Simple workaround

For the purpose of the computation given in the title one could simply use variables p, u and sets q = u^2 - p. However, this is not satisfactory as all the troubles mentioned in this issue are worth being fixed!

videlec commented 1 year ago

@mezzarobba You might have insight on how to make that works out of the box?

miguelmarco commented 1 year ago

you can get the simplification like this:

sage:  R.<p,q,u> = QQ[]
sage:  P.<z> = PowerSeriesRing(R)
sage:  S = (1 - (2 * q / u * z).exp()) / (1 - (- 2 * p / u * z).exp())
sage: S[0]
(-14134776518227074636666380005943348126619871175004951664972849610340958208*q)/(14134776518227074636666380005943348126619871175004951664972849610340958208*p)
sage: S[0].factor()
(-1) * p^-1 * q
sage: S[0].factor().prod()
(-q)/p

I guess the decision of not doing the simplification automatically was taken because of performance reasons (in this case, it is two order of magnitude slower than just computing the non-simplified expression).

videlec commented 1 year ago

you can get the simplification like this:
sage:  R.<p,q,u> = QQ[]
sage:  P.<z> = PowerSeriesRing(R)
sage:  S = (1 - (2 * q / u * z).exp()) / (1 - (- 2 * p / u * z).exp())
sage: S[0]
(-14134776518227074636666380005943348126619871175004951664972849610340958208*q)/(14134776518227074636666380005943348126619871175004951664972849610340958208*p)
sage: S[0].factor()
(-1) * p^-1 * q
sage: S[0].factor().prod()
(-q)/p
I guess the decision of not doing the simplification automatically was taken because of performance reasons (in this case, it is two order of magnitude slower than just computing the non-simplified expression).

Your comment discusses one aspect of the question : the coefficients. And, concerning this aspect, I don't find it satisfactory. Here the coefficients of your version of S[n] grows extremly fast and makes the computation unfeasible at large n (let say n=100) by first making the series expansion.

sage: list_plot([max([c.height().ndigits() for c in S[k].numerator().coefficients()]) for k in range(1,16)])

coeff_growth

For example

sage: P.<z> = PowerSeriesRing(R, default_prec=40)
sage: %time S = (1 - (2 * q / u * z).exp()) / (1 - (- 2 * p / u * z).exp())
CPU times: user 9.61 s, sys: 44.6 ms, total: 9.66 s
Wall time: 9.66 s

But note also that the question contains another aspect : the simplifications not automatically happening when taking $u = \sqrt{p + q}$.

mezzarobba commented 1 year ago

@mezzarobba You might have insight on how to make that works out of the box?

Well, the Sage implementation of fraction fields is a bit of a mess.

What is happening here as far as I understand is that fractions of multivariate polynomials use a generic implementation of fraction fields, which calls the gcd method of the numerator/denominator to simplify fractions. That is why constant common factors are removed from fractions of polynomials over ZZ but not over QQ.

For fractions of univariate polynomials, we have a dedicated subclass that does a bit better. But even that dedicated class normalizes fractions in a slightly strange way: for instance

sage: P.<t> = QQ[]
sage: 1/(2*t)
1/2/t

where 1/(2*t) would probably be more reasonable in most cases.

We certainly need some kind of mechanism to reduce Frac(Frac(R)[x,y,...]) to Frac(R[x,y,...]), but I am not sure exactly how that should work. I also wonder if the current normalization for univariate polynomials over fields needs to remain available as an alternative option even for polynomials over fraction fields.

As a shorter term fix, it shouldn't be hard to to something similar to what we do in the univariate case (that is, force the leading coefficient of the denominator to one) but for multivariate polynomials. It would be useful in any case for working with multivariate polynomials over fields that are not fraction fields.

Finally, in the specific case of multivariate polyomials over ℚ, we could also interface FLINT's (formerly Calcium's) fmpz_mpoly_q.

Related tickets include #16268, #26339.

sagemath / sage