Trac #27958: enhance the integration by using also giac and sympy

URL: https://trac.sagemath.org/27958
Reported by: chapoton
Ticket author(s): Frédéric Chapoton
Reviewer(s): Thierry Monteil

src/sage/calculus/functional.py

@@ -207,12 +207,12 @@ def integral(f, *args, **kwds):
         0
         sage: restore('x,y')   # restore the symbolic variables x and y
 
-    Sage is unable to do anything with the following integral::
+    Sage is now (:trac:`27958`) able to compute the following integral::
 
-        sage: integral( exp(-x^2)*log(x), x )
-        integrate(e^(-x^2)*log(x), x)
+        sage: integral(exp(-x^2)*log(x), x)
+        1/2*sqrt(pi)*erf(x)*log(x) - x*hypergeometric((1/2, 1/2), (3/2, 3/2), -x^2)
 
-    Note, however, that::
+    and its value::
 
         sage: integral( exp(-x^2)*ln(x), x, 0, oo)
         -1/4*sqrt(pi)*(euler_gamma + 2*log(2))
@@ -222,7 +222,7 @@ def integral(f, *args, **kwds):
         sage: integral( ln(x)/x, x, 1, 2)
         1/2*log(2)^2
 
-    Sage can't do this elliptic integral (yet)::
+    Sage cannot do this elliptic integral (yet)::
 
         sage: integral(1/sqrt(2*t^4 - 3*t^2 - 2), t, 2, 3)
         integrate(1/sqrt(2*t^4 - 3*t^2 - 2), t, 2, 3)

src/sage/calculus/tests.py

@@ -122,7 +122,7 @@
     1/16*sqrt(pi)*((I + 1)*sqrt(2)*erf((1/2*I + 1/2)*sqrt(2)*x) + (I - 1)*sqrt(2)*erf((1/2*I - 1/2)*sqrt(2)*x) - (I - 1)*sqrt(2)*erf(sqrt(-I)*x) + (I + 1)*sqrt(2)*erf((-1)^(1/4)*x))
 
     sage: integrate((1-x^2)^n,x)
-    integrate((-x^2 + 1)^n, x)
+    x*hypergeometric((1/2, -n), (3/2,), x^2*exp_polar(2*I*pi))
     sage: integrate(x^x,x)
     integrate(x^x, x)
     sage: integrate(1/(x^3+1),x)
@@ -139,24 +139,20 @@
     sqrt(pi)/sqrt(c)
     sage: forget()
 
-The following are a bunch of examples of integrals that Mathematica
-can do, but Sage currently can't do::
+Other examples that now (:trac:`27958`) work::
 
-    sage: integrate(log(x)*exp(-x^2), x)        # todo -- Mathematica can do this
-    integrate(e^(-x^2)*log(x), x)
-
-Todo - Mathematica can do this and gets `\pi^2/15`.
-
-::
+    sage: integrate(log(x)*exp(-x^2), x)
+    1/2*sqrt(pi)*erf(x)*log(x) - x*hypergeometric((1/2, 1/2), (3/2, 3/2), -x^2)
 
     sage: integrate(log(1+sqrt(1+4*x)/2)/x, x, 0, 1)
     Traceback (most recent call last):
     ...
     ValueError: Integral is divergent.
 
-::
+The following is an example of integral that Mathematica
+can do, but Sage currently cannot do::
 
-    sage: integrate(ceil(x^2 + floor(x)), x, 0, 5)    # todo: Mathematica can do this
+    sage: integrate(ceil(x^2 + floor(x)), x, 0, 5, algorithm='maxima')
     integrate(ceil(x^2) + floor(x), x, 0, 5)
 
 MAPLE: The basic differentiation and integration examples in the
@@ -203,8 +199,8 @@
     1/3*sqrt(3)*arctan(1/3*sqrt(3)*(2*x + 1)) - 1/6*log(x^2 + x + 1) + 1/3*log(x - 1)
     sage: integrate(exp(-x^2), x)
     1/2*sqrt(pi)*erf(x)
