Repository: pearu/functional_algorithms

Issue:
Add complex arctan
Per https://github.com/pearu/complex_function_validation/tree/main/mpmath_jax_results_dev, JAX arctan accuracy is validated as BAD.

Discussion:

Patch:
diff --git a/README.md b/README.md
index 3be09a9..c223277 100644
--- a/README.md
+++ b/README.md
@@ -39,6 +39,8 @@ following math functions:
   al](https://dl.acm.org/doi/10.1145/275323.275324) algorithm for
   complex `asin`,
 - `asinh(z: complex | float)`, using the relation `asinh(z) = -I * asin(z * I)`,
+- `atan(z: complex | float)`, using the relation `atan(z) = -I * atanh(z * I)`,
+- `atanh(z: complex | float)`, using a custom algorithm,
 - `hypot(x: float, y: float)` and `absolute(z: complex)`,
 - `square(z: complex | float)`.
 
diff --git a/functional_algorithms/__init__.py b/functional_algorithms/__init__.py
index ae61989..1b6b924 100644
--- a/functional_algorithms/__init__.py
+++ b/functional_algorithms/__init__.py
@@ -1,4 +1,4 @@
-__all__ = ["Constext", "Expr", "targets"]
+__all__ = ["Constext", "Expr", "targets", "TextImage"]
 
 try:
     from ._version import __version__
@@ -13,4 +13,6 @@ except ImportError:
 
 from .expr import Expr
 from .context import Context
+from . import algorithms
 from . import targets
+from .textimage import TextImage
diff --git a/functional_algorithms/algorithms.py b/functional_algorithms/algorithms.py
index 4647038..d27be54 100644
--- a/functional_algorithms/algorithms.py
+++ b/functional_algorithms/algorithms.py
@@ -779,3 +779,177 @@ def log1p(ctx, z: complex | float):
     if z.is_complex:
         return complex_log1p(ctx, z)
     return ctx.log1p(z)
+
+
+def atanh_imag_is_half_pi(ctx, x: float):
+    """Return smallest positive x such that imag(atanh(I*x)) == pi/2
+
+    Using `largest` to detect the floating point type: float16,
+    float32, or float64.
+    """
+    largest = ctx.constant("largest", x)
+    # see tools/tune_atanh.py for computing the following constants:
+    fp64 = ctx.constant(5805358775541310.0, x)
+    fp32 = ctx.constant(62919776.0, x)
+    fp16 = ctx.constant(1028.0, x)
+    return ctx(ctx.select(largest > 1e308, fp64, ctx.select(largest > 1e38, fp32, fp16)))
+
+
+def inverse_fp_negeps(ctx, x: float):
+    """Return smallest positive x such that x + 1.0 == x
+
+    Using `largest` to detect the floating point type: float16,
+    float32, or float64.
+    """
+    import numpy
+
+    largest = ctx.constant("largest", x)
+
+    def get_value(dtype):
+        return float(numpy.nextafter(dtype(1 / numpy.finfo(dtype).epsneg), dtype(numpy.inf)))
+
+    fp64 = ctx.constant(get_value(numpy.float64), x)
+    fp32 = ctx.constant(get_value(numpy.float32), x)
+    fp16 = ctx.constant(get_value(numpy.float16), x)
+    return ctx(ctx.select(largest > 1e308, fp64, ctx.select(largest > 1e38, fp32, fp16)))
+
+
+def complex_atanh(ctx, z: complex):
+    """Inverse hyperbolic tangent on complex inputs:
+
+    atanh(z) = (log1p(z) - log1p(-z)) / 2
+
+    Algorithm derivation
+    --------------------
+
+    We have
+
+      log1p(x + I * y) = log((x + 1)**2 + y ** 2) / 2 + I * atan2(y, x + 1)
+
+    where
+
+      x and y are real and imaginary parts of the input to log1p, and
+      I is imaginary unit.
+
+    The real part
+    -------------
+
+    We have
+
+      real(atanh(x + I * y)) = (log((x + 1)**2 + y ** 2) - log((x - 1)**2 + y ** 2)) / 4
+        = log(((x + 1)**2 + y ** 2) / ((x - 1)**2 + y ** 2)) / 4
+        = log1p(((x + 1)**2 + y ** 2) / ((x - 1)**2 + y ** 2) - 1) / 4
+        = log1p(4 * x / ((x - 1)**2 + y ** 2)) / 4
+        = sign(x) * log1p(4 * abs(x) / ((abs(x) - 1)**2 + y ** 2)) / 4      Eq 2.3.
+
+    However, when abs(x) is large so that (~= denotes equality between
+    floating point numbers)
+
+      abs(x) + 1 ~= abs(x)
+
+    or equivalently,
+
+      abs(x) > 1 / epsneg ** 2
+
+    where
+
+      epsneg = 1 - nextafter(1, 0)
+
+    we'll find
+
+      real(atanh(x + I * y))
+        = sign(x) * log1p(4 * x / (x * y * (x / y + y / x))) / 4
+        = log1p(4 / y / (x / y + y / x)) / 4                                Eq 5.
+
+    If abs(y) is not large, so that
+
+      x / y + y / x ~= x / y
+
+    or equivalently,
+
+      (x / y) ** 2 + 1 ~= (x / y) ** 2
+      (x / y) ** 2 > 1 / epsneg ** 2
+      abs(y) < abs(x) * epsneg
+
+    then we'll use
+
+      real(atanh(x + I * y)) = log1p(4 / x) / 4                             Eq 4.
+
+    The imaginary part
+    ------------------
+
+      imag(atanh(x + I * y)) = (atan2(y, x + 1) - atan2(y, x - 1)) / 2
+        = atan2(y * (1/(x+1) - 1/(x-1)), 1 + y ** 2 / (x ** 2 - 1))
+        = atan2(-2 * y, x ** 2 + y ** 2 - 1)
+        = atan2(-2 * y, (x + 1) * (x - 1) + y * y)
+
+    For large values of abs(x) or abs(y), that is, when
+
+      x ** 2 + y ** 2 > 1 / epsneg ** 2
+
+    we have
+
+        imag(atanh(x + I * y)) = sign(y) * pi / 2
+    """
+    x = z.real
+    y = z.imag
+
+    zero = ctx.constant(0.0, x)
+    one = ctx.constant(1.0, x)
+    four = ctx.constant(4.0, x)
+    two = ctx.constant(2.0, x)
+    half = ctx.constant(0.5, x)
+    quarter = ctx.constant(0.25, x)
+    pi = ctx.constant("pi", x)
+    inv_negeps = inverse_fp_negeps(ctx, x)
+    safe_max = inv_negeps * inv_negeps
+
+    ax = abs(x)
+    ay = abs(y)
+    naxm1 = one - ax
+    y2 = y * y
+    sx = ctx.select(x >= 0, one, -one)
+    sy = ctx.select(y >= 0, one, -one)
+    in_safe_region = ctx.And(ax < safe_max, ay < safe_max)
+
+    # real part
+    arg23 = ax / (naxm1 * naxm1 + y2)  # Eq 2.3
+    arg4 = one / ax  # Eq 4
+    arg5 = one / (ax / y + y / ax) / y  # Eq 5.
+    arg5 = ctx.select(ctx.Or(x.is_inf, y.is_inf), zero, arg5)
+    arg = ctx.select(in_safe_region, arg23, ctx.select(ay * inv_negeps < ax, arg4, arg5))
+    real = sx * ctx.log1p(four * arg) * quarter
+
+    # imaginary part
+    imag = ctx.select(in_safe_region, ctx.atan2(y + y, naxm1 * (one + ax) - y2), sy * pi) * half
+
+    return ctx.complex(real, imag)
+
+
+def atanh(ctx, z: complex):
+    """Inverse hyperbolic tangent on real and complex inputs:
+
+    See complex_atanh for more information.
+    """
+    if not z.is_complex:
+        return ctx.atanh(z)
+    return complex_atanh(ctx, z)
+
+
+def complex_atan(ctx, z: complex):
+    """Arcus tangent on complex inputs:
+
+    atan(z) = -I * atanh(I * z)
+    """
+    w = complex_atanh(ctx, ctx.complex(-z.imag, z.real))
+    return ctx.complex(w.imag, -w.real)
+
+
+def atan(ctx, z: complex | float):
+    """Arcus tangent on complex and real inputs.
+
+    See complex_atan for more information.
+    """
+    if not z.is_complex:
+        return ctx.atan(z)
+    return complex_atan(ctx, z)
diff --git a/functional_algorithms/expr.py b/functional_algorithms/expr.py
index 29e0adf..bf711a2 100644
--- a/functional_algorithms/expr.py
+++ b/functional_algorithms/expr.py
@@ -1,7 +1,6 @@
 import struct
 import warnings
 from .utils import UNSPECIFIED
-from . import algorithms
 from .typesystem import Type
 from .rewrite import rewrite
 
@@ -105,7 +104,7 @@ def toidentifier(value):
             return "neg" + str(-value)
         return str(value)
     elif isinstance(value, float):
-        return "f" + hex(struct.unpack("<I", struct.pack("<f", value))[0])[1:]
+        return "f" + hex(struct.unpack("<Q", struct.pack("<d", value))[0])[1:]
     elif isinstance(value, complex):
         return "c" + toidentifier(value.real) + toidentifier(value.imag)
     elif isinstance(value, str):
@@ -194,7 +193,7 @@ class Expr:
             assert len(operands) == 2 and isinstance(operands[0], str) and isinstance(operands[1], Type), operands
         elif kind == "constant":
             assert len(operands) == 2
-            assert isinstance(operands[0], (int, float, bool, complex, str, Expr))
+            assert isinstance(operands[0], (int, float, bool, complex, str, Expr)), type(operands[0])
             if isinstance(operands[0], Expr):
                 assert operands[0].context is context.alt, operands
             if context.alt is not None and not isinstance(operands[0], Expr):
diff --git a/functional_algorithms/rewrite.py b/functional_algorithms/rewrite.py
index 696ce22..caea6b6 100644
--- a/functional_algorithms/rewrite.py
+++ b/functional_algorithms/rewrite.py
@@ -280,7 +280,6 @@ class Rewriter:
 
     def _compare(self, expr, relop):
         x, y = expr.operands
-
         if x.kind == "constant" and y.kind == "constant":
             xvalue, xlike = x.operands
             yvalue, ylike = y.operands
@@ -296,6 +295,12 @@ class Rewriter:
                 ctx = expr.context
                 return ctx.constant(r, ctx.symbol(None, "boolean"))
 
+            if isinstance(xvalue, str):
+                if xvalue == "largest":
+                    # cannot rewrite as largest depends on dtype
+                    if isinstance(yvalue, (int, float)):
+                        return
+
             self._todo(expr)
 
     def ge(self, expr):
diff --git a/functional_algorithms/targets/numpy.py b/functional_algorithms/targets/numpy.py
index fcaf5da..ce127c2 100644
--- a/functional_algorithms/targets/numpy.py
+++ b/functional_algorithms/targets/numpy.py
@@ -107,6 +107,7 @@ kind_to_target = dict(
     ge="({0}) >= ({1})",
     eq="numpy.equal({0}, {1}, dtype=numpy.bool_)",
     ne="({0}) != ({1})",
+    nextafter="numpy.nextafter({0}, {1})",
 )
 
 constant_to_target = dict(
diff --git a/functional_algorithms/targets/stablehlo.py b/functional_algorithms/targets/stablehlo.py
index 776b816..91b0400 100644
--- a/functional_algorithms/targets/stablehlo.py
+++ b/functional_algorithms/targets/stablehlo.py
@@ -98,6 +98,7 @@ kind_to_target = dict(
     ge=None,
     eq=None,
     ne=None,
+    nextafter="CHLO_NextAfterOp",
 )
 
 
diff --git a/functional_algorithms/textimage.py b/functional_algorithms/textimage.py
new file mode 100644
index 0000000..f05662e
--- /dev/null
+++ b/functional_algorithms/textimage.py
@@ -0,0 +1,193 @@
+import numpy
+from collections import defaultdict
+
+
+class TextImage:
+    """Image for text terminals."""
+
+    def __init__(self, rows=0, cols=0, bg=" "):
+        """
+        Parameters
+        ----------
+        rows: int
+          Initial number of image rows, corresponds to image height.
+        cols: int
+          Initial number of columns, corresponds to image width.
+        bg: str
+          Single character used as the background of the image.
+        """
+        self.bg = bg
+        self.table = numpy.full((rows, cols), bg, dtype="c")
+
+    def _allocate(self, row, col):
+        """Re-allocate table so that given table cell is accessible.
+
+        This corresponds to expanding image to right or down.
+
+        Notice that for an originally empty table, the shape of
+        resulting table will be (row + 1, col + 1).
+        """
+        if row >= self.table.shape[0]:
+            self.table = numpy.vstack(
+                (self.table, numpy.full((row - self.table.shape[0] + 1, self.table.shape[1]), self.bg, dtype=self.table.dtype))
+            )
+        if col >= self.table.shape[1]:
+            self.table = numpy.hstack(
+                (self.table, numpy.full((self.table.shape[0], col - self.table.shape[1] + 1), self.bg, dtype=self.table.dtype))
+            )
+
+    def _handle_negative_location(self, row, col):
+        while row < 0:
+            row += self.table.shape[0]
+        while col < 0:
+            col += self.table.shape[1]
+        return row, col
+
+    def _shift_by_loc(self, row, col, loc, shape):
+        drow = dcol = 0
+        if loc in {"ul", "l", "ll"}:
+            pass
+        elif loc in {"uc", "c", "lc"}:
+            dcol = -(shape[1] // 2)
+        elif loc in {"ur", "r", "lr"}:
+            dcol = -shape[1]
+        else:
+            assert 0, loc  # unreachable
+        if loc in {"ul", "uc", "ur"}:
+            pass
+        elif loc in {"l", "c", "r"}:
+            drow = -(shape[0] // 2)
+        elif loc in {"ll", "lc", "lr"}:
+            drow = -shape[0]
+        else:
+            assert 0, loc  # unreachable
+        return row + drow, col + dcol
+
+    @staticmethod
+    def _value_to_text(value):
+        if isinstance(value, numpy.ndarray):
+            dtype = value.dtype
+            assert value.shape == (), value.shape
+        elif isinstance(value, numpy.floating):
+            dtype = type(value)
+        else:
+            assert 0, type(value)  # not impl
+        f = numpy.finfo(dtype)
+        if numpy.isnan(value):
+            return "nan"
+        if numpy.isposinf(value):
+            return "+inf"
+        if numpy.isneginf(value):
+            return "-inf"
+        if value == 0:
+            return "0"
+        if value == f.max:
+            return "max"
+        if value == f.min:
+            return "min"
+        if value == f.tiny:
+            return "tiny"
+        if value == -f.tiny:
+            return "-tiny"
+        return f"{value:1.0e}".replace("e-0", "e-").replace("e+0", "e").replace("e+", "e").replace("e0", "")
+
+    @classmethod
+    def fromstring(cls, text):
+        image = cls()
+        row = 0
+        col = 0
+        for i, line in enumerate(text.splitlines()):
+            image._allocate(row + i, col + len(line) - 1)
+            image.table[row + i : row + i + 1, col : col + len(line)] = line
+        return image
+
+    @classmethod
+    def fromseq(cls, seq, colstart="", colsep=" ", colend="", mintextwidth=None, maxtextwidth=None):
+        image = cls()
+        col_widths = defaultdict(int)
+        data = []
+        for i, row in enumerate(seq):
+            d = []
+            for j, value in enumerate(row):
+                text = cls._value_to_text(value)
+                col_widths[j] = max(col_widths[j], len(text))
+                d.append(text)
+            data.append(d)
+
+        if maxtextwidth is None:
+            maxtextwidth = max(col_widths.values())
+        if mintextwidth is None:
+            mintextwidth = 0
+
+        offsets = [0]
+        for j, w in sorted(col_widths.items()):
+            offsets.append(
+                offsets[-1] + (len(colstart) if j > 0 else 0) + min(max(w, mintextwidth), maxtextwidth) + len(colsep)
+            )
+
+        image._allocate(len(data) - 1, offsets[-1] + len(colend) - len(colsep))
+
+        for i, row in enumerate(data):
+            for j, text in enumerate(row):
+                prefix = colstart if j == 0 else ""
+                suffix = colsep if j < len(row) - 1 else colend
+                t = text[:maxtextwidth]
+                if len(t) < mintextwidth:
+                    t += " " * (mintextwidth - len(t))
+                image.insert(i, offsets[j], prefix + t + suffix)
+
+        return image
+
+    def tostring(self):
+        lst = []
+        for row in self.table:
+            lst.append((b"".join(row)).decode().rstrip())
+        return "\n".join(lst)
+
+    def __str__(self):
+        return self.tostring()
+
+    def append(self, row, col, *args, **kwargs):
+        row, col = self._handle_negative_location(row, col)
+        self.insert(row + 1, col, *args, **kwargs)
+
+    def insert(self, row, col, subimage, loc="ul"):
+        """Insert subimage to image.
+
+        Insertion will override existing cells of the current text image.
+
+        Parameters
+        ----------
+        row, col: (int, int)
+          The coordinates of loc within the subimage.
+        subimage: {TextImage, ndarray, str}
+          Subimage.
+        loc: {'ul', 'uc', 'ur', 'l', 'c', 'r', 'll', 'lc', 'lr'}
+          The location within subimage.
+        """
+        if isinstance(subimage, str):
+            subimage = TextImage.fromstring(subimage)
+
+        if isinstance(subimage, TextImage):
+            subimage = subimage.table
+
+        assert isinstance(subimage, numpy.ndarray), type(subimage)
+
+        self.get_view(row, col, subimage.shape, loc=loc)[:] = subimage
+
+    def get_view(self, row, col, shape, loc="ul"):
+        """Return subimage."""
+        row, col = self._handle_negative_location(row, col)
+        row, col = self._shift_by_loc(row, col, loc, shape)
+        self._allocate(row + shape[0] - 1, col + shape[1] - 1)
+        return self.table[row : row + shape[0], col : col + shape[1]]
+
+    def fill(self, row, col, mask, loc="ul", symbol="x"):
+        """Fill subimage with symbol per mask."""
+        assert mask.dtype == numpy.bool_, mask.dtype
+        textmask = self.get_view(row, col, mask.shape, loc=loc)
+        textmask[numpy.where(mask)] = symbol
+
+        # textmask = numpy.full(mask.shape, self.bg, dtype='c')
+        # textmask[numpy.where(mask)] = symbol
+        # self.insert(row, col, textmask)
diff --git a/functional_algorithms/utils.py b/functional_algorithms/utils.py
index d6872ff..879ca2b 100644
--- a/functional_algorithms/utils.py
+++ b/functional_algorithms/utils.py
@@ -171,6 +171,9 @@ class vectorize_with_mpmath(numpy.vectorize):
                         return -ctx.zero if sign else ctx.zero
                     if exp + bc > self.float_maxexp[fp_format]:
                         return ctx.ninf if sign else ctx.inf
+                    # conversion from mpz to longdouble is buggy, see
+                    # aleaxit/gmpy#507, hence we convert mpz to int:
+                    man = int(man)
                     man = dtype(-man if sign else man)
                     r = numpy.ldexp(man, exp)
                     assert numpy.isfinite(r), (x, r, x._mpf_, man)
@@ -544,6 +547,48 @@ class mpmath_array_api:
             return self.arctan2(x.imag, x.real)
         return ctx.nan
 
+    def arctanh(self, x):
+        ctx = x.context
+        if isinstance(x, ctx.mpc):
+            pi = ctx.pi
+            inf = ctx.inf
+            zero = ctx.zero
+            if ctx.isinf(x.real) or ctx.isinf(x.imag):
+                if x.imag < 0:
+                    return ctx.make_mpc((zero._mpf_, (-pi / 2)._mpf_))
+                return ctx.make_mpc((zero._mpf_, (pi / 2)._mpf_))
+        r = ctx.atanh(x)
+        if isinstance(x, ctx.mpc):
+            # The following if-block ensures compatibiliy with
+            # numpy.arctanh.
+            if x.imag == 0 and x.real > 1:
+                return r.conjugate()
+        if isinstance(x, ctx.mpf) and isinstance(r, ctx.mpc):
+            # otherwise, mpmath.atanh would return complex value
+            return ctx.nan
+        return r
+
+    def arctan(self, x):
+        ctx = x.context
+        if isinstance(x, ctx.mpc):
+            pi = ctx.pi
+            inf = ctx.inf
+            zero = ctx.zero
+            if ctx.isinf(x.real) or ctx.isinf(x.imag):
+                if x.real < 0:
+                    return ctx.make_mpc(((-pi / 2)._mpf_, zero._mpf_))
+                return ctx.make_mpc(((pi / 2)._mpf_, zero._mpf_))
+        r = ctx.atan(x)
+        if isinstance(x, ctx.mpc):
+            # The following if-block ensures compatibiliy with
+            # numpy.arctan:
+            if x.real == 0 and x.imag < -1:
+                return ctx.make_mpc(((-r.real)._mpf_, r.imag._mpf_))
+        if isinstance(x, ctx.mpf) and isinstance(r, ctx.mpc):
+            # otherwise, mpmath.atan would return complex value
+            return ctx.nan
+        return r
+
 
 class numpy_with_mpmath:
     """Namespace of universal functions on numpy arrays that use mpmath
@@ -556,7 +601,9 @@ class numpy_with_mpmath:
         self.params = params
 
     def __getattr__(self, name):
-        name = dict(asinh="arcsinh", acos="arccos", asin="arcsin", acosh="arccosh").get(name, name)
+        name = dict(asinh="arcsinh", acos="arccos", asin="arcsin", acosh="arccosh", atan="arctan", atanh="arctanh").get(
+            name, name
+        )
         if name in self._vfunc_cache:
             return self._vfunc_cache[name]
         if hasattr(mpmath_array_api, name):
@@ -681,7 +728,15 @@ def extra_samples(name, dtype):
                 values.append(complex(numpy.nextafter(x, fdtype(1)), y))
 
             values.append(-0.24246988 - 0.49096385j)
-
+        if name == "atanh":
+            values.extend(
+                [
+                    -0.9939759 - 0.00042062782j,
+                    0.9939759 - 0.00085302145j,
+                    -0.73795176 + 0.0877259j,
+                    315387600000000 + 36880000000000j,
+                ]
+            )
     return numpy.array(values, dtype=dtype)
 
 
@@ -694,6 +749,8 @@ def real_samples(
     include_nan=False,
     include_huge=True,
     nonnegative=False,
+    min_value=None,
+    max_value=None,
 ):
     """Return a 1-D array of real line samples.
 
@@ -717,17 +774,49 @@ def real_samples(
       maximal value.
     nonnegative: bool
       When True, finite samples are all non-negative.
+    min_value, max_value: dtype, None
+      When min_value or max_value is specified, use for constructing
+      uniform samples. Parameters include_infinity, include_nan,
+      include_huge, and nonnegative are silently ignored.
     """
     if isinstance(dtype, str):
         dtype = getattr(numpy, dtype)
     assert dtype in {numpy.float32, numpy.float64}, dtype
     utype = {numpy.float32: numpy.uint32, numpy.float64: numpy.uint64}[dtype]
     fi = numpy.finfo(dtype)
-    num = size // 2 if not nonnegative else size
+    user_specified_bounds = min_value is not None or max_value is not None
+    num = size // 2 if not nonnegative and not user_specified_bounds else size
     if include_infinity:
         num -= 1
-    min_value = dtype(fi.smallest_subnormal if include_subnormal else fi.smallest_normal)
-    max_value = dtype(fi.max)
+
+    min_pos_value = dtype(fi.smallest_subnormal if include_subnormal else fi.smallest_normal)
+    if min_value is None:
+        min_value = min_pos_value
+    else:
+        min_value = dtype(min_value)
+
+    if max_value is None:
+        max_value = dtype(fi.max)
+    else:
+        max_value = dtype(max_value)
+
+    if min_value >= max_value:
+        raise ValueError(f"minimal value (={min_value}) cannot be greater than maximal value (={max_value})")
+
+    if user_specified_bounds:
+        if min_value > 0:
+            return numpy.linspace(min_value.view(utype), max_value.view(utype), num=num, dtype=utype).view(dtype)
+        if max_value < 0:
+            return -numpy.linspace(abs(max_value).view(utype), abs(min_value).view(utype), num=num, dtype=utype).view(dtype)
+        neg_diff = diff_ulp(abs(min_value), min_pos_value)
+        pos_diff = diff_ulp(abs(max_value), min_pos_value)
+        neg_num = int(neg_diff * num / (neg_diff + pos_diff))
+        pos_num = num - neg_num
+        neg_part = -numpy.linspace(abs(min_value).view(utype), min_pos_value.view(utype), num=neg_num, dtype=utype).view(dtype)
+        pos_part = numpy.linspace(min_pos_value.view(utype), max_value.view(utype), num=pos_num, dtype=utype).view(dtype)
+        if include_zero:
+            return numpy.unique(numpy.concatenate([neg_part, numpy.array([0], dtype=dtype), pos_part]))
+        return numpy.unique(numpy.concatenate([neg_part, pos_part]))
     if 1:
         # The following method gives a sample distibution that is
         # uniform with respect to ULP distance between positive
@@ -778,6 +867,10 @@ def complex_samples(
     include_nan=False,
     include_huge=True,
     nonnegative=False,
+    min_real_value=None,
+    max_real_value=None,
+    min_imag_value=None,
+    max_imag_value=None,
 ):
     """Return a 2-D array of complex plane samples.
 
@@ -792,6 +885,7 @@ def complex_samples(
     include_subnormal: bool
     include_nan: bool
     nonnegative: bool
+    min_real_value, max_real_value, min_imag_value, max_imag_value: dtype
     """
     if isinstance(dtype, str):
         dtype = getattr(numpy, dtype)
@@ -805,6 +899,8 @@ def complex_samples(
         include_nan=include_nan,
         include_huge=include_huge,
         nonnegative=nonnegative,
+        min_value=min_real_value,
+        max_value=max_real_value,
     )
     im = real_samples(
         size[1],
@@ -815,6 +911,8 @@ def complex_samples(
         include_nan=include_nan,
         include_huge=include_huge,
         nonnegative=nonnegative,
+        min_value=min_imag_value,
+        max_value=max_imag_value,
     )
     complex_dtype = {numpy.float32: numpy.complex64, numpy.float64: numpy.complex128}[dtype]
     real_part = re.reshape((-1, re.size)).repeat(im.size, 0).astype(complex_dtype)
@@ -1192,7 +1290,8 @@ def function_validation_parameters(func_name, dtype):
     elif func_name == "log1p":
         extra_prec_multiplier = 10  # remove when mpmath#803 becomes available
         max_bound_ulp_width = dict(complex64=3, complex128=3).get(dtype_name, max_bound_ulp_width)
-
+    elif func_name in {"atanh", "atan"}:
+        extra_prec_multiplier = 20
     return dict(
         extra_prec_multiplier=extra_prec_multiplier,
         max_valid_ulp_count=max_valid_ulp_count,
diff --git a/results/README.md b/results/README.md
index 55d1008..9df4e3a 100644
--- a/results/README.md
+++ b/results/README.md
@@ -56,6 +56,14 @@ in complex plane where `ulp_width` needs to be larger to pass the
 | asinh | float64 | 825453 | 174482 | 66 | - | - |
 | asinh | complex64 | 807415 | 193174 | 1320 | 92 | - |
 | asinh | complex128 | 687179 | 313978 | 844 | - | - |
+| atan | float32 | 957469 | 42532 | - | - | - |
+| atan | float64 | 992077 | 7924 | - | - | - |
+| atan | complex64 | 850741 | 146392 | 4784 | 84 | - |
+| atan | complex128 | 936337 | 65264 | 400 | - | - |
+| atanh | float32 | 986261 | 13740 | - | - | - |
+| atanh | float64 | 995895 | 4106 | - | - | - |
+| atanh | complex64 | 850741 | 146392 | 4784 | 84 | - |
+| atanh | complex128 | 936337 | 65264 | 400 | - | - |
 | square | float32 | 997347 | 2654 | - | - | - |
 | square | float64 | 999593 | 408 | - | - | - |
 | square | complex64 | 976809 | 25192 | - | - | - |
diff --git a/results/estimate_accuracy.py b/results/estimate_accuracy.py
index 754540d..78d6227 100644
--- a/results/estimate_accuracy.py
+++ b/results/estimate_accuracy.py
@@ -31,6 +31,14 @@ def get_inputs():
         ("asinh", np.float64, {}),
         ("asinh", np.complex64, {}),
         ("asinh", np.complex128, {}),
+        ("atan", np.float32, {}),
+        ("atan", np.float64, {}),
+        ("atan", np.complex64, {}),
+        ("atan", np.complex128, {}),
+        ("atanh", np.float32, {}),
+        ("atanh", np.float64, {}),
+        ("atanh", np.complex64, {}),
+        ("atanh", np.complex128, {}),
         ("square", np.float32, {}),
         ("square", np.float64, {}),
         ("square", np.complex64, {}),
diff --git a/tools/tune_atanh.py b/tools/tune_atanh.py
new file mode 100644
index 0000000..f5695b7
--- /dev/null
+++ b/tools/tune_atanh.py
@@ -0,0 +1,520 @@
+"""Tools for tuning atanh algorithm.
+"""
+
+import functional_algorithms as fa
+import numpy as np
+
+
+def unarray(func):
+
+    def wrapper(*args):
+        return func(*args)[()]
+
+    return wrapper
+
+
+def get_reference_func(name, dtype):
+    params = fa.utils.function_validation_parameters(name, dtype)
+    extra_prec_multiplier = params["extra_prec_multiplier"]
+    return getattr(fa.utils.numpy_with_mpmath(extra_prec_multiplier=extra_prec_multiplier), name)
+
+
+def bisect(pred, a, b, dtype=None, extra_args=()):
+    """find x such that
+
+      pred(prev(x)) == pred(x) and pred(x) != pred(next(x))
+
+    or
+
+     pred(prev(x)) != pred(x) and pred(x) == pred(next(x))
+
+    within interval [a, b].
+
+    pred(a) != pred(b) must hold
+    """
+    if dtype is None:
+        dtype = type(a)
+    if a > b:
+        a, b = b, a
+    pa = pred(a, *extra_args)
+    pb = pred(b, *extra_args)
+    assert pa != pb, (a, b, pa, pb)
+    while True:
+        m = dtype((a + b) / 2)
+        if m == a or m == b:
+            return m
+        pm = pred(m, *extra_args)
+        if pa != pm:
+            b = m
+            pb = pm
+        elif pb != pm:
+            a = m
+            pa = pm
+        else:
+            assert 0, (a, b, m)
+
+
+def maximal_bisect(pred, a, b, mn, mx, dtype):
+    """find x such that
+
+       bisect(pred, mn, mx, extra_args=(x,))
+
+    result has a maximal possible value whereas a < x and x < b.
+
+    The maximal value must exist in range [mn, mx].
+    """
+    a, b = dtype(a), dtype(b)
+    fa = bisect(pred, mn, mx, dtype, extra_args=(a,))
+    fb = bisect(pred, mn, mx, dtype, extra_args=(b,))
+
+    while True:
+        m1 = a + (b - a) / 3
+        m2 = dtype((m1 + b) / 2)
+        m1 = dtype(m1)
+        f1 = bisect(pred, mn, mx, dtype, extra_args=(m1,))
+        if m1 == m2:
+            return f1
+        f2 = bisect(pred, mn, mx, dtype, extra_args=(m2,))
+        if f1 > f2:
+            a, b = a, m2
+            fa, fb = fa, f2
+        else:
+            a, b = m1, b
+            fa, fb = f1, fb
+        mn = min(fa, fb)
+
+
+def main_imag():
+    """The imaginary part of complex atanh(z).
+
+    Observe the mask of fp(abs(imag(atanh(z)))) == fp(pi/2) on the
+    first quarted of complex plane (0<=z.real<inf, 0<=z.imag<inf):
+
+    ^  log(imag)
+    |
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+    xxxxxxxxxxxxxxxxxxxZxxxxxxxxxxxx
+                       xxxxxxxxxxxxx
+                      xxxxxxxxxxxxxx
+                      xxxxxxxxxxxxxx
+                     xxxxxxxxxxxxxxx
+                     xxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+                    xxxxxxxxxxxxxxxx
+    o               xxxxxxxxxxxxxxxx  -> log(real)
+
+    This shows that there exists parameters X and Y (or Z=X+I*Y as
+    marked in the mask representation above) such that
+
+      atanh(z) == pi / 2
+
+    for any z such that z.imag > Y or if z.real > X where the equality
+    `==` is defined as the equality of floating point numbers.
+
+    In the following we'll find the parameters X and Y for different
+    floating point systems:
+
+      float16: X, Y=(514.0, 1028.0)
+      float32: X, Y=(31459888.0, 62919776.0)
+      float64: X, Y=(2902679387770655.0, 5805358775541310.0)
+      longdouble: X, Y=(1.716601425640972592e+19, 3.433202851281945184e+19)
+
+    Notice the following relations:
+
+      Y is approximately 2 * X
+
+      1 / eps < Y < 8 / eps
+
+    for all fp types where
+
+      eps = nextafter(fp(1), fp(2)) - fp(1)
+    """
+
+    dtype = np.complex64
+    fdtype = np.float32
+
+    params = fa.utils.function_validation_parameters("atanh", dtype)
+    max_valid_ulp_count = params["max_valid_ulp_count"]
+    max_bound_ulp_width = params["max_bound_ulp_width"]
+    extra_prec_multiplier = params["extra_prec_multiplier"]
+    func = getattr(
+        fa.utils.numpy_with_mpmath(extra_prec_multiplier=extra_prec_multiplier, flush_subnormals=True),
+        "atanh",
+    )
+
+    size = 31
+    samples = fa.utils.complex_samples(
+        (size, size), dtype=dtype, include_huge=False, include_subnormal=False, nonnegative=True
+    )[::-1]
+
+    values = func.call(samples)
+
+    mask = np.float32(abs(values.imag)) == np.float32(np.pi / 2)
+
+    timage = fa.TextImage()
+    timage.fill(0, 0, mask, symbol="x")
+    print(timage)
+
+    # compute X and Y for
+    for fdtype in [np.float16, np.float32, np.float64, np.float128]:
+        dtype = {np.float16: np.complex64, np.float32: np.complex64, np.float64: np.complex128, np.float128: np.complex256}[
+            fdtype
+        ]
+
+        largest = np.finfo(fdtype).max / 8
+        target = fdtype("1.5707963267948966192313216916397514420985846996875529104874722")  # pi / 2
+
+        def f2(x, y):
+            z = np.empty((), dtype=dtype)
+            z.real = x
+            z.imag = y
+            r = func(z).imag
+            imag = fdtype(r)
+            return imag == target
+
+        def f1(y):
+            return f2(0, y)
+
+        Y = bisect(f1, fdtype(100), largest, fdtype)
+        X = maximal_bisect(f2, 0.1, Y * 0.9, 1, Y, fdtype)
+
+        print(f"{fdtype.__name__}: {X, Y=}", "8 / eps = %e" % (8 / np.finfo(fdtype).eps))
+
+
+def functional_nextafter(func, x, d):
+    """Return x1 such that
+
+    func(nextafter(x1, x)) == func(x)
+    func(x1) != func(x)
+    x1 > x if right else x1 < x
+    """
+    dtype = type(x)
+    f = func(x)
+    x1 = x
+    oo = dtype(d)
+    while func(x1) == f:
+        x1 = np.nextafter(x1, oo)
+    assert func(np.nextafter(x1, x)) == f
+    assert func(x1) != f
+    return x1
+
+
+def main_real():
+    """The real part of complex atanh(z)."""
+
+    dtype = np.complex128
+    fdtype = np.float64
+
+    params = fa.utils.function_validation_parameters("atanh", dtype)
+    max_valid_ulp_count = params["max_valid_ulp_count"]
+    max_bound_ulp_width = params["max_bound_ulp_width"]
+    extra_prec_multiplier = params["extra_prec_multiplier"]
+
+    atanh = get_reference_func("atanh", dtype)
+    log1p = get_reference_func("log1p", dtype)
+    log = get_reference_func("log", dtype)
+
+    def atanh1(z):
+        """Naive approach. The real part of the single precision float atanh
+              on the first quarter of complex plane compared to reference
+              atanh(z).real. Legend: (=) - exact match, (1) - match with 1
+              ULP difference, (2) - match with 2 ULPs difference, (c) -
+              match with less that 100 ULPs difference, ( ) - (no)match with
+              more than ULPs difference.
+
+        +inf
+        max
+        8e35
+        2e33
+        4e30
+        1e28
+        2e25
+        5e22
+        1e20
+        3e17    =============================
+        7e14    =========================
+        2e12    =====================
+        4e9     ================
+        8e6     ============
+        2e4     ========
+        4e1     ====                          cccc
+        1e-1    ==                           c2=cc
+        2e-4    ==                           c21cc
+        5e-7    =1                           c1=cc
+        1e-9    =1                           c1=cc
+        3e-12   =1                           c1=cc
+        6e-15   =1                           c1=cc
+        1e-17   =1                           c1=cc
+        4e-20   =1                           c1=cc
+        8e-23   =1                           c1=cc
+        2e-25   =1                           c1=cc
+        4e-28   =1                           c1=cc
+        1e-30   =1                           c1=cc
+        2e-33   =1                           c1=cc
+        5e-36   =1                           c1=cc
+        tiny    =1                           c1=cc
+        0       =1                           c1=cc
+                0     3e-32 1e-24 5e-17 3e-9  1e-1  5e6   2e14  9e21  4e29  2e37
+
+        """
+        x, y = z.real, z.imag
+        fptype = x.dtype.type if isinstance(x, np.ndarray) else type(x)
+
+        one = fptype(1)
+        quarter = fptype(0.25)
+        return log(((x + one) * (x + one) + y * y) / ((x - one) * (x - one) + y * y)) * quarter
+
+    def atanh2(z):
+        """
+        +inf    ==============================================================
+        max     ==============================================================
+        8e35    ==========================================================   =
+        2e33    ======================================================       =
+        4e30    ==================================================           =
+        1e28    =============================================                =
+        2e25    =========================================                    =
+        5e22    =====================================                        =
+        1e20    =================================                            =
+        3e17    =====================================1=======1               =
+        7e14    ===========================1========1======11=               =
+        2e12    ====================1==1====1====1==========1=               =
+        4e9     ===================1===11=1=1====111====1=111=               =
+        8e6     ============1==111=111=111==111=11==111===111=               =
+        2e4     =========1========111===11=1=1=====11=====111=               =
+        4e1     ========1==1=1==1=11==1=1==11=======1=====111=               =
+        1e-1    =11111===1==1111111=111===1======1==1=====111=               =
+        2e-4    =111111111=111=11111111=111==1===1==1=====111=               =
+        5e-7    =========================11=111==1==1=====111=               =
+        1e-9    =========================11=111==1==1=====111=               =
+        3e-12   =========================11=111==1==1=====111=               =
+        6e-15   =========================11=111==1==1=====111=               =
+        1e-17   =========================11=111==1==1=====111=               =
+        4e-20   =========================11=111==1==1=====111=               =
+        8e-23   =========================11=111==1==1=====111=               =
+        2e-25   =========================11=111==1==1=====111=               =
+        4e-28   =========================11=111==1==1=====111=               =
+        1e-30   =========================11=111==1==1=====111=               =
+        2e-33   =========================11=111==1==1=====111=               =
+        5e-36   =========================11=111==1==1=====111=               =
+        tiny    =========================11=111==1==1=====111=               =
+        0       =========================11=111==1==1=====111=               =
+                0     3e-32 1e-24 5e-17 3e-9  1e-1  5e6   2e14  9e21  4e29  2e37
+        """
+        x, y = z.real, z.imag
+        fptype = x.dtype.type if isinstance(x, np.ndarray) else type(x)
+
+        one = fptype(1)
+        four = fptype(4)
+        quarter = fptype(0.25)
+        return log1p(x / ((x - one) * (x - one) + y * y) * four) * quarter
+
+    def atanh3(z):
+        x, y = z.real, z.imag
+        fptype = x.dtype.type if isinstance(x, np.ndarray) else type(x)
+
+        one = fptype(1)
+        four = fptype(4)
+        quarter = fptype(0.25)
+        return -log1p(-x / ((x + one) * (x + one) + y * y) * four) * quarter
+
+    def atanh4(z):
+        """
+        +inf                                                                 ==
+        max                                                                  ==
+        8e35                                                                 ==
+        2e33                                                               c===
+        4e30                                                             c=====
+        1e28                                                           c=======
+        2e25                                                         c=========
+        5e22                                                       c===========
+        1e20                                                     c=============
+        3e17                                                   c===============
+        7e14                                                 c=================
+        2e12                                               1===================
+        4e9                                              1=====================
+        8e6                                            1=======================
+        2e4                                          1=========================
+        4e1                                        cc==========================
+        1e-1                                       cc==========================
+        2e-4                                       cc==========================
+        5e-7                                       cc==========================
+        1e-9                                       cc==========================
+        3e-12                                      cc==========================
+        6e-15                                      cc==========================
+        1e-17                                      cc==========================
+        4e-20                                      cc==========================
+        8e-23                                      cc==========================
+        2e-25                                      cc==========================
+        4e-28                                      cc==========================
+        1e-30                                      cc==========================
+        2e-33                                      cc==========================
+        5e-36                                      cc==========================
+        tiny                                       cc==========================
+        0                                          cc==========================
+                0     3e-32 1e-24 5e-17 3e-9  1e-1  5e6   2e14  9e21  4e29  2e37
+        """
+        x, y = z.real, z.imag
+        fptype = x.dtype.type if isinstance(x, np.ndarray) else type(x)
+
+        one = fptype(1)
+        four = fptype(4)
+        quarter = fptype(0.25)
+        return log1p(four / x) * quarter
+
+    def atanh5(z):
+        """
+        +inf    ==============================================================
+        max     ===============================================================
+        8e35    ==========================================================1=1==
+        2e33    =====================================================11========
+        4e30    =================================================1111=====1====
+        1e28    =============================================1=11=====1=1111===
+        2e25    ===========================================1==1===1==1==11=====
+        5e22    =================================================1===111111====
+        1e20    ====================================1=11==1=11==1==1=1==11=====
+        3e17    ==============================1=1111=====1===1==11===1=1=1==1==
+        7e14    ==========================1111111===1=1==21===1===1=====111====
+        2e12    ====================1====1===1===1==1==111===11=====1=1==11====
+        4e9     ===================11===1====1==1=11===11==1=1===1===1=1=1==1==
+        8e6     ===============1=1=1====1===11===11=21=========1==11=1=1=1=====
+        2e4     ========1===1==11==11===1=====1=2c cc11===11=1=====11111==1====
+        4e1     ====                               cc==1=1=1==1===1=11==1==1===
+        1e-1    =                                  cc====1=1==1===111=111=1====
+        2e-4    =                                  cc=1=1==1=====11===1===1  ==
+        5e-7    =                                  cc===1=1==1=1=1=1=11=1    ==
+        1e-9    =                                  cc11=1=1===1==1=1=1       ==
+        3e-12   =                                  cc11===11======1=         ==
+        6e-15   =                                  cc===1======1=1           ==
+        1e-17   =                                  cc===========             ==
+        4e-20   =                                  cc11======1               ==
+        8e-23   =                                  cc==1111=                 ==
+        2e-25   =                                  cc===11                   ==
+        4e-28   =                                  cc1=1                     ==
+        1e-30   =                                  cc1                       ==
+        2e-33   =                                  c                         ==
+        5e-36   =                                                            ==
+        tiny    =                                                            ==
+        0
+                0     3e-32 1e-24 5e-17 3e-9  1e-1  5e6   2e14  9e21  4e29  2e37
+        """
+        x, y = z.real, z.imag
+        fptype = x.dtype.type if isinstance(x, np.ndarray) else type(x)
+
+        four = fptype(4)
+        quarter = fptype(0.25)
+        return log1p(four / (x / y + y / x) / y) * quarter
+
+    def atanh6(z):
+        """
+        +inf     =============================================================
+        max      =============================================================
+        8e35     =========================================================   =
+        2e33     =====================================================       =
+        4e30     =================================================           =
+        1e28     ============================================                =
+        2e25     ========================================                    =
+        5e22     ====================================                        =
+        1e20     ================================                            =
+        3e17     ============================                                =
+        7e14     ========================                                    =
+        2e12     ====================                                        =
+        4e9      ===============                                             =
+        8e6      ===========                                                 =
+        2e4      =======                                                     =
+        4e1      ===                                                         =
+        1e-1     =                                                           =
+        2e-4     =                                                           =
+        5e-7     1                                                           =
+        1e-9     1                                                           =
+        3e-12    1                                                           =
+        6e-15    1                                                           =
+        1e-17    1                                                           =
+        4e-20    1                                                           =
+        8e-23    1                                                           =
+        2e-25    1                                                           =
+        4e-28    1                                                           =
+        1e-30    1                                                           =
+        2e-33    1                                                           =
+        5e-36    1                                                           =
+        tiny     1                                                           =
+        0        1                                                           =
+                0     3e-32 1e-24 5e-17 3e-9  1e-1  5e6   2e14  9e21  4e29  2e37
+        """
+        x, y = z.real, z.imag
+        fptype = x.dtype.type if isinstance(x, np.ndarray) else type(x)
+        quarter = fptype(0.25)
+        return log(x / x) * quarter  # it is effectively zero
+
+    size = 51
+
+    atanh_approx = atanh4
+
+    if 1:
+        min_real_value = None
+        max_real_value = None
+        min_imag_value = None
+        max_imag_value = None
+    else:
+        # study the surrounding of the given point:
+        z, dz = -1 - 0.00042062782j, 0.5 + 1e3j
+        z, dz = 315387600000000 + 36880000000000j, 1e8 + 1e8j
+        min_real_value = z.real - dz.real
+        max_real_value = z.real + dz.real
+        min_imag_value = z.imag - dz.imag
+        max_imag_value = z.imag + dz.imag
+
+    samples = fa.utils.complex_samples(
+        (2 * size, size),
+        dtype=dtype,
+        include_huge=False,
+        include_subnormal=False,
+        nonnegative=True,
+        min_real_value=min_real_value,
+        max_real_value=max_real_value,
+        min_imag_value=min_imag_value,
+        max_imag_value=max_imag_value,
+    )[::-1]
+
+    ref = atanh.call(samples).real
+
+    diff = fa.utils.diff_ulp(ref, atanh_approx(samples))
+
+    timage = fa.TextImage()
+    timage.fill(0, 10, diff == 0, symbol="=")
+    timage.fill(0, 10, diff == 1, symbol="1")
+    timage.fill(0, 10, diff == 2, symbol="2")
+    timage.fill(0, 10, (diff > 2) & (diff < 100), symbol="c")
+    real_axis = samples[0:1].real
+    imag_axis = samples[:, 0:1].imag
+    timage.insert(0, 2, fa.TextImage.fromseq(imag_axis))
+    timage.append(-1, 10, fa.TextImage.fromseq(real_axis[:, ::6], mintextwidth=5, maxtextwidth=5))
+    print(timage)
+
+    print(np.finfo(fdtype))
+
+
+if __name__ == "__main__":
+    # main_imag()
+    main_real()


Test patch:
diff --git a/functional_algorithms/tests/test_algorithms.py b/functional_algorithms/tests/test_algorithms.py
index 183318d..453723a 100644
--- a/functional_algorithms/tests/test_algorithms.py
+++ b/functional_algorithms/tests/test_algorithms.py
@@ -16,12 +16,18 @@ def dtype_name(request):
     return request.param
 
 
-@pytest.fixture(scope="function", params=["absolute", "acos", "acosh", "asin", "asinh", "hypot", "square", "sqrt", "angle"])
+@pytest.fixture(
+    scope="function",
+    params=["absolute", "acos", "acosh", "asin", "asinh", "hypot", "square", "sqrt", "angle", "atan", "atanh"],
+)
 def func_name(request):
     return request.param
 
 
-@pytest.fixture(scope="function", params=["absolute", "acos", "acosh", "asin", "asinh", "square", "sqrt", "angle", "log1p"])
+@pytest.fixture(
+    scope="function",
+    params=["absolute", "acos", "acosh", "asin", "asinh", "square", "sqrt", "angle", "log1p", "atan", "atanh"],
+)
 def unary_func_name(request):
     return request.param