Repository: pearu/functional_algorithms

Issue:
JAX acosh accuracy and correctness
number of samples: 2000400

```
  -inf    11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  -3e35   11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  -3e32   11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  -3e29   11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  -3e26   11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  -3e23   11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  -3e20   11111112111111111111111111111111111111111111111111111111111111111111111111111111111111111111211111111
  -3e17   11111112122111111111111111111111111111111111111111111111111111111111111111111111111111111121211111111
  -3e14   11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  -3e11   11111111111111==11111111111111111111111111111111111=========================111=111111=11111111111111
  -2e8    11111111111111====111111111111111111111111111111111=======================1111112111===11111111111111
  -2e5    11111111111111======11111111111111111111111111111111111111111111111111111111122221=====11111111111111
  -2e2    11111111111111======111112231111111111111111111111111111111111111111111112232311=======11111111111111
  -2e-1   11111111111111======11111222222111111111111111111111111111111111111111222222221========11111111111111
  -2e-4   11111111111111======11111222221111111111111111111111111111111111111111122222221========11111111111111
  -2e-7   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -2e-10  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -2e-13  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -2e-16  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -2e-19  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -2e-22  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -2e-25  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -2e-28  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -2e-31  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -1e-34  11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  -1e-37  EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE2==========11111111111111
  1e-37   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  1e-34   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  2e-31   11111111111111======1111111111111111111111111111111===================1111122211=======11111111111111
  2e-28   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  2e-25   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  2e-22   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  2e-19   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  2e-16   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  2e-13   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  2e-10   11111111111111======1111111111111111111111111111111===================111112221========11111111111111
  2e-7    11111111111111======1111111111111111111111111111111===================111122221========11111111111111
  2e-4    11111111111111======11111222222111111111111111111111111111111111111111222222221========11111111111111
  2e-1    11111111111111======111112222221111111111111111111111111111111111111112223333211=======11111111111111
  3e2     11111111111111======1111111111111111111111111111111=====================1111122221=====11111111111111
  3e5     11111111111111====111111111111111111111111111111111=======================1111122111===11111111111111
  3e8     11111111111111==11111111111111111111111111111111111=========================1111111111=11111111111111
  3e11    11111111112111111111111111111111111111111111111111111111111111111111111111111111111111111221211111111
  3e14    11111111112111111111111111111111111111111111111111111111111111111111111111111111111111111121211111111
  4e17    11111112111111111111111111111111111111111111111111111111111111111111111111111111111111111111211111111
  4e20    11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  4e23    11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  4e26    11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  4e29    11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  5e32    11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
  5e35    11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
          -inf  -3e29 -2e20 -1e11 -1e2  -8e-8 -6e-1 -4e-2 -3e-3 1e-33 1e-24 2e-15 3e-6  4e3   5e12  8e21  1e31
```

| ULP-difference |                             z |          jax:acosh(z) |       mpmath:acosh(z) |            numpy:acosh(z) |               fa:acosh(z) |
| -------------- | ----------------------------- | --------------------- | --------------------- | ------------------------- | ------------------------- |
|     2157060022 | (-1.9861786e+20-5.88902e-39j) | (47.43106+3.1415927j) | (47.43106-3.1415927j) |     (47.43106-3.1415927j) |     (47.43106-3.1415927j) |
|     2141945163 |     (-0.0989222-5.88902e-39j) |            1.6698805j |           -1.6698806j | (5.918046e-39-1.6698806j) | (5.918046e-39-1.6698806j) |
|     2140332882 |   (-0.002984855-5.88902e-39j) |            1.5737811j |           -1.5737811j | (5.889046e-39-1.5737811j) | (5.889046e-39-1.5737813j) |
|     2140284314 |  (-8.994973e-05-5.88902e-39j) |            1.5708863j |           -1.5708863j |  (5.88902e-39-1.5708864j) |  (5.88902e-39-1.5708864j) |
|     2140282850 | (-2.7069607e-06-5.88902e-39j) |             1.570799j |            -1.570799j |  (5.88902e-39-1.5707991j) |  (5.88902e-39-1.5707991j) |
|     2140282806 |  (-8.134351e-08-5.88902e-39j) |            1.5707964j |           -1.5707964j |  (5.88902e-39-1.5707965j) |  (5.88902e-39-1.5707965j) |
|     2140282806 |  (1.4992649e-24-5.88902e-39j) |            1.5707964j |           -1.5707964j |  (5.88902e-39-1.5707964j) |  (5.88902e-39-1.5707964j) |
|     2140282805 | (-4.8165032e-23-5.88902e-39j) |            1.5707963j |           -1.5707964j |  (5.88902e-39-1.5707964j) |  (5.88902e-39-1.5707964j) |
|     2140282804 |  (8.7455625e-08-5.88902e-39j) |            1.5707963j |           -1.5707963j |  (5.88902e-39-1.5707963j) |  (5.88902e-39-1.5707963j) |
|     2140282756 |  (2.9121127e-06-5.88902e-39j) |            1.5707934j |           -1.5707934j |  (5.88902e-39-1.5707935j) |  (5.88902e-39-1.5707935j) |
|     2140281180 |   (9.682065e-05-5.88902e-39j) |            1.5706995j |           -1.5706995j |  (5.88902e-39-1.5706996j) |  (5.88902e-39-1.5706996j) |
|     2140228874 |   (0.0032145183-5.88902e-39j) |            1.5675818j |           -1.5675818j |  (5.88905e-39-1.5675818j) |  (5.88905e-39-1.5675819j) |
|     2138491206 |    (0.106584825-5.88902e-39j) |            1.4640087j |           -1.4640087j | (5.922759e-39-1.4640087j) | (5.922759e-39-1.4640087j) |