-    sage: integrate(exp(-x^2)*log(x), x)       # todo: maple can compute this exactly.
-    integrate(e^(-x^2)*log(x), x)
+    sage: integrate(exp(-x^2)*log(x), x)
+    1/2*sqrt(pi)*erf(x)*log(x) - x*hypergeometric((1/2, 1/2), (3/2, 3/2), -x^2)
     sage: f = exp(-x^2)*log(x)
     sage: f.nintegral(x, 0, 999)
     (-0.87005772672831..., 7.5584...e-10, 567, 0)

src/sage/functions/exp_integral.py

@@ -436,7 +436,8 @@ def __init__(self):
             log_integral(3)
             sage: log_integral(x)._sympy_()
             li(x)
-
+            sage: log_integral(x)._fricas_init_()
+            'li(x)'
         """
         BuiltinFunction.__init__(self, "log_integral", nargs=1,
                                  latex_name=r'log_integral',
@@ -793,13 +794,16 @@ def __init__(self):
             sin_integral(1)
             sage: sin_integral(x)._sympy_()
             Si(x)
-
+            sage: sin_integral(x)._fricas_init_()
+            'Si(x)'
+            sage: sin_integral(x)._giac_()
+            Si(x)
         """
         BuiltinFunction.__init__(self, "sin_integral", nargs=1,
                                  latex_name=r'\operatorname{Si}',
                                  conversions=dict(maxima='expintegral_si',
                                                   sympy='Si',
-                                                  fricas='Si'))
+                                                  fricas='Si', giac='Si'))
 
     def _eval_(self, z):
         """
@@ -909,7 +913,8 @@ class Function_cos_integral(BuiltinFunction):
     Compare ``cos_integral(3.0)`` to the definition of the value using
     numerical integration::
 
-        sage: N(euler_gamma + log(3.0) + integrate((cos(x)-1)/x, x, 0, 3.0) - cos_integral(3.0)) < 1e-14
+        sage: a = numerical_integral((cos(x)-1)/x, 0, 3)[0]
+        sage: abs(N(euler_gamma + log(3)) + a - N(cos_integral(3.0))) < 1e-14
         True
 
     Arbitrary precision and complex arguments are handled::
@@ -965,13 +970,16 @@ def __init__(self):
             cos_integral(1)
             sage: cos_integral(x)._sympy_()
             Ci(x)
-
+            sage: cos_integral(x)._fricas_init_()
+            'Ci(x)'
+            sage: cos_integral(x)._giac_()
+            Ci(x)
         """
         BuiltinFunction.__init__(self, "cos_integral", nargs=1,
                                  latex_name=r'\operatorname{Ci}',
                                  conversions=dict(maxima='expintegral_ci',
                                                   sympy='Ci',
-                                                  fricas='Ci'))
+                                                  fricas='Ci', giac='Ci'))
 
     def _evalf_(self, z, parent=None, algorithm=None):
         """
@@ -1038,7 +1046,8 @@ class Function_sinh_integral(BuiltinFunction):
     Compare ``sinh_integral(3.0)`` to the definition of the value using
     numerical integration::
 
-        sage: N(integrate((sinh(x))/x, x, 0, 3.0) - sinh_integral(3.0)) < 1e-14
+        sage: a = numerical_integral(sinh(x)/x, 0, 3)[0]
+        sage: abs(a - N(sinh_integral(3))) < 1e-14
         True
 
     Arbitrary precision and complex arguments are handled::
@@ -1197,8 +1206,8 @@ class Function_cosh_integral(BuiltinFunction):
     Compare ``cosh_integral(3.0)`` to the definition of the value using
     numerical integration::
 
-        sage: N(euler_gamma + log(3.0) + integrate((cosh(x)-1)/x, x, 0, 3.0) -
-        ....:   cosh_integral(3.0)) < 1e-14
+        sage: a = numerical_integral((cosh(x)-1)/x, 0, 3)[0]
+        sage: abs(N(euler_gamma + log(3)) + a - N(cosh_integral(3.0))) < 1e-14
         True
 
     Arbitrary precision and complex arguments are handled::

src/sage/functions/generalized.py

@@ -391,10 +391,14 @@ class FunctionSignum(BuiltinFunction):
 
     TESTS:
 
-    Check if conversion to sympy works :trac:`11921`::
+    Check if conversions to sympy and others work (:trac:`11921`)::
 
         sage: sgn(x)._sympy_()
         sign(x)
+        sage: sgn(x)._fricas_init_()
+        'sign(x)'
+        sage: sgn(x)._giac_()
+        sign(x)
 
     REFERENCES:
 
@@ -419,7 +423,8 @@ def __init__(self):
             sign(x)
         """
         BuiltinFunction.__init__(self, "sgn", latex_name=r"\mathrm{sgn}",
-                conversions=dict(maxima='signum',mathematica='Sign',sympy='sign'),
+                conversions=dict(maxima='signum', mathematica='Sign',
+                                 sympy='sign', giac='sign', fricas='sign'),
                 alt_name="sign")
 
     def _eval_(self, x):

