Trac #10070: make heaviside and step play nicely together.

The heaviside built-in function and the step symbolic method are not the
same at zero.
{{{
sage: heaviside(0)
heaviside(0)
sage: SR(0).step()
1/2
}}}
And the documentation seems to indicate that `heaviside` should be
undefined (?) at zero, though it's not definitive.

In addition to reconciling these, we should probably unify notation or
something.

URL: https://trac.sagemath.org/10070
Reported by: kcrisman
Ticket author(s): Ralf Stephan
Reviewer(s): Travis Scrimshaw, Marcelo Forets

Author: Release Manager

src/sage/functions/generalized.py

@@ -51,7 +51,7 @@
 #
 ##############################################################################
 
-from sage.symbolic.function import BuiltinFunction
+from sage.symbolic.function import (BuiltinFunction, GinacFunction)
 from sage.rings.all import ComplexIntervalField, ZZ
 
 class FunctionDiracDelta(BuiltinFunction):
@@ -166,7 +166,7 @@ def _evalf_(self, x, **kwds):
 
 dirac_delta = FunctionDiracDelta()
 
-class FunctionHeaviside(BuiltinFunction):
+class FunctionHeaviside(GinacFunction):
     r"""
     The Heaviside step function, `H(x)` (``heaviside(x)``).
 
@@ -180,6 +180,8 @@ class FunctionHeaviside(BuiltinFunction):
 
         `H(x) = 0` for `x < 0` and `H(x) = 1` for `x > 0`
 
+    .. SEEALSO:: :func:`unit_step()<sage.functions.generalized.FunctionUnitStep>`
+
     EXAMPLES::
 
         sage: heaviside(-1)
@@ -191,10 +193,29 @@ class FunctionHeaviside(BuiltinFunction):
         sage: heaviside(x)
         heaviside(x)
 
+        sage: heaviside(-1/2)
+        0
+        sage: heaviside(exp(-1000000000000000000000))
+        1
+
     TESTS::
 
         sage: heaviside(x)._sympy_()
         Heaviside(x)
+        sage: heaviside(x).subs(x=1)
+        1
+        sage: heaviside(x).subs(x=-1)
+        0
+
+    ::
+
+        sage: ex = heaviside(x)+1
+        sage: t = loads(dumps(ex)); t
+        heaviside(x) + 1
+        sage: bool(t == ex)
+        True
+        sage: t.subs(x=1)
+        2
 
     REFERENCES:
 
@@ -225,74 +246,16 @@ def __init__(self):
             Heaviside(x)
             sage: heaviside(x)._giac_()
             Heaviside(x)
+            sage: h(x) = heaviside(x)
+            sage: h(pi).numerical_approx()
+            1.00000000000000
         """
