ArcTanh, ArcCoth: Indeterminate where the error covers the pole, as in WMA

mathics/builtin/numbers/hyperbolic.py

@@ -103,6 +103,7 @@ class ArcCoth(MPMathFunction):
     summary_text = "inverse hyperbolic cotangent function"
     sympy_name = "acoth"
     mpmath_name = "acoth"
+    indeterminate_at = (1, -1)
 
     @staticmethod
     def condition_numbers(x: mpmath.mpf, value: mpmath.mpf) -> tuple[mpmath.mpf]:
@@ -245,6 +246,10 @@ class ArcTanh(MPMathFunction):
      = 0
     >> ArcTanh[1]
      = Infinity
+
+    An inexact argument whose error covers the pole gives 'Indeterminate':
+    >> ArcTanh[1.0000001`5]
+     = Indeterminate
     >> ArcTanh[0]
      = 0
     >> ArcTanh[.5 + 2 I]
@@ -260,6 +265,7 @@ def condition_numbers(x: mpmath.mpf, value: mpmath.mpf) -> tuple[mpmath.mpf]:
         return (abs(x / ((1 - x**2) * value)),)
 
     numpy_name = "arctanh"
+    indeterminate_at = (1, -1)
 
     rules = {
         "ArcTanh[Undefined]": "Undefined",

mathics/core/builtin.py

@@ -95,6 +95,7 @@
 )
 from mathics.core.systemsymbols import (
     SymbolDefault,
+    SymbolIndeterminate,
     SymbolLessEqual,
     SymbolMessageName,
     SymbolRange,
@@ -103,7 +104,7 @@
 )
 from mathics.eval.arithmetic import eval_mpmath_function
 from mathics.eval.numbers.numbers import cancel
-from mathics.eval.numeric import eval_mpmath_function_by_condition
+from mathics.eval.numeric import covers_point, eval_mpmath_function_by_condition
 from mathics.eval.numerify import numerify
 from mathics.eval.scoping import dynamic_scoping
 from mathics.eval.sympy import eval_sympy
@@ -825,6 +826,9 @@ class MPMathFunction(SympyFunction):
     # no larger than the bound.
     has_poles: bool = False
     value_bound: Optional[float] = None
+    # Real points where, as in WMA, an inexact argument whose error covers
+    # one of them gives Indeterminate, like ArcTanh at 1.
+    indeterminate_at: tuple[int, ...] = ()
 
     def get_mpmath_function(self, args):
         if self.mpmath_name is None or len(args) not in self.nargs:
@@ -856,6 +860,9 @@ def eval(self, z, evaluation: Evaluation):
         if mpmath_function is None:
             return
 
+        if len(args) == 1 and covers_point(args[0], self.indeterminate_at):
+            return SymbolIndeterminate
+
         if self.condition_numbers is not None:
             result = eval_mpmath_function_by_condition(
                 mpmath_function,

mathics/eval/numeric.py

@@ -705,6 +705,26 @@ def eval_mpmath_function_by_condition(
     return _precision_real(value, precision10, accuracy10)
 
 
+def covers_point(number: Number, points: tuple[int, ...]) -> bool:
+    """
+    True if ``number`` is a nonzero inexact real number and one of
+    ``points`` lies within its error. As in WMA, ArcTanh[1.0000001`5]
+    is Indeterminate: its error 10^-5 covers the pole at 1, where f goes
+    to Infinity on one side and has an imaginary part on the other. A
+    MachineReal has no error, so only the point itself is covered.
+    """
+    if isinstance(number, PrecisionReal):
+        if number.is_zero:
+            return False
+        with mpmath.workprec(_stored_bits((number,)) + 10):
+            x = _full_mpmath(number)
+            error = abs(x) * mpmath.mpf(10) ** -number.precision10
+            return any(abs(x - point) <= error for point in points)
+    if isinstance(number, Real):
+        return number.to_mpmath() in points
+    return False
+
+
 def _change_within(
     mpmath_function: Callable, x: mpmath.mpf, delta: mpmath.mpf, digits: float
 ) -> Optional[tuple[mpmath.mpf, mpmath.mpf]]:

test/builtin/atomic/test_numbers.py

@@ -716,6 +716,19 @@ def test_arctrig_real_valued_precision(str_expr, str_expected, msg):
         ("Precision[ArcCosh[1.00000001`20]]", "12.301", "next to the end"),
         ("Precision[ArcSech[0.99999999`20]]", "12.301", None),
         ("Precision[ArcTanh[0.99999999`20]]", "13.2813", "next to a pole"),
+        ("ArcTanh[1.`20]", "Indeterminate", "an inexact pole"),
+        ("ArcTanh[-1.`20]", "Indeterminate", None),
+        ("ArcTanh[1.]", "Indeterminate", None),
+        ("ArcTanh[1.0000001`5]", "Indeterminate", "the error covers the pole"),
+        ("ArcTanh[0.9999999`5]", "Indeterminate", None),
+        ("ArcTanh[1.000001`6]", "Indeterminate", None),
+        ("Precision[ArcTanh[0.999999`6.1]]", "1.26163", "the error misses it"),
+        ("ArcTanh[1]", "Infinity", "an exact pole"),
+        ("ArcCoth[1.`20]", "Indeterminate", None),
+        ("ArcCoth[-1.]", "Indeterminate", None),
+        ("ArcCoth[1.0000001`5]", "Indeterminate", None),
+        ("Precision[ArcCoth[1.000001`6.1]]", "1.26163", None),
+        ("ArcCoth[-1]", "-Infinity", None),
         ("Accuracy[Sinh[1000.`2]]", "-434.993", "no digit left, no bound"),
     ],
 )