Bug in factorization of simple symbolic expressions

[egourgoulhon]

As reported in this ask.sagemath question, we have currently (Sage 8.7.beta3):

sage: factor(2*exp(x) + exp(-x))
3*e^x

as well as

sage: factor(x*exp(-x) + exp(-x))
(x + 1)*e^x

Another example of erroneous result exhibited in https://groups.google.com/d/msg/sage-devel/ytLqIb4soLw/c14ZKGqcAAAJ is

sage: var('t u')
(t, u)
sage: L = (u + t)/(u - t)
sage: factor(L.substitute(t=sqrt(u)))
(u + 1)/(u - 1)

This bug is there since at least Sage 8.4.beta4. It is not in Sage 8.3, hence it must have been introduced between 8.4.beta0 and 8.4.beta4. It seems to have been introduced in #23835 (see comment:5).

CC: @rwst @sagetrac-bpage

Component: symbolics

Keywords: factor, exponential, pynac

Author: Benjamin Hackl

Branch/Commit: 7c2432c

Reviewer: Emmanuel Charpentier

Issue created by migration from https://trac.sagemath.org/ticket/27304

[EmmanuelCharpentier]

comment:2

in the same vein :

sage: e^(x+2*I*pi)
cosh(-x) - sinh(-x)

which is true, but largely not what we seek...

sage: e^(x+2*I*pi)==e^x
cosh(-x) - sinh(-x) == e^x
sage: bool(e^(x+2*I*pi)==e^x)
False

which is false. But :

sage: e^(x+2*I*pi).maxima_methods().exponentialize()
cosh(-x) - sinh(-x)
sage: (e^(x+2*I*pi)).maxima_methods().exponentialize()
e^x

Deserves to be critical, if not blocker.

[egourgoulhon]

comment:3

A difference between Sage 8.3 (where the factorization works) and Sage 8.4.beta4 is that in Sage 8.3 the factorization is performed by Maxima: cf. these lines of src/sage/symbolic/expression.pyx

    def factor(self, dontfactor=[]):
        from sage.calculus.calculus import symbolic_expression_from_maxima_string, symbolic_expression_from_string
        if len(dontfactor) > 0:
            m = self._maxima_()
            name = m.name()
            varstr = ','.join(['_SAGE_VAR_'+str(v) for v in dontfactor])
            cmd = 'block([dontfactor:[%s]],factor(%s))'%(varstr, name)
            return symbolic_expression_from_maxima_string(cmd)
        else:
            try:
                from sage.rings.all import QQ
                f = self.polynomial(QQ)
                w = repr(f.factor())
                return symbolic_expression_from_string(w)
            except (TypeError, NotImplementedError):
                pass
            return self.parent()(self._maxima_().factor())

while in Sage 8.4.beta4 and later, it is performed by pynac (via the function g_factor):

     def factor(self, dontfactor=[]):
        from sage.calculus.calculus import symbolic_expression_from_maxima_string
        cdef GEx x
        cdef bint b
        if dontfactor:
            m = self._maxima_()
            name = m.name()
            varstr = ','.join(['_SAGE_VAR_' + str(v) for v in dontfactor])
            cmd = 'block([dontfactor:[%s]],factor(%s))' % (varstr, name)
            return symbolic_expression_from_maxima_string(cmd)
        sig_on()
        try:
            b = g_factor(self._gobj, x)
        finally:
            sig_off()
        if b:
            return new_Expression_from_GEx(self._parent, x)
        else:
            return self

Maxima is correct there: in Sage 8.7.beta3, we have

./sage -maxima
(%i1) factor(2*exp(x) + exp(-x));     
                                - x      2 x
(%o1)                         %e    (2 %e    + 1)

So I would say it is a pynac bug.

[egourgoulhon]

Changed keywords from factor, exponential to factor, exponential, pynac

[egourgoulhon]

comment:5

The change from Maxima factorization to Pynac one, which triggered the bug, was introduced in #23835, which has been merged in Sage 8.4.beta3.

[EmmanuelCharpentier]

comment:7

This sage-support post might be related...

[embray]

comment:8

Moving all blocker/critical issues from 8.7 to 8.8.

[behackl]

comment:12

Would it be an option to let pynac handle factorization only when called on a rational function and call maxima (just as it was before #23835) otherwise?

[EmmanuelCharpentier]

comment:13

Replying to @behackl:

Would it be an option to let pynac handle factorization only when called on a rational function and call maxima (just as it was before #23835) otherwise?

Seconded :

If that means reverting the changes brought in by #23835, and if that does not break other changes intervening since, this would be a reasonable option (even if bad for performance), and would give people who know what they're talking about time to poner our next move...

[egourgoulhon]

comment:14

Replying to @EmmanuelCharpentier:

Replying to @behackl:

Would it be an option to let pynac handle factorization only when called on a rational function and call maxima (just as it was before #23835) otherwise?

Seconded :

Idem, +1. This sounds a good strategy for the short term. It would indeed be a pity if Sage 8.8 is shipped with such a bug.
Benjamin, do you intend to take this in charge?

[behackl]

comment:15

Replying to @egourgoulhon:

Idem, +1. This sounds a good strategy for the short term. It would indeed be a pity if Sage 8.8 is shipped with such a bug.
Benjamin, do you intend to take this in charge?

I already thought a bit about it, so yes, I will prepare a branch with the (short-term) bugfix I proposed above.

[egourgoulhon]

comment:16

Replying to @behackl:

I already thought a bit about it, so yes, I will prepare a branch with the (short-term) bugfix I proposed above.

Great!

[behackl]

Commit: 7c2432c

[behackl]

comment:17

This is what I would suggest. I'm entirely open for suggestions regarding the name as well as the code for is_rational_expression.

Locally, I only checked the doctests in expression.pyx, maybe a patchbot finds further occurrences of the symbolic factor method where the output has to be changed.

---

New commits:

0618fd3	helper method: is_rational_expression
131cd84	let maxima handle factorization of non-rational expressions
7c2432c	fix doctest; add more tests

[behackl]

Author: Benjamin Hackl

[behackl]

Branch: u/behackl/simple-factorization-rational

[EmmanuelCharpentier]

Reviewer: Emmanuel Charpentier

[EmmanuelCharpentier]

comment:18

On top of 8.8.beta5+#27738, builds and passes ptestlong exept for a transient (and already reported twice):

----------------------------------------------------------------------
sage -t --long --warn-long 151.5 src/sage/tests/books/computational-mathematics-with-sagemath/graphique_doctest.py  # 1 doctest failed
----------------------------------------------------------------------

which passes when ran standalone.

==> positive_review.

[egourgoulhon]

comment:19

Thank you Benjamin!

[vbraun]

Changed branch from u/behackl/simple-factorization-rational to 7c2432c