-        BuiltinFunction.__init__(self, "heaviside", latex_name="H",
+        GinacFunction.__init__(self, "heaviside", latex_name="H",
                                  conversions=dict(maxima='hstep',
                                                   mathematica='HeavisideTheta',
                                                   sympy='Heaviside',
                                                   giac='Heaviside'))
 
-    def _eval_(self, x):
-        """
-        INPUT:
-
-        -  ``x`` - a real number or a symbolic expression
-
-        EXAMPLES::
-
-            sage: heaviside(-1/2)
-            0
-            sage: heaviside(1)
-            1
-            sage: heaviside(0)
-            heaviside(0)
-            sage: heaviside(x)
-            heaviside(x)
-            sage: heaviside(exp(-1000000000000000000000))
-            1
-
-        Evaluation test::
-
-            sage: heaviside(x).subs(x=1)
-            1
-            sage: heaviside(x).subs(x=-1)
-            0
-
-        ::
-
-            sage: ex = heaviside(x)+1
-            sage: t = loads(dumps(ex)); t
-            heaviside(x) + 1
-            sage: bool(t == ex)
-            True
-            sage: t.subs(x=1)
-            2
-        """
-        try:
-            return self._evalf_(x)
-        except (TypeError,ValueError):      # x is symbolic
-            pass
-        return None
-
-    def _evalf_(self, x, **kwds):
-        """
-        TESTS::
-
-            sage: h(x) = heaviside(x)
-            sage: h(pi).numerical_approx()
-            1.00000000000000
-        """
-        approx_x = ComplexIntervalField()(x)
-        if bool(approx_x.imag() == 0):      # x is real
-            if bool(approx_x.real() == 0):  # x is zero
-                return None
-            # Now we have a non-zero real
-            if bool((approx_x**(0.5)).imag() == 0): # Check: x > 0
-                return 1
-            else:
-                return 0
-        raise ValueError("Numeric evaluation of symbolic expression")
-
     def _derivative_(self, x, diff_param=None):
         """
         Derivative of Heaviside step function
@@ -306,7 +269,7 @@ def _derivative_(self, x, diff_param=None):
 
 heaviside = FunctionHeaviside()
 
-class FunctionUnitStep(BuiltinFunction):
+class FunctionUnitStep(GinacFunction):
     r"""
     The unit step function, `\mathrm{u}(x)` (``unit_step(x)``).
 
@@ -320,6 +283,8 @@ class FunctionUnitStep(BuiltinFunction):
 
         `\mathrm{u}(x) = 0` for `x < 0` and `\mathrm{u}(x) = 1` for `x \geq 0`
 
+    .. SEEALSO:: :func:`heaviside()<sage.functions.generalized.FunctionHeaviside>`
+
     EXAMPLES::
 
         sage: unit_step(-1)
@@ -330,6 +295,18 @@ class FunctionUnitStep(BuiltinFunction):
         1
         sage: unit_step(x)
         unit_step(x)
+        sage: unit_step(-exp(-10000000000000000000))
+        0
+
+    TESTS::
+
+        sage: unit_step(x).subs(x=1)
+        1
+        sage: unit_step(x).subs(x=0)
+        1
+        sage: h(x) = unit_step(x)
+        sage: h(pi).numerical_approx()
+        1.00000000000000
     """
     def __init__(self):
         r"""
@@ -359,60 +336,9 @@ def __init__(self):
             sage: t.subs(x=0)
             2
         """
-        BuiltinFunction.__init__(self, "unit_step", latex_name=r"\mathrm{u}",
+        GinacFunction.__init__(self, "unit_step", latex_name=r"\mathrm{u}",
                                    conversions=dict(mathematica='UnitStep'))
 
-    def _eval_(self, x):
-        """
-        INPUT:
-
-        -  ``x`` - a real number or a symbolic expression
-
-        EXAMPLES::
-
-            sage: unit_step(-1)
-            0
-            sage: unit_step(1)
-            1
-            sage: unit_step(0)
-            1
-            sage: unit_step(x)
-            unit_step(x)
-            sage: unit_step(-exp(-10000000000000000000))
-            0
-
-        Evaluation test::
-
-            sage: unit_step(x).subs(x=1)
-            1
-            sage: unit_step(x).subs(x=0)
-            1
-        """
-        try:
-            return self._evalf_(x)
-        except (TypeError,ValueError):      # x is symbolic
-            pass
-        return None
-
-    def _evalf_(self, x, **kwds):
-        """
-        TESTS::
-
-            sage: h(x) = unit_step(x)
-            sage: h(pi).numerical_approx()
-            1.00000000000000
-        """
-        approx_x = ComplexIntervalField()(x)
-        if bool(approx_x.imag() == 0):      # x is real
-            if bool(approx_x.real() == 0):  # x is zero
-                return 1
-            # Now we have a non-zero real
-            if bool((approx_x**(0.5)).imag() == 0): # Check: x > 0
-                return 1
-            else:
-                return 0
-        raise ValueError("Numeric evaluation of symbolic expression")
-
     def _derivative_(self, x, diff_param=None):
         """
         Derivative of unit step function

src/sage/symbolic/expression.pyx

@@ -7165,9 +7165,13 @@ cdef class Expression(CommutativeRingElement):
 
     def step(self, hold=False):
         """
-        Return the value of the Heaviside step function, which is 0 for
+        Return the value of the unit step function, which is 0 for
         negative x, 1 for 0, and 1 for positive x.
 
+        .. SEEALSO::
+
+            :class:`sage.functions.generalized.FunctionUnitStep`
+
         EXAMPLES::
 
             sage: x = var('x')