src/sage/functions/min_max.py

@@ -249,12 +249,14 @@ def _evalf_(self, *args, **kwds):
         ::
 
             sage: f = max_symbolic(sin(x), cos(x))
-            sage: r = integral(f, x, 0, 1)
+            sage: r = integral(f, x, 0, 1); r
+            sqrt(2) - cos(1)
             sage: r.n()
-            0.8739124411567263
+            0.873911256504955
         """
         return max_symbolic(args)
 
+
 max_symbolic = MaxSymbolic()
 
 

src/sage/symbolic/integration/integral.py

@@ -61,7 +61,9 @@ def __init__(self):
         # in the given order. This is an attribute of the class instead of
         # a global variable in this module to enable customization by
         # creating a subclasses which define a different set of integrators
-        self.integrators = [external.maxima_integrator]
+        self.integrators = [external.maxima_integrator,
+                            external.giac_integrator,
+                            external.sympy_integrator]
 
         BuiltinFunction.__init__(self, "integrate", nargs=2, conversions={'sympy': 'Integral',
                                                                           'giac': 'integrate'})
@@ -86,12 +88,21 @@ def _eval_(self, f, x):
                 return None
 
         # we try all listed integration algorithms
+        A = None
         for integrator in self.integrators:
             try:
-                return integrator(f, x)
-            except NotImplementedError:
+                A = integrator(f, x)
+            except (NotImplementedError, TypeError):
                 pass
-        return None
+            except ValueError:
+                # maxima is telling us something
+                raise
+            else:
+                if not hasattr(A, 'operator'):
+                    return A
+                elif not isinstance(A.operator(), IndefiniteIntegral):
+                    return A
+        return A
 
     def _tderivative_(self, f, x, diff_param=None):
         """
@@ -151,7 +162,9 @@ def __init__(self):
         # in the given order. This is an attribute of the class instead of
         # a global variable in this module to enable customization by
         # creating a subclasses which define a different set of integrators
-        self.integrators = [external.maxima_integrator]
+        self.integrators = [external.maxima_integrator,
+                            external.giac_integrator,
+                            external.sympy_integrator]
 
         BuiltinFunction.__init__(self, "integrate", nargs=4, conversions={'sympy': 'Integral',
                                                                           'giac': 'integrate'})
@@ -178,12 +191,21 @@ def _eval_(self, f, x, a, b):
         args = (f, x, a, b)
 
         # we try all listed integration algorithms
+        A = None
         for integrator in self.integrators:
             try:
-                return integrator(*args)
-            except NotImplementedError:
+                A = integrator(*args)
+            except (NotImplementedError, TypeError):
                 pass
