diff --git a/tests/helpers.py b/tests/helpers.py
index 7dc7fcea..27155442 100644
--- a/tests/helpers.py
+++ b/tests/helpers.py
@@ -191,6 +191,18 @@ def numbers_close(x, y, tolerance=1e-6):
         return abs(x or y) < tolerance
 
 
+def nominal_and_std_dev_close(x, y, tolerance=1e-6):
+    """
+    Tests if two numbers with uncertainties are close, NOT as random
+    variables. Checks whether the magnitude of the nominal
+    values and standard deviations are close.
+
+    The tolerance is applied to both the nominal value and the
+    standard deviation of the difference between the numbers.
+    """
+    return numbers_close(x.n, y.n, tolerance) and numbers_close(x.s, y.s, tolerance)
+
+
 def ufloats_close(x, y, tolerance=1e-6):
     """
     Tests if two numbers with uncertainties are close, as random
diff --git a/tests/test_ufloatnumpy.py b/tests/test_ufloatnumpy.py
new file mode 100644
index 00000000..75eafd71
--- /dev/null
+++ b/tests/test_ufloatnumpy.py
@@ -0,0 +1,247 @@
+from uncertainties import umath, ufloat
+from helpers import nominal_and_std_dev_close
+import numpy as np
+import pytest
+
+a = ufloat(1, 0.1)
+b = ufloat(2, 0.2)
+
+
+class TestArithmetic:
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(3.0, 0.223606797749979)),
+            (a, a, ufloat(2.0, 0.2)),
+        ],
+    )
+    def test_add(self, first, second, expected):
+        result = first + second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.add(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(-1.00, 0.223606797749979)),
+            (a, a, ufloat(0.0, 0.0)),
+        ],
+    )
+    def test_subtact(self, first, second, expected):
+        result = first - second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.subtract(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(2.0, 0.28284271247461906)),
+            (a, a, ufloat(1.0, 0.2)),
+        ],
+    )
+    def test_multiply(self, first, second, expected):
+        result = first * second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.multiply(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(0.5, 0.07071067811865477)),
+            (a, a, ufloat(1.0, 0.0)),
+        ],
+    )
+    def test_divide(self, first, second, expected):
+        result = first / second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.true_divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(0.0, 0.0)),
+            (a, a, ufloat(1.0, 0.0)),
+        ],
+    )
+    def test_floor_divide(self, first, second, expected):
+        result = first // second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.floor_divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+
+class TestComparative:
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, True),
+        ],
+    )
+    def test_equal(self, first, second, expected):
+        result = first == second
+        assert result == expected
+
+        result = np.equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, False),
+        ],
+    )
+    def test_not_equal(self, first, second, expected):
+        result = first != second
+        assert result == expected
+
+        result = np.not_equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, False),
+        ],
+    )
+    def test_less(self, first, second, expected):
+        result = first < second
+        assert result == expected
+
+        result = np.less(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, True),
+        ],
+    )
+    def test_less_equal(self, first, second, expected):
+        result = first <= second
+        assert result == expected
+
+        result = np.less_equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, False),
+        ],
+    )
+    def test_greater(self, first, second, expected):
+        result = first > second
+        assert result == expected
+
+        result = np.greater(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, True),
+        ],
+    )
+    def test_greater_equal(self, first, second, expected):
+        result = first >= second
+        assert result == expected
+
+        result = np.greater_equal(first, second)
+        assert result == expected
+
+
+class TestUfuncs:
+    zero = ufloat(0.0, 0.1)
+    one = ufloat(1.0, 0.1)
+    pi_4 = ufloat(0.7853981633974483, 0.1)  # pi/4
+    pi_2 = ufloat(1.5707963267948966, 0.1)  # pi/2
+
+    @pytest.mark.parametrize(
+        "numpy_func, umath_func, arg, expected",
+        [
+            ("cos", "cos", zero, ufloat(1.0, 0.0)),
+            ("cos", "cos", pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
+            ("cos", "cos", pi_2, ufloat(6.123233995736766e-17, 0.1)),
+            ("cosh", "cosh", zero, ufloat(1.0, 0.0)),
+            ("cosh", "cosh", pi_4, ufloat(1.324609089252006, 0.08686709614860096)),
+            ("cosh", "cosh", pi_2, ufloat(2.5091784786580567, 0.2301298902307295)),
+            ("sin", "sin", zero, ufloat(0.0, 0.1)),
+            ("sin", "sin", pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
+            ("sin", "sin", pi_2, ufloat(1.0, 6.123233995736766e-18)),
+            ("sinh", "sinh", zero, ufloat(0.0, 0.1)),
+            ("sinh", "sinh", pi_4, ufloat(0.8686709614860095, 0.1324609089252006)),
+            ("sinh", "sinh", pi_2, ufloat(2.3012989023072947, 0.2509178478658057)),
+            ("tan", "tan", zero, ufloat(0.0, 0.1)),
+            ("tan", "tan", pi_4, ufloat(0.9999999999999999, 0.19999999999999998)),
+            ("tan", "tan", pi_2, ufloat(1.633123935319537e16, 2.6670937881135717e31)),
+            ("tanh", "tanh", zero, ufloat(0.0, 0.1)),
+            ("tanh", "tanh", pi_4, ufloat(0.6557942026326724, 0.05699339637933774)),
+            ("tanh", "tanh", pi_2, ufloat(0.9171523356672744, 0.015883159318006324)),
+            ("arccos", "acos", zero, ufloat(1.5707963267948966, 0.1)),
+            ("arccos", "acos", one, ufloat(0.0, float("nan"))),
+            ("arccosh", "acosh", one, ufloat(0.0, float("nan"))),
+            ("arcsin", "asin", zero, ufloat(0.0, 0.1)),
+            ("arcsin", "asin", one, ufloat(1.5707963267948966, float("nan"))),
+            ("arcsinh", "asinh", zero, ufloat(0.0, 0.1)),
+            ("arcsinh", "asinh", one, ufloat(0.8813735870195429, 0.07071067811865475)),
+            ("arctan", "atan", zero, ufloat(0.0, 0.1)),
+            ("arctan", "atan", one, ufloat(0.7853981633974483, 0.05)),
+            ("arctanh", "atanh", zero, ufloat(0.0, 0.1)),
+            ("exp", "exp", zero, ufloat(1.0, 0.1)),
+            ("exp", "exp", one, ufloat(2.718281828459045, 0.27182818284590454)),
+            ("exp2", None, zero, ufloat(1.0, 0.06931471805599453)),
+            ("exp2", None, one, ufloat(2.0, 0.13862943611198905)),
+            ("expm1", "expm1", zero, ufloat(0.0, 0.1)),
+            ("expm1", "expm1", one, ufloat(1.718281828459045, 0.27182818284590454)),
+            ("log10", "log10", one, ufloat(0.0, 0.04342944819032518)),
+            ("log1p", "log1p", zero, ufloat(0.0, 0.1)),
+            ("log1p", "log1p", one, ufloat(0.6931471805599453, 0.05)),
+            ("degrees", "degrees", zero, ufloat(0.0, 5.729577951308233)),
+            ("degrees", "degrees", one, ufloat(57.29577951308232, 5.729577951308233)),
+            ("radians", "radians", zero, ufloat(0.0, 0.0017453292519943296)),
+            (
+                "radians",
+                "radians",
+                one,
+                ufloat(0.017453292519943295, 0.0017453292519943296),
+            ),
+            ("rad2deg", "degrees", zero, ufloat(0.0, 5.729577951308233)),
+            ("rad2deg", "degrees", one, ufloat(57.29577951308232, 5.729577951308233)),
+            ("deg2rad", "radians", zero, ufloat(0.0, 0.0017453292519943296)),
+            (
+                "deg2rad",
+                "radians",
+                one,
+                ufloat(0.017453292519943295, 0.0017453292519943296),
+            ),
+            ("sqrt", "sqrt", zero, ufloat(0.0, float("nan"))),
+            ("sqrt", "sqrt", one, ufloat(1.0, 0.05)),
+        ],
+    )
+    def test_single_arg(self, numpy_func, umath_func, arg, expected):
+        func = getattr(np, numpy_func)
+        result = func(arg)
+        assert nominal_and_std_dev_close(result, expected)
+
+        if umath_func:
+            func = getattr(umath, umath_func)
+            result = func(arg)
+            assert nominal_and_std_dev_close(result, expected)
diff --git a/uncertainties/core.py b/uncertainties/core.py
index fdd64db7..82e9e849 100644
--- a/uncertainties/core.py
+++ b/uncertainties/core.py
@@ -39,6 +39,8 @@ def isinfinite(x):
     modified_operators,
     modified_ops_with_reflection,
 )
+from .ufloatnumpy import UFloatNumpy
+
 
 # Attributes that are always exported (some other attributes are
 # exported only if the NumPy module is available...):
@@ -332,7 +334,15 @@ def __setstate__(self, state):
         (self.linear_combo,) = state
 
 
-class AffineScalarFunc(object):
+try:
+    import numpy
+except ImportError:
+    parent_classes = []
+else:
+    parent_classes = [UFloatNumpy]
+
+
+class AffineScalarFunc(*parent_classes):
     """
     Affine functions that support basic mathematical operations
     (addition, etc.).  Such functions can for instance be used for
@@ -657,6 +667,8 @@ def __setstate__(self, data_dict):
 ops.add_arithmetic_ops(AffineScalarFunc)
 ops.add_comparative_ops(AffineScalarFunc)
 to_affine_scalar = AffineScalarFunc._to_affine_scalar
+AffineScalarFunc._add_numpy_arithmetic_ufuncs()
+AffineScalarFunc._add_numpy_comparative_ufuncs()
 
 # Nicer name, for users: isinstance(ufloat(...), UFloat) is
 # True. Also: isinstance(..., UFloat) is the test for "is this a
diff --git a/uncertainties/ufloatnumpy.py b/uncertainties/ufloatnumpy.py
new file mode 100644
index 00000000..5373d535
--- /dev/null
+++ b/uncertainties/ufloatnumpy.py
@@ -0,0 +1,236 @@
+import numpy as np
+import math
+
+# ufuncs are listed at https://numpy.org/doc/stable/reference/ufuncs.html
+from . import ops
+# from .umath_core import log_der0,_deriv_copysign, _deriv_fabs, _deriv_pow_0, _deriv_pow_1
+
+from .ops import nan_if_exception
+
+
+def log_der0(*args):
+    """
+    Derivative of math.log() with respect to its first argument.
+
+    Works whether 1 or 2 arguments are given.
+    """
+    if len(args) == 1:
+        return 1 / args[0]
+    else:
+        return 1 / args[0] / math.log(args[1])  # 2-argument form
+
+    # The following version goes about as fast:
+
+    ## A 'try' is used for the most common case because it is fast when no
+    ## exception is raised:
+    # try:
+    #    return log_1arg_der(*args)  # Argument number check
+    # except TypeError:
+    #    return 1/args[0]/math.log(args[1])  # 2-argument form
+
+
+def _deriv_copysign(x, y):
+    if x >= 0:
+        return math.copysign(1, y)
+    else:
+        return -math.copysign(1, y)
+
+
+def _deriv_fabs(x):
+    if x >= 0:
+        return 1
+    else:
+        return -1
+
+
+def _deriv_pow_0(x, y):
+    if y == 0:
+        return 0.0
+    elif x != 0 or y % 1 == 0:
+        return y * math.pow(x, y - 1)
+    else:
+        return float("nan")
+
+
+def _deriv_pow_1(x, y):
+    if x == 0 and y > 0:
+        return 0.0
+    else:
+        return math.log(x) * math.pow(x, y)
+
+
+def is_upcast_type(t):
+    # This can be used to allow downstream modules to overide operations; see pint
+    # TODO add upcast_type list or dict to a public interface
+    return False
+
+
+def implements(numpy_func_string, func_type):
+    """Register an __array_function__/__array_ufunc__ implementation for UArray
+    objects.
+    """
+    print(numpy_func_string, func_type)
+
+    def decorator(func):
+        if func_type == "function":
+            HANDLED_FUNCTIONS[numpy_func_string] = func
+        elif func_type == "ufunc":
+            HANDLED_UFUNCS[numpy_func_string] = func
+        else:
+            raise ValueError(f"Invalid func_type {func_type}")
+        return func
+
+    return decorator
+
+
+HANDLED_UFUNCS = {}
+HANDLED_FUNCTIONS = {}
+
+
+def apply_func_elementwise(func, inputs, kwargs, result_dtype="object"):
+    if len(inputs) == 1:
+        result = func(*inputs, **kwargs)
+    elif isinstance(inputs[0], np.ndarray):
+        result = np.empty_like(inputs[0], dtype=result_dtype)
+        for index, x in np.ndenumerate(inputs[0]):
+            inputs_ = [x if i == 0 else inputs[i] for i in range(len(inputs))]
+            result[index] = func(*inputs_, **kwargs)
+        # unpack the result of operations with ndim=0 arrays
+        if inputs[0].ndim == 0:
+            result = result.item()
+    elif isinstance(inputs[1], np.ndarray):
+        result = np.empty_like(inputs[1], dtype=result_dtype)
+        for index, x in np.ndenumerate(inputs[1]):
+            inputs_ = [x if i == 1 else inputs[i] for i in range(len(inputs))]
+            result[index] = func(*inputs_, **kwargs)
+        # unpack the result of operations with ndim=0 arrays
+        if inputs[1].ndim == 0:
+            result = result.item()
+    else:
+        result = func(*inputs, **kwargs)
+    return result
+
+
+def numpy_wrap(func_type, func, args, kwargs, types):
+    """Return the result from a NumPy function/ufunc as wrapped by uncertainties."""
+
+    if func_type == "function":
+        handled = HANDLED_FUNCTIONS
+        # Need to handle functions in submodules
+        name = ".".join(func.__module__.split(".")[1:] + [func.__name__])
+    elif func_type == "ufunc":
+        handled = HANDLED_UFUNCS
+        # ufuncs do not have func.__module__
+        name = func.__name__
+    else:
+        raise ValueError(f"Invalid func_type {func_type}")
+
+    if name not in handled or any(is_upcast_type(t) for t in types):
+        return NotImplemented
+    return handled[name](*args, **kwargs)
+
+
+class UFloatNumpy(object):
+    # NumPy function/ufunc support
+    __array_priority__ = 17
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        if method != "__call__":
+            # Only handle ufuncs as callables
+            return NotImplemented
+
+        # Replicate types from __array_function__
+        types = {
+            type(arg)
+            for arg in list(inputs) + list(kwargs.values())
+            if hasattr(arg, "__array_ufunc__")
+        }
+
+        return numpy_wrap("ufunc", ufunc, inputs, kwargs, types)
+
+    def __array_function__(self, func, types, args, kwargs):
+        return numpy_wrap("function", func, args, kwargs, types)
+
+    @classmethod
+    def _add_numpy_arithmetic_ufuncs(cls):
+        def implement_ufunc(func_str, derivatives):
+            func = getattr(np, func_str)
+
+            @implements(func_str, "ufunc")
+            def implementation(*inputs, **kwargs):
+                return apply_func_elementwise(
+                    ops._wrap(cls, func, derivatives), inputs, kwargs
+                )
+
+            return implementation
+
+        ufunc_derivatives = {
+            "add": [lambda x, y: 1.0, lambda x, y: 1.0],
+            "subtract": [lambda x, y: 1.0, lambda x, y: -1.0],
+            "multiply": [lambda x, y: y, lambda x, y: x],
+            "divide": [lambda x, y: 1.0 / y, lambda x, y: -x / y**2],
+            "true_divide": [lambda x, y: 1.0 / y, lambda x, y: -x / y**2],
+            "floor_divide": [lambda x, y: 0.0, lambda x, y: 0.0],
+            "arccos": [nan_if_exception(lambda x: -1 / math.sqrt(1 - x**2))],
+            "arccosh": [nan_if_exception(lambda x: 1 / math.sqrt(x**2 - 1))],
+            "arcsin": [nan_if_exception(lambda x: 1 / math.sqrt(1 - x**2))],
+            "arcsinh": [nan_if_exception(lambda x: 1 / math.sqrt(1 + x**2))],
+            "arctan": [nan_if_exception(lambda x: 1 / (1 + x**2))],
+            "arctan2": [
+                nan_if_exception(lambda y, x: x / (x**2 + y**2)),  # Correct for x == 0
+                nan_if_exception(lambda y, x: -y / (x**2 + y**2)),
+            ],  # Correct for x == 0
+            "arctanh": [nan_if_exception(lambda x: 1 / (1 - x**2))],
+            "cos": [lambda x: -math.sin(x)],
+            "cosh": [math.sinh],
+            "sin": [math.cos],
+            "sinh": [math.cosh],
+            "tan": [nan_if_exception(lambda x: 1 + math.tan(x) ** 2)],
+            "tanh": [nan_if_exception(lambda x: 1 - math.tanh(x) ** 2)],
+            "exp": [math.exp],
+            "exp2": [lambda y: _deriv_pow_1(2, y)],
+            "expm1": [math.exp],
+            "log10": [nan_if_exception(lambda x: 1 / x / math.log(10))],
+            "log1p": [nan_if_exception(lambda x: 1 / (1 + x))],
+            "degrees": [lambda x: math.degrees(1)],
+            "rad2deg": [lambda x: math.degrees(1)],
+            "radians": [lambda x: math.radians(1)],
+            "deg2rad": [lambda x: math.radians(1)],
+            "power": [_deriv_pow_0, _deriv_pow_1],
+            "sqrt": [nan_if_exception(lambda x: 0.5 / math.sqrt(x))],
+            "hypot": [
+                lambda x, y: x / math.hypot(x, y),
+                lambda x, y: y / math.hypot(x, y),
+            ],
+        }
+        # TODO: test arctan2, power, hypot
+        for func_str, derivatives in ufunc_derivatives.items():
+            implement_ufunc(func_str, derivatives)
+
+    @classmethod
+    def _add_numpy_comparative_ufuncs(cls):
+        def recursive_to_affine_scalar(arr):
+            if isinstance(arr, (list, tuple)):
+                return type(arr)([recursive_to_affine_scalar(i) for i in arr])
+            if isinstance(arr, np.ndarray):
+                return np.array([recursive_to_affine_scalar(i) for i in arr], "object")
+            return cls._to_affine_scalar(arr)
+
+        def implement_ufunc(func_str, func):
+            @implements(func_str, "ufunc")
+            def implementation(*inputs, **kwargs):
+                inputs = recursive_to_affine_scalar(inputs)
+                return apply_func_elementwise(func, inputs, kwargs, result_dtype=bool)
+
+            return implementation
+
+        ufunc_comparatives = {
+            "equal": ops.eq_on_aff_funcs,
+            "not_equal": ops.ne_on_aff_funcs,
+            "less": ops.lt_on_aff_funcs,
+            "less_equal": ops.le_on_aff_funcs,
+            "greater": ops.gt_on_aff_funcs,
+            "greater_equal": ops.ge_on_aff_funcs,
+        }
+        for func_str, func in ufunc_comparatives.items():
+            implement_ufunc(func_str, func)