Hints from issue comments:
@pearu: After disabling samples with subnormals, all tests pass on CUDA.@pearu: Closing via https://github.com/pearu/functional_algorithms/pull/56 that disables subnormal samples for JAX

Solution patch:
diff --git a/functional_algorithms/__init__.py b/functional_algorithms/__init__.py
index 288d3b6..6338d98 100644
--- a/functional_algorithms/__init__.py
+++ b/functional_algorithms/__init__.py
@@ -17,3 +17,4 @@
 from . import targets
 from .textimage import TextImage
 from . import rewrite
+from .signatures import filter_signatures
diff --git a/functional_algorithms/signatures.py b/functional_algorithms/signatures.py
new file mode 100644
index 0000000..27f2ec6
--- /dev/null
+++ b/functional_algorithms/signatures.py
@@ -0,0 +1,108 @@
+"""Provides signatures of math functions.
+"""
+
+
+def all_signatures():
+    yield from _signature_iter(_array_api_element_wise_functions)
+
+
+def filter_signatures(atypes="C", rtypes="C"):
+    atypes = tuple(atypes)
+    rtypes = tuple(rtypes)
+    for sig in all_signatures():
+        if sig["atypes"] == atypes and sig["rtypes"] == rtypes:
+            yield sig
+
+
+def _signature_iter(text):
+    for line in text.splitlines():
+        line = line.strip()
+        if line.startswith("#") or not line:
+            continue
+        left, rtype = line.split("->", 1)
+        rtypes = tuple(r.strip() for r in rtype.strip().split(","))
+        name = line[: line.find("(")].strip()
+        atypes = []
+        anames = []
+        for arg in left.rstrip()[left.find("(") + 1 : -1].split(","):
+            arg = arg.strip()
+            if arg == "/":
+                break
+            atype, aname = arg.split()
+            atype = atype.strip()
+            aname = aname.strip()
+            atypes.append(atype)
+            anames.append(anames)
+        yield dict(name=name, atypes=tuple(atypes), anames=tuple(anames), rtypes=rtypes)
+
+
+_array_api_element_wise_functions = """
+abs(F x, /) -> F
+abs(C x, /) -> F
+acos(F x, /) -> F
+acos(C x, /) -> C
+acosh(F x, /) -> F
+acosh(C x, /) -> C
+add(F x1, F x2, /) -> F
+add(C x1, C x2, /) -> C
+add(F x1, C x2, /) -> C
+add(C x1, V x2, /) -> C
+asin(F x, /) -> F
+asin(C x, /) -> C
+asinh(F x, /) -> F
+asinh(C x, /) -> C
+atan(F x, /) -> F
+atan(C x, /) -> C
+atan2(F x1, F x2, /) -> F
+atanh(F x, /) -> F
+atanh(C x, /) -> C
+conj(C x, /) -> C
+cos(F x, /) -> F
+cos(C x, /) -> C
+cosh(F x, /) -> F
+cosh(C x, /) -> C
+exp(F x, /) -> F
+exp(C x, /) -> C
+expm1(F x, /) -> F
+expm1(C x, /) -> C
+hypot(F x1, F x2, /) -> F
+imag(C x, /) -> F
+log(F x, /) -> F
+log(C x, /) -> C
+log1p(F x, /) -> F
+log1p(C x, /) -> C
+log2(F x, /) -> F
+log2(C x, /) -> C
+log10(F x, /) -> F
+log10(C x, /) -> C
+logaddexp(F x1, F x2, /) -> F
+multiply(F x1, F x2, /) -> F
+multiply(C x1, C x2, /) -> C
+multiply(F x1, C x2, /) -> C
+multiply(C x1, F x2, /) -> C
+negative(F x, /) -> F
+negative(C x, /) -> C
+positive(F x, /) -> F
+positive(C x, /) -> C
+pow(F x1, F x2, /) -> F
+pow(C x1, C x2, /) -> C
+pow(C x1, F x2, /) -> C
+pow(F x1, C x2, /) -> C
+real(C x, /) -> F
+sin(F x, /) -> F
+sin(C x, /) -> C
+sinh(F x, /) -> F
+sinh(C x, /) -> C
+square(F x, /) -> F
+square(C x, /) -> C
+sqrt(F x, /) -> F
+sqrt(C x, /) -> C
+subtract(F x1, F x2, /) -> F
+subtract(C x1, C x2, /) -> C
+subtract(F x1, C x2, /) -> C
+subtract(C x1, V x2, /) -> C
+tan(F x, /) -> F
+tan(C x, /) -> C
+tanh(F x, /) -> F
+tanh(C x, /) -> C
+"""
diff --git a/functional_algorithms/targets/numpy.py b/functional_algorithms/targets/numpy.py
index 95c5f48..80b0ced 100644
--- a/functional_algorithms/targets/numpy.py
+++ b/functional_algorithms/targets/numpy.py
@@ -1,4 +1,3 @@
-import black
 import numpy
 import sys
 import warnings
diff --git a/functional_algorithms/utils.py b/functional_algorithms/utils.py
index 5ad825e..be5c8fd 100644
--- a/functional_algorithms/utils.py
+++ b/functional_algorithms/utils.py
@@ -32,7 +32,108 @@ def __str__(self):
 default_flush_subnormals = False
 
 
-class vectorize_with_mpmath(numpy.vectorize):
+class vectorize_with_backend(numpy.vectorize):
+
+    @classmethod
+    def backend_is_available(cls, device):
+        return device == "cpu"
+
+    def __init__(self, *args, **kwargs):
+        self.device = kwargs.pop("device", "cpu")
+        super().__init__(*args, **kwargs)
+
+    @property
+    def backend_types(self):
+        raise NotImplemented(f"{type(self).__name__}.backend_types")
+
+    def backend_context(self, context):
+        raise NotImplemented(f"{type(self).__name__}.backend_context")
+
+    def numpy_to_backend(self, arr):
+        raise NotImplemented(f"{type(self).__name__}.numpy_to_backend")
+
+    def numpy_from_backend(self, obj):
+        raise NotImplemented(f"{type(self).__name__}.numpy_from_backend")
+
+    def context_from_backend(self, obj):
+        raise NotImplemented(f"{type(self).__name__}.context_from_backend")
+
+    def _call_eval(self, sample):
+        context = self.context_from_backend(sample)
+        sample = self.numpy_to_backend(sample)
+        with self.backend_context(context):
+            with warnings.catch_warnings(action="ignore"):
+                if isinstance(sample, tuple):
+                    result = super().__call__(*sample)
+                else:
+                    result = super().__call__(sample)
+        return self.numpy_from_backend(result)
+
+    def call(self, samples, workers=None, enable_progressbar=False):
+        """Apply function to samples using concurrency with the given number
+        of workers. When workers is None, use os.cpu_count() workers.
+
+        Using this method is unnecessary if the backend already
+        supports multithreading functions evaluations.
+        """
+        if workers is None:
+            workers = min(os.cpu_count(), len(samples))
+
+        if enable_progressbar:
+
+            def progress(iter, total):
+                yield from tqdm(iter, total=total)
+
+        else:
+
+            def progress(iter, total):
+                return iter
+
+        results = []
+        if workers > 1:
+            # cannot use the default `fork` start method with JAX:
+            ctx = multiprocessing.get_context("spawn" if isinstance(self, vectorize_with_jax) else "fork")
+            with ctx.Pool(workers) as p:
+                for result in p.starmap(self._call_eval, [(s,) for s in samples]):
+                    results.append(result)
+        else:
+            for sample in progress(samples, total=len(samples)):
+                results.append(self._call_eval(sample))
+        if results and isinstance(results[0], (numpy.number, numpy.ndarray)):
+            return numpy.array(results)
+        return results
+
+    def __call__(self, *args, **kwargs):
+        mp_args = []
+        context = None
+        for a in args:
+            if isinstance(a, (numpy.ndarray, numpy.floating, numpy.complexfloating)):
+                if context is None:
+                    context = self.context_from_backend(a)  # to be used as a __call__ context below
+                with self.backend_context(self.context_from_backend(a)):
+                    a = self.numpy_to_backend(a)
+                assert self.device.upper() in str(getattr(a, "device", "cpu")).upper()
+            mp_args.append(a)
+
+        with self.backend_context(context):
+            with warnings.catch_warnings(action="ignore"):
+                result = super().__call__(*mp_args, **kwargs)
+
+        if isinstance(result, tuple):
+            lst = []
+            for r in result:
+                if (isinstance(r, numpy.ndarray) and r.dtype.kind == "O") or isinstance(r, self.backend_types):
+                    r = self.numpy_from_backend(r)
+                lst.append(r)
+            return tuple(lst)
+
+        if (isinstance(result, numpy.ndarray) and result.dtype.kind == "O") or isinstance(result, self.backend_types):
+            return self.numpy_from_backend(result)
+
+        return result
+
+
+class vectorize_with_mpmath(vectorize_with_backend):
     """Same as numpy.vectorize but using mpmath backend for function evaluation."""
 
     map_float_to_complex = dict(
@@ -60,6 +161,7 @@ def __init__(self, *args, **kwargs):
         self._contexts = None
         self._contexts_inv = None
         super().__init__(*args, **kwargs)
+        assert self.device == "cpu", self.device
 
     def __getstate__(self):
         state = self.__dict__.copy()
@@ -122,7 +224,7 @@ def nptomp(self, x):
                 return ctx.make_mpf(mpmath.libmp.fnan)
             elif numpy.isfinite(x):
                 mantissa, exponent = numpy.frexp(x)
-                man = int(numpy.ldexp(mantissa, prec))
+                man = int(mpmath.ldexp(mantissa, prec))
                 exp = int(exponent - prec)
                 r = ctx.make_mpf(mpmath.libmp.from_man_exp(man, exp, prec, rounding))
                 assert ctx.isfinite(r), r._mpf_
@@ -187,78 +289,61 @@ def mptonp(self, x):
 
         raise NotImplementedError(f"convert {type(x)} instance to numpy floating point type")
 
-    def _eval(self, sample, extraprec):
-        context = self.get_context(sample)
-        sample = self.nptomp(sample)
-        with context.extraprec(extraprec):
-            if isinstance(sample, tuple):
-                result = super().__call__(*sample)
-            else:
-                result = super().__call__(sample)
-        return self.mptonp(result)
+    def context_from_backend(self, sample):
+        return self.get_context(sample)
 
-    def call(self, samples, workers=None, enable_progressbar=False):
-        """Apply function to samples using concurrency with the given number
-        of workers. When workers is None, use os.cpu_count() workers.
+    def backend_context(self, context):
+        extraprec = int(context.prec * self.extra_prec_multiplier) + self.extra_prec
+        return context.extraprec(extraprec)
 
-        """
-        if workers is None:
-            workers = min(os.cpu_count(), len(samples))
+    def numpy_to_backend(self, arr):
+        return self.nptomp(arr)
 
-        if enable_progressbar:
+    def numpy_from_backend(self, obj):
+        return self.mptonp(obj)
 
-            def progress(iter, total):
-                yield from tqdm(iter, total=total)
+    @property
+    def backend_types(self):
+        return (mpmath.ctx_mp.mpnumeric,)
 
-        else:
 
-            def progress(iter, total):
-                return iter
+class vectorize_with_jax(vectorize_with_backend):
 
-        context = self.get_context(samples[0])
-        extraprec = int(context.prec * self.extra_prec_multiplier) + self.extra_prec
-        samples_extraprec = [(s, extraprec) for s in samples]
-        results = []
-        if workers > 1:
-            with multiprocessing.Pool(workers) as p:
-                for result in p.starmap(self._eval, samples_extraprec):
-                    results.append(result)
-        else:
-            for sample, extraprec in progress(samples_extraprec, total=len(samples)):
-                results.append(self._eval(sample, extraprec))
-        if results and isinstance(results[0], (numpy.number, numpy.ndarray)):
-            return numpy.array(results)
-        return results
+    @classmethod
+    def backend_is_available(cls, device):
+        import jax
 
-    def __call__(self, *args, **kwargs):
-        mp_args = []
-        context = None
-        for a in args:
-            if isinstance(a, (numpy.ndarray, numpy.floating, numpy.complexfloating)):
-                mp_args.append(self.nptomp(a))
-                if context is None:
-                    context = self.get_context(a)
-                else:
-                    assert context is self.get_context(a), (context, self.get_context(a), a.dtype)
-            else:
-                mp_args.append(a)
+        try:
+            return bool(jax.devices(device))
+        except RuntimeError:
+            return False
 
-        extra_prec = int(context.prec * self.extra_prec_multiplier) + self.extra_prec
-        with context.extraprec(extra_prec):
-            result = super().__call__(*mp_args, **kwargs)
+    def __init__(self, *args, **kwargs):
+        self.device = kwargs.pop("device", "cpu")
+        super().__init__(*args, **kwargs)
 
-        if isinstance(result, tuple):
-            lst = []
-            for r in result:
-                if (isinstance(r, numpy.ndarray) and r.dtype.kind == "O") or isinstance(r, mpmath.ctx_mp.mpnumeric):
-                    r = self.mptonp(r)
-                lst.append(r)
-            return tuple(lst)
+    def backend_context(self, context):
+        return context
 
-        if (isinstance(result, numpy.ndarray) and result.dtype.kind == "O") or isinstance(result, mpmath.ctx_mp.mpnumeric):
-            return self.mptonp(result)
+    def numpy_to_backend(self, arr):
+        import jax
 
-        return result
+        return jax.numpy.array(arr)
+
+    def numpy_from_backend(self, obj):
+        return numpy.array(obj)
+
+    @property
+    def backend_types(self):
+        import jaxlib
+
+        return (jaxlib.xla_extension.ArrayImpl,)
+
+    def context_from_backend(self, sample):
+        import jax
+
+        jax.config.update("jax_enable_x64", sample.dtype.type in {numpy.float64, numpy.complex128})
+        return jax.default_device(jax.devices(self.device)[0])
 
 
 class mpmath_array_api:
@@ -266,6 +351,9 @@ class mpmath_array_api:
     mpmath bugs.
     """
 
+    def square(self, x):
+        return x * x
+
     def hypot(self, x, y):
         return x.context.hypot(x, y)
 
@@ -314,6 +402,8 @@ def cosh(self, x):
     def conjugate(self, x):
         return x.context.conjugate(x)
 
+    conj = conjugate
+
     def sign(self, x):
         return x.context.sign(x)
 
@@ -678,6 +768,82 @@ def worker(ctx, s, e, r, v):
             assert 0  # unreachable
 
 
+class numpy_with_jax:
+    """Namespace of universal functions on numpy arrays that use jax
+    backend for evaluation and return numpy arrays as outputs.
+    """
+
+    _vfunc_cache = dict()
+
+    def __init__(self, **params):
+        self.params = params
+
+    def __getattr__(self, 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]
+        import jax
+
+        vfunc = vectorize_with_jax(getattr(jax.numpy, name), **self.params)
+        self._vfunc_cache[name] = vfunc
+        return vfunc
+
+
+class numpy_with_numpy:
+    """Namespace of universal functions on numpy arrays that use numpy
+    backend for evaluation and return numpy arrays as outputs.
+    """
+
+    _vfunc_cache = dict()
+
+    def __init__(self, **params):
+        self.params = params
+
+    def __getattr__(self, 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]
+        import numpy
+
+        vfunc = numpy.vectorize(getattr(numpy, name), **self.params)
+        self._vfunc_cache[name] = vfunc
+        return vfunc
+
+
+class numpy_with_algorithms:
+    """Namespace of universal functions on numpy arrays that use
+    functional_algorithms definitions of algorithms.
+    """
+
+    _vfunc_cache = dict()
+
+    def __init__(self, **params):
+        self.params = params
+
+    def __getattr__(self, name):
+        name = dict(arcsinh="asinh", arccos="acos", arcsin="asin", arccosh="acosh", arctan="atan", arctanh="atanh").get(
+            name, name
+        )
+        if name in self._vfunc_cache:
+            return self._vfunc_cache[name]
+
+        import functional_algorithms as fa
+        import numpy
+
+        dtype = self.params["dtype"]
+        ctx = fa.Context(paths=[fa.algorithms])
+        graph = ctx.trace(getattr(fa.algorithms, name), dtype)
+        graph2 = graph.rewrite(fa.targets.numpy, fa.rewrite)
+        func = fa.targets.numpy.as_function(graph2, debug=0)
+        vfunc = numpy.vectorize(func)
+        self._vfunc_cache[name] = vfunc
+        return vfunc
+
+
 def extra_samples(name, dtype):
     """Return a list of samples that are special to a given function.
 
@@ -814,7 +980,7 @@ def real_samples(
         raise ValueError(f"minimal value (={min_value}) cannot be greater than maximal value (={max_value})")
 
     if user_specified_bounds:
-        if min_value > 0:
+        if min_value >= 0:
             start, end = min_value.view(utype), max_value.view(utype)
             step = int(end - start)
             r = (start + numpy.array([i // (num - 1) for i in range(0, num * step, step)], dtype=utype)).view(dtype)
@@ -1337,6 +1503,12 @@ def function_validation_parameters(func_name, dtype):
     else:
         dtype_name = dtype.__name__
 
+    # If a function has symmetries, exclude superfluous samples by
+    # specifying a region of the function domain:
+    samples_limits = dict(
+        min_real_value=-numpy.inf, max_real_value=numpy.inf, min_imag_value=-numpy.inf, max_imag_value=numpy.inf
+    )
+
     # ulp_diff(func(sample), reference(sample)) <= max_valid_ulp_count
     max_valid_ulp_count = 3
 
@@ -1354,7 +1526,7 @@ def function_validation_parameters(func_name, dtype):
     # precision of mpmath to improve the accurancy of reference
     # values:
     extra_prec_multiplier = 1
-    if func_name in {"acos", "asin", "asinh", "acosh"}:
+    if func_name in {"acos", "asin", "asinh", "acosh", "atan", "atanh"}:
         extra_prec_multiplier = 20
     elif func_name == "sqrt":
         extra_prec_multiplier = 2
@@ -1369,6 +1541,7 @@ def function_validation_parameters(func_name, dtype):
         extra_prec_multiplier=extra_prec_multiplier,
         max_valid_ulp_count=max_valid_ulp_count,
         max_bound_ulp_width=max_bound_ulp_width,
+        samples_limits=samples_limits,
     )
 
 
@@ -1423,3 +1596,68 @@ def warn_once(msg, stacklevel=1):
         return
     _warn_once_cache.add(msg)
     warnings.warn(msg, stacklevel=stacklevel + 1)
+
+
+class Clusters:
+
+    class Cluster:
+
+        def __init__(self):
+            self.points = set()
+
+        def contains(self, point):
+            for x, y in self.points:
+                for dx in {-1, 0, 1}:
+                    for dy in {-1, 0, 1}:
+                        if point == (x + dx, y + dy):
+                            return True
+            return False
+
+        def add(self, point):
+            self.points.add(point)
+
+        def merge_from(self, other):
+            self.points.update(other.points)
+
+        def __repr__(self):
+            return f"{type(self).__name__}({self.points})"
+
+        def center(self):
+            sx, sy = 0, 0
+            for x, y in self.points:
+                sx += x
+                sy += y
+            return (sx / len(self.points), sy / len(self.points))
+
+        def center_point(self):
+            cx, cy = self.center()
+            lst = []
+            for x, y in self.points:
+                lst.append((abs(x - cx) + abs(y - cy), (x, y)))
+            return sorted(lst)[0][1]
+
+    def __init__(self):
+        self.clusters = []
+
+    def __repr__(self):
+        return f"{type(self).__name__}({self.clusters})"
+
+    def add(self, point):
+        matching_clusters = []
+        for index, cluster in enumerate(self.clusters):
+            if cluster.contains(point):
+                matching_clusters.append(index)
+
+        if len(matching_clusters) == 0:
+            self.clusters.append(Clusters.Cluster())
+            self.clusters[-1].add(point)
+        elif len(matching_clusters) == 1:
+            self.clusters[matching_clusters[0]].add(point)
+        else:
+            cluster = self.clusters[matching_clusters[0]]
+            for index in reversed(matching_clusters[1:]):
+                cluster.merge_from(self.clusters[index])
+                del self.clusters[index]
+
+    def __iter__(self):
+        return iter(self.clusters)


Test patch:
diff --git a/functional_algorithms/tests/test_accuracy.py b/functional_algorithms/tests/test_accuracy.py
new file mode 100644
index 0000000..0b3d31c
--- /dev/null
+++ b/functional_algorithms/tests/test_accuracy.py
@@ -0,0 +1,155 @@
+import numpy
+import functional_algorithms as fa
+import pytest
+import warnings
+
+
+@pytest.fixture(scope="function", params=["jax"])
+def backend(request):
+    return request.param
+
+
+@pytest.fixture(scope="function", params=["cpu", "cuda"])
+def device(request):
+    return request.param
+
+
+@pytest.fixture(scope="function", params=list(sig["name"] for sig in fa.filter_signatures("C", "C")))
+def unary_func_name(request):
+    return request.param
+
+
+@pytest.fixture(scope="function", params=[numpy.float32, numpy.float64, numpy.complex64, numpy.complex128])
+def dtype(request):
+    return request.param
+
+
+def test_unary(unary_func_name, backend, device, dtype):
+    pytest.importorskip(backend)
+
+    func = getattr(getattr(fa.utils, f"numpy_with_{backend}")(device=device), unary_func_name)
+
+    if not func.backend_is_available(device):
+        pytest.skip(f"{device} support is unavailable")
+    params = fa.utils.function_validation_parameters(unary_func_name, dtype)
+
+    max_valid_ulp_count = params["max_valid_ulp_count"]
+    extra_prec_multiplier = params["extra_prec_multiplier"]
+    samples_limits = params["samples_limits"]
+    include_subnormal = True
+
+    mpmath = fa.utils.numpy_with_mpmath(extra_prec_multiplier=extra_prec_multiplier)
+
+    reference = getattr(mpmath, unary_func_name)
+    npy_reference = getattr(fa.utils.numpy_with_numpy(), unary_func_name)
+
+    re_blocks, im_blocks = 101, 51
+    re_blocksize, im_blocksize = 20, 20
+    re_size, im_size = re_blocks * re_blocksize, im_blocks * im_blocksize
+
+    if dtype in {numpy.complex64, numpy.complex128}:
+        samples = fa.utils.complex_samples(
+            (re_size, im_size), dtype=dtype, include_subnormal=include_subnormal, **samples_limits
+        )
+    else:
+        pytest.skip(f"{dtype} support not implemented")
+    assert samples.shape == (im_size, re_size)
+
+    expected = reference.call(samples)
+    result = func(samples)
+
+    ulp = fa.utils.diff_ulp(result, expected, flush_subnormals=not include_subnormal)
+
+    if numpy.all(ulp == 0):
+        return
+    if numpy.all(ulp <= max_valid_ulp_count):
+        print(f"maximal ULP difference: {ulp.max()}")
+        return
+
+    bsamples = numpy.zeros((im_blocks, re_blocks), dtype=samples.dtype)
+    bulp = numpy.zeros((im_blocks, re_blocks), dtype=ulp.dtype)
+    for j, blocks in enumerate(numpy.split(ulp, im_blocks, axis=0)):
+        for i, block in enumerate(numpy.split(blocks, re_blocks, axis=1)):
+            ind = numpy.unravel_index(numpy.argmax(block, axis=None), block.shape)
+            j_, i_ = j * im_blocksize + ind[0], i * re_blocksize + ind[1]
+            bsamples[j, i] = samples[j_, i_]
+            bulp[j, i] = ulp[j_, i_]
+
+    try:
+        fa_reference = getattr(fa.utils.numpy_with_algorithms(dtype=dtype), unary_func_name)
+    except Exception as msg:
+        print(f"disabling functional_algorithms output: {msg}")
+        fa_reference = None
+
+    timage = fa.TextImage()
+
+    timage.fill(0, 10, bulp == 0, symbol="=")
+    for i in range(1, max_valid_ulp_count + 1):
+        timage.fill(0, 10, bulp == i, symbol=str(i))
+
+    timage.fill(0, 10, bulp > max_valid_ulp_count, symbol="!")
+    timage.fill(0, 10, bulp > 100 * max_valid_ulp_count, symbol="E")
+
+    real_axis = samples[0:1, ::im_blocksize].real
+    imag_axis = samples[::re_blocksize, 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()
+    print(timage)
+    print()
+
+    if fa_reference is not None:
+        rows = [
+            (
+                "ULP-difference",
+                "z",
+                f"{backend}:{unary_func_name}(z)",
+                f"mpmath:{unary_func_name}(z)",
+                f"numpy:{unary_func_name}(z)",
+                f"fa:{unary_func_name}(z)",
+            )
+        ]
+    else:
+        rows = [
+            (
+                "ULP-difference",
+                "z",
+                f"{backend}:{unary_func_name}(z)",
+                f"mpmath:{unary_func_name}(z)",
+                f"numpy:{unary_func_name}(z)",
+            )
+        ]
+
+    samples = bsamples
+    ulp = bulp
+
+    max_rows = 20
+    for value in reversed(sorted(set(ulp.flatten()))):
+        if value <= max_valid_ulp_count:
+            break
+        clusters = fa.utils.Clusters()
+        for re, im in zip(*numpy.where(ulp == value)):
+            clusters.add((re, im))
+        for cluster in clusters:
+            re, im = cluster.center_point()
+            with warnings.catch_warnings(action="ignore"):
+                np_value = npy_reference(samples[re, im])
+            r = func(samples[re, im])
+            e = reference(samples[re, im])
+            if fa_reference is not None:
+                fa_value = fa_reference(samples[re, im])
+                rows.append((value, samples[re, im], r, e, np_value, fa_value))
+            else:
+                rows.append((value, samples[re, im], r, e, np_value))
+            if len(rows) > max_rows:
+                break
+        if len(rows) > max_rows:
+            rows.append(("...",) * len(rows[0]))
+            break
+    col_widths = [max(len(str(row[col])) for row in rows) for col in range(len(rows[0]))]
+    rows.insert(1, tuple("-" * w for w in col_widths))
+    col_fmt = "| " + " | ".join([f"{{{i}:>{w}}}" for i, w in enumerate(col_widths)]) + " |"
+    print("\n".join([col_fmt.format(*map(str, row)) for row in rows]))
+
+    pytest.xfail("inaccurate or incorrect results")
diff --git a/functional_algorithms/tests/test_functional_algorithms.py b/functional_algorithms/tests/test_functional_algorithms.py
index 246164d..21ec47d 100644
--- a/functional_algorithms/tests/test_functional_algorithms.py
+++ b/functional_algorithms/tests/test_functional_algorithms.py
@@ -240,16 +240,16 @@ def test_readme_square_numpy_debug_0():
     assert py == utils.format_python(
         """\
 def readme_square(z: numpy.complex128) -> numpy.complex128:
-    with warnings.catch_warnings(action="ignore"):
-        z = numpy.complex128(z)
-        real_z: numpy.float64 = (z).real
-        x: numpy.float64 = numpy.abs(real_z)
-        y: numpy.float64 = numpy.abs((z).imag)
-        real_part: numpy.float64 = ((x) - (y)) * ((y) + (y))
-        real: numpy.float64 = (numpy.float64(0)) if (numpy.equal(x, y, dtype=numpy.bool_)) else (real_part)
-        imag: numpy.float64 = (numpy.float64(2)) * ((x) * (y))
-        result = make_complex(real, imag)
-        return result"""
+  with warnings.catch_warnings(action="ignore"):
+    z = numpy.complex128(z)
+    real_z: numpy.float64 = (z).real
+    x: numpy.float64 = numpy.abs(real_z)
+    y: numpy.float64 = numpy.abs((z).imag)
+    real_part: numpy.float64 = ((x) - (y)) * ((y) + (y))
+    real: numpy.float64 = (numpy.float64(0)) if (numpy.equal(x, y, dtype=numpy.bool_)) else (real_part)
+    imag: numpy.float64 = (numpy.float64(2)) * ((x) * (y))
+    result = make_complex(real, imag)
+    return result"""
     )
 
 
diff --git a/functional_algorithms/tests/test_utils.py b/functional_algorithms/tests/test_utils.py
index 7b5cd08..7068d11 100644
--- a/functional_algorithms/tests/test_utils.py
+++ b/functional_algorithms/tests/test_utils.py
@@ -6,6 +6,11 @@
 
 
 @pytest.fixture(scope="function", params=[numpy.float32, numpy.float64])
+def real_dtype(request):
+    return request.param
+
+
+@pytest.fixture(scope="function", params=[numpy.float32, numpy.float64, numpy.complex64, numpy.complex128])
 def dtype(request):
     return request.param
 
@@ -109,35 +114,35 @@ def _iter_samples_parameters():
         )
 
 
-def test_real_samples(dtype):
+def test_real_samples(real_dtype):
     for size in range(6, 20):
         for params in _iter_samples_parameters():
-            r = utils.real_samples(size=size, dtype=dtype, **params)
-            assert r.dtype == dtype
+            r = utils.real_samples(size=size, dtype=real_dtype, **params)
+            assert r.dtype == real_dtype
             _check_real_samples(r, **params)
 
 
-def test_complex_samples(dtype):
+def test_complex_samples(real_dtype):
     for size in [(6, 6), (6, 7), (7, 6), (13, 13), (13, 15), (15, 13)]:
         for params in _iter_samples_parameters():
-            r = utils.complex_samples(size=size, dtype=dtype, **params)
+            r = utils.complex_samples(size=size, dtype=real_dtype, **params)
             re, im = r.real, r.imag
-            assert re.dtype == dtype
-            assert im.dtype == dtype
+            assert re.dtype == real_dtype
+            assert im.dtype == real_dtype
             for i in range(r.shape[0]):
                 _check_real_samples(re[i], **params)
             for j in range(r.shape[1]):
                 _check_real_samples(im[:, j], **params)
 
 
-def test_real_pair_samples(dtype):
+def test_real_pair_samples(real_dtype):
     for size in [(6, 6), (6, 7), (7, 6), (13, 13), (13, 15), (15, 13)]:
         for params in _iter_samples_parameters():
-            s1 = utils.real_samples(size=size[0], dtype=dtype, **params).size
-            s2 = utils.real_samples(size=size[1], dtype=dtype, **params).size
-            r1, r2 = utils.real_pair_samples(size=size, dtype=dtype, **params)
-            assert r1.dtype == dtype
-            assert r2.dtype == dtype
+            s1 = utils.real_samples(size=size[0], dtype=real_dtype, **params).size
+            s2 = utils.real_samples(size=size[1], dtype=real_dtype, **params).size
+            r1, r2 = utils.real_pair_samples(size=size, dtype=real_dtype, **params)
+            assert r1.dtype == real_dtype
+            assert r2.dtype == real_dtype
             assert r1.size == s1 * s2
             assert r2.size == s1 * s2
             for r in r1.reshape(s2, s1):
@@ -146,17 +151,17 @@ def test_real_pair_samples(dtype):
                 _check_real_samples(r, **params)
 
 
-def test_complex_pair_samples(dtype):
+def test_complex_pair_samples(real_dtype):
     for size1 in [(6, 6), (6, 7)]:
         for size2 in [(6, 6), (7, 6)]:
             for params in _iter_samples_parameters():
-                s1 = utils.complex_samples(size=size1, dtype=dtype, **params).shape
-                s2 = utils.complex_samples(size=size2, dtype=dtype, **params).shape
-                r1, r2 = utils.complex_pair_samples(size=(size1, size2), dtype=dtype, **params)
+                s1 = utils.complex_samples(size=size1, dtype=real_dtype, **params).shape
+                s2 = utils.complex_samples(size=size2, dtype=real_dtype, **params).shape
+                r1, r2 = utils.complex_pair_samples(size=(size1, size2), dtype=real_dtype, **params)
                 re1, im1 = r1.real, r1.imag
                 re2, im2 = r2.real, r2.imag
-                assert re1.dtype == dtype
-                assert re2.dtype == dtype
+                assert re1.dtype == real_dtype
+                assert re2.dtype == real_dtype
                 assert r1.shape == (s1[0] * s2[0], s1[1] * s2[1])
                 assert r2.shape == (s1[0] * s2[0], s1[1] * s2[1])
 
@@ -179,12 +184,12 @@ def test_complex_pair_samples(dtype):
                             _check_real_samples(r[:, j1], **params)
 
 
-def test_periodic_samples(dtype):
+def test_periodic_samples(real_dtype):
     for period in [1, 3.14, numpy.pi / 2]:
         for size in [10, 11, 51]:
             for periods in [4, 5, 13]:
-                samples = utils.periodic_samples(period=period, size=size, dtype=dtype, periods=periods)
-                assert samples.dtype == dtype
+                samples = utils.periodic_samples(period=period, size=size, dtype=real_dtype, periods=periods)
+                assert samples.dtype == real_dtype
                 assert samples.size == size * periods
                 assert numpy.diff(samples).min() >= 0
                 assert numpy.diff(samples).max() > 0
@@ -206,3 +211,67 @@ def test_periodic_samples(dtype):
 
                     p = samples_per_period[-1] - samples_per_period[0]
                     assert p <= period * size / (size - 1) and p >= period * (size - 1) / size
+
+
+@pytest.fixture(scope="function", params=["mpmath", "jax"])
+def backend(request):
+    return request.param
+
+
+@pytest.fixture(scope="function", params=["cpu", "cuda"])
+def device(request):
+    return request.param
+
+
+def _square(x):
+    return x * x
+
+
+def test_vectorize_with_backend(backend, dtype, device):
+    pytest.importorskip(backend)
+
+    assert_equal = numpy.testing.assert_equal
+
+    vectorize_with_backend = getattr(utils, f"vectorize_with_{backend}")
+
+    if not vectorize_with_backend.backend_is_available(device):
+        pytest.skip(f"{device} support is unavailable")
+
+    func = vectorize_with_backend(_square, device=device)
+
+    arr = numpy.array([1, 2, 3], dtype=dtype)
+    expected = _square(arr)
+
+    def test(result):
+        assert isinstance(result, numpy.ndarray)
+        assert result.dtype == expected.dtype
+        assert_equal(result, expected)
+
+    test(func(arr))
+    # Using call for JAX backend is slow and unreasonable, so
+    # multiprocessing is disabled for JAX:
+    test(func.call(arr, workers=1 if backend == "jax" else None))
+
+
+def test_numpy_with_backend(backend, dtype, device):
+    pytest.importorskip(backend)
+
+    assert_equal = numpy.testing.assert_equal
+
+    numpy_with_backend = getattr(utils, f"numpy_with_{backend}")(device=device)
+
+    func = numpy_with_backend.square
+
+    if not func.backend_is_available(device):
+        pytest.skip(f"{device} support is unavailable")
+
+    arr = numpy.array([1, 2, 3], dtype=dtype)
+
+    expected = numpy.square(arr)
+
+    def test(result):
+        assert isinstance(result, numpy.ndarray)
+        assert result.dtype == expected.dtype
+        assert_equal(result, expected)
+
+    test(func(arr))