-        return None
+            except ValueError:
+                # maxima is telling us something
+                raise
+            else:
+                if not hasattr(A, 'operator'):
+                    return A
+                elif not isinstance(A.operator(), DefiniteIntegral):
+                    return A
+        return A
 
     def _evalf_(self, f, x, a, b, parent=None, algorithm=None):
         """
@@ -552,20 +574,28 @@ def integrate(expression, v=None, a=None, b=None, algorithm=None, hold=False):
         sage: integral(e^(-x^2),(x, 0, 0.1))
         0.05623145800914245*sqrt(pi)
 
-    An example of an integral that fricas can integrate, but the
-    default integrator cannot::
+    An example of an integral that fricas can integrate::
 
         sage: f(x) = sqrt(x+sqrt(1+x^2))/x
         sage: integrate(f(x), x, algorithm="fricas")      # optional - fricas
         2*sqrt(x + sqrt(x^2 + 1)) - 2*arctan(sqrt(x + sqrt(x^2 + 1))) - log(sqrt(x + sqrt(x^2 + 1)) + 1) + log(sqrt(x + sqrt(x^2 + 1)) - 1)
 
-    The following definite integral is not found with the
-    default integrator::
+    where the default integrator obtains another answer::
+
+        sage: integrate(f(x), x)
+        1/8*sqrt(x)*gamma(1/4)*gamma(-1/4)^2*hypergeometric((-1/4, -1/4, 1/4), (1/2, 3/4), -1/x^2)/(pi*gamma(3/4))
+
+    The following definite integral is not found by maxima::
 
         sage: f(x) = (x^4 - 3*x^2 + 6) / (x^6 - 5*x^4 + 5*x^2 + 4)
-        sage: integrate(f(x), x, 1, 2)
+        sage: integrate(f(x), x, 1, 2, algorithm='maxima')
         integrate((x^4 - 3*x^2 + 6)/(x^6 - 5*x^4 + 5*x^2 + 4), x, 1, 2)
 
+    but is nevertheless computed::
+
+        sage: integrate(f(x), x, 1, 2)
+        -1/2*pi + arctan(8) + arctan(5) + arctan(2) + arctan(1/2)
+
     Both fricas and sympy give the correct result::
 
         sage: integrate(f(x), x, 1, 2, algorithm="fricas")  # optional - fricas
@@ -577,6 +607,8 @@ def integrate(expression, v=None, a=None, b=None, algorithm=None, hold=False):
 
         sage: integrate(abs(cos(x)), x, 0, 2*pi, algorithm='giac')
         4
+        sage: integrate(abs(cos(x)), x, 0, 2*pi)
+        4
 
     ALIASES: integral() and integrate() are the same.
 
@@ -718,12 +750,12 @@ def integrate(expression, v=None, a=None, b=None, algorithm=None, hold=False):
         sage: error.numerical_approx() # abs tol 1e-10
         0
 
-    We will not get an evaluated answer here, which is better than
+    We will get a correct answer here, which is better than
     the previous (wrong) answer of zero. See :trac:`10914`::
 
         sage: f = abs(sin(x))
-        sage: integrate(f, x, 0, 2*pi)  # long time (4s on sage.math, 2012)
-        integrate(abs(sin(x)), x, 0, 2*pi)
+        sage: integrate(f, x, 0, 2*pi)
+        4
 
     Another incorrect integral fixed upstream in Maxima, from
     :trac:`11233`::
@@ -815,6 +847,43 @@ def integrate(expression, v=None, a=None, b=None, algorithm=None, hold=False):
         sage: f = (x^2)*exp(x) / (1+exp(x))^2
         sage: integrate(f, (x, -infinity, infinity))
         1/3*pi^2
+
+    Some integrals are now working (:trac:`27958`, using giac or sympy)::
+
+        sage: integrate(1/(1 + abs(x)), x)
+        log(abs(x*sgn(x) + 1))/sgn(x)
+
+        sage: integrate(cos(x + abs(x)), x)
+        sin(x*sgn(x) + x)/(sgn(x) + 1)
+
+        sage: integrate(abs(x^2 - 1), x, -2, 2)
+        4
+
+        sage: f = sqrt(x + 1/x^2)
+        sage: integrate(f, x)
+        1/3*(2*sqrt(x^3 + 1) - log(sqrt(x^3 + 1) + 1) + log(abs(sqrt(x^3 + 1) - 1)))*sgn(x)
+
+        sage: g = abs(sin(x)*cos(x))
+        sage: g.integrate(x, 0, 2*pi)
+        2
+
+        sage: integrate(1/sqrt(abs(x)), x)
+        2*sqrt(x*sgn(x))/sgn(x)
+
+        sage: integrate(sgn(x) - sgn(1-x), x)
+        x*(sgn(x) - sgn(-x + 1)) + sgn(-x + 1)
+
+        sage: integrate(1 / (1 + abs(x-5)), x, -5, 6)
+        log(11) + log(2)
+
+        sage: integrate(1/(1 + abs(x)), x)
+        log(abs(x*sgn(x) + 1))/sgn(x)
+
+        sage: integrate(cos(x + abs(x)), x)
+        sin(x*sgn(x) + x)/(sgn(x) + 1)
+
+        sage: integrate(abs(x^2 - 1), x, -2, 2)
+        4
     """
     expression, v, a, b = _normalize_integral_input(expression, v, a, b)
     if algorithm is not None: