diff --git a/INSTALL.txt b/INSTALL.txt
index 1aa398cc..3408197f 100644
--- a/INSTALL.txt
+++ b/INSTALL.txt
@@ -15,14 +15,7 @@ or, if additional access rights are needed (Unix):
 
   sudo python setup.py install
 
-* The tests programs (test_*.py) are meant to be run through the Nose
-testing framework.  This can be achieved for instance with a command
+* The tests programs (test_*.py) are meant to be run through pytest.  This can be achieved for instance with a command
 like
 
-  nosetests -sv uncertainties/
-
-or simply
-
-  nosetests uncertainties/
-
-(for a less verbose output).
+  pytest ./tests
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 6b8043cb..5798fcaa 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -57,3 +57,6 @@ Changelog = "https://github.com/lmfit/uncertainties/blob/master/CHANGES.rst"
 
 [project.optional-dependencies]
 optional = ["numpy"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
\ No newline at end of file
diff --git a/tests/helpers.py b/tests/helpers.py
new file mode 100644
index 00000000..eec53dff
--- /dev/null
+++ b/tests/helpers.py
@@ -0,0 +1,389 @@
+import random
+from math import isnan, isinf
+
+import uncertainties.core as uncert_core
+from uncertainties.core import ufloat, AffineScalarFunc, ufloat_fromstr
+
+def power_all_cases(op):
+    '''
+    Checks all cases for the value and derivatives of power-like
+    operator op (op is typically the built-in pow(), or math.pow()).
+
+    Checks only the details of special results like 0, 1 or NaN).
+
+    Different cases for the value of x**p and its derivatives are
+    tested by dividing the (x, p) plane with:
+
+    - x < 0, x = 0, x > 0
+    - p integer or not, p < 0, p = 0, p > 0
+
+    (not all combinations are distinct: for instance x > 0 gives
+    identical formulas for all p).
+    '''
+
+    zero = ufloat(0, 0.1)
+    zero2 = ufloat(0, 0.1)
+    one = ufloat(1, 0.1)
+    positive = ufloat(0.3, 0.01)
+    positive2 = ufloat(0.3, 0.01)
+    negative = ufloat(-0.3, 0.01)
+    integer = ufloat(-3, 0)
+    non_int_larger_than_one = ufloat(3.1, 0.01)
+    positive_smaller_than_one = ufloat(0.3, 0.01)
+
+    ## negative**integer
+
+    result = op(negative, integer)
+    assert not isnan(result.derivatives[negative])
+    assert isnan(result.derivatives[integer])
+
+    # Limit cases:
+    result = op(negative, one)
+    assert result.derivatives[negative] == 1
+    assert isnan(result.derivatives[one])
+
+    result = op(negative, zero)
+    assert result.derivatives[negative] == 0
+    assert isnan(result.derivatives[zero])
+
+    ## negative**non-integer
+
+    ## zero**...
+
+    result = op(zero, non_int_larger_than_one)
+    assert isnan(result.derivatives[zero])
+    assert result.derivatives[non_int_larger_than_one] == 0
+
+    # Special cases:
+    result = op(zero, one)
+    assert result.derivatives[zero] == 1
+    assert result.derivatives[one] == 0
+
+    result = op(zero, 2*one)
+    assert result.derivatives[zero] == 0
+    assert result.derivatives[one] == 0
+
+    result = op(zero, positive_smaller_than_one)
+    assert isnan(result.derivatives[zero])
+    assert result.derivatives[positive_smaller_than_one] == 0
+
+    result = op(zero, zero2)
+    assert result.derivatives[zero] == 0
+    assert isnan(result.derivatives[zero2])
+
+    ## positive**...: this is a quite regular case where the value and
+    ## the derivatives are all defined.
+
+    result = op(positive, positive2)
+    assert not isnan(result.derivatives[positive])
+    assert not isnan(result.derivatives[positive2])
+
+    result = op(positive, zero)
+    assert result.derivatives[positive] == 0
+    assert not isnan(result.derivatives[zero])
+
+    result = op(positive, negative)
+    assert not isnan(result.derivatives[positive])
+    assert not isnan(result.derivatives[negative])
+
+
+def power_special_cases(op):
+    '''
+    Checks special cases of the uncertainty power operator op (where
+    op is typically the built-in pow or uncertainties.umath.pow).
+
+    The values x = 0, x = 1 and x = NaN are special, as are null,
+    integral and NaN values of p.
+    '''
+
+    zero = ufloat(0, 0)
+    one = ufloat(1, 0)
+    p = ufloat(0.3, 0.01)
+
+    assert op(0, p) == 0
+    assert op(zero, p) == 0
+
+    # The outcome of 1**nan and nan**0 was undefined before Python
+    # 2.6 (http://docs.python.org/library/math.html#math.pow):
+    assert op(float('nan'), zero) == 1.0
+    assert op(one, float('nan')) == 1.0
+
+    # …**0 == 1.0:
+    assert op(p, 0) == 1.0
+    assert op(zero, 0) == 1.0
+    assert op((-p), 0) == 1.0
+    # …**zero:
+    assert op((-10.3), zero) == 1.0
+    assert op(0, zero) == 1.0
+    assert op(0.3, zero) == 1.0
+    assert op((-p), zero) == 1.0
+    assert op(zero, zero) == 1.0
+    assert op(p, zero) == 1.0
+
+    # one**… == 1.0
+    assert op(one, -3) == 1.0
+    assert op(one, -3.1) == 1.0
+    assert op(one, 0) == 1.0
+    assert op(one, 3) == 1.0
+    assert op(one, 3.1) == 1.0
+
+    # … with two numbers with uncertainties:
+    assert op(one, (-p)) == 1.0
+    assert op(one, zero) == 1.0
+    assert op(one, p) == 1.0
+    # 1**… == 1.0:
+    assert op(1., (-p)) == 1.0
+    assert op(1., zero) == 1.0
+    assert op(1., p) == 1.0
+
+def power_wrt_ref(op, ref_op):
+    '''
+    Checks special cases of the uncertainty power operator op (where
+    op is typically the built-in pow or uncertainties.umath.pow), by
+    comparing its results to the reference power operator ref_op
+    (which is typically the built-in pow or math.pow).
+    '''
+
+    # Negative numbers with uncertainty can be exponentiated to an
+    # integral power:
+    assert op(ufloat(-1.1, 0.1), -9).nominal_value == ref_op(-1.1, -9)
+
+    # Case of numbers with no uncertainty: should give the same result
+    # as numbers with uncertainties:
+    assert op(ufloat(-1, 0), 9) == ref_op(-1, 9)
+    assert op(ufloat(-1.1, 0), 9) == ref_op(-1.1, 9)
+
+
+###############################################################
+# TODO: move to uncertainties/testing.py
+###############################################################################
+
+# Utilities for unit testing
+
+def numbers_close(x, y, tolerance=1e-6):
+    """
+    Returns True if the given floats are close enough.
+
+    The given tolerance is the relative difference allowed, or the absolute
+    difference, if one of the numbers is 0.
+
+    NaN is allowed: it is considered close to itself.
+    """
+
+    # !!! Python 3.5+ has math.isclose(): maybe it could be used here.
+
+    # Instead of using a try and ZeroDivisionError, we do a test,
+    # NaN could appear silently:
+
+    if x != 0 and y != 0:
+        if isinf(x):
+            return isinf(y)
+        elif isnan(x):
+            return isnan(y)
+        else:
+            # Symmetric form of the test:
+            return 2*abs(x-y)/(abs(x)+abs(y)) < tolerance
+
+    else:  # Either x or y is zero
+        return abs(x or y) < tolerance
+
+def ufloats_close(x, y, tolerance=1e-6):
+    '''
+    Tests if two numbers with uncertainties are close, as random
+    variables: this is stronger than testing whether their nominal
+    value and standard deviation are close.
+
+    The tolerance is applied to both the nominal value and the
+    standard deviation of the difference between the numbers.
+    '''
+
+    diff = x-y
+    return (numbers_close(diff.nominal_value, 0, tolerance)
+            and numbers_close(diff.std_dev, 0, tolerance))
+
+class DerivativesDiffer(Exception):
+    pass
+
+
+def compare_derivatives(func, numerical_derivatives,
+                         num_args_list=None):
+    """
+    Checks the derivatives of a function 'func' (as returned by the
+    wrap() wrapper), by comparing them to the
+    'numerical_derivatives' functions.
+
+    Raises a DerivativesDiffer exception in case of problem.
+
+    These functions all take the number of arguments listed in
+    num_args_list.  If num_args is None, it is automatically obtained.
+
+    Tests are done on random arguments.
+    """
+
+    try:
+        funcname = func.name
+    except AttributeError:
+        funcname = func.__name__
+
+    # print "Testing", func.__name__
+
+    if not num_args_list:
+
+        # Detecting automatically the correct number of arguments is not
+        # always easy (because not all values are allowed, etc.):
+
+        num_args_table = {
+            'atanh': [1],
+            'log': [1, 2]  # Both numbers of arguments are tested
+            }
+        if funcname in num_args_table:
+            num_args_list = num_args_table[funcname]
+        else:
+
+            num_args_list = []
+
+            # We loop until we find reasonable function arguments:
+            # We get the number of arguments by trial and error:
+            for num_args in range(10):
+                try:
+                    #! Giving integer arguments is good for preventing
+                    # certain functions from failing even though num_args
+                    # is their correct number of arguments
+                    # (e.g. math.ldexp(x, i), where i must be an integer)
+                    func(*(1,)*num_args)
+                except TypeError:
+                    pass  # Not the right number of arguments
+                else:  # No error
+                    # num_args is a good number of arguments for func:
+                    num_args_list.append(num_args)
+
+            if not num_args_list:
+                raise Exception("Can't find a reasonable number of arguments"
+                                " for function '%s'." % funcname)
+
+    for num_args in num_args_list:
+
+        # Argument numbers that will have a random integer value:
+        integer_arg_nums = set()
+
+        if funcname == 'ldexp':
+            # The second argument must be an integer:
+            integer_arg_nums.add(1)
+
+        while True:
+            try:
+
+                # We include negative numbers, for more thorough tests:
+                args = []
+                for arg_num in range(num_args):
+                    if arg_num in integer_arg_nums:
+                        args.append(random.choice(range(-10, 10)))
+                    else:
+                        args.append(
+                            uncert_core.Variable(random.random()*4-2, 0))
+
+                # 'args', but as scalar values:
+                args_scalar = [uncert_core.nominal_value(v)
+                               for v in args]
+
+                func_approx = func(*args)
+
+                # Some functions yield simple Python constants, after
+                # wrapping in wrap(): no test has to be performed.
+                # Some functions also yield tuples...
+                if isinstance(func_approx, AffineScalarFunc):
+
+                    # We compare all derivatives:
+                    for (arg_num, (arg, numerical_deriv)) in (
+                        enumerate(zip(args, numerical_derivatives))):
+
+                        # Some arguments might not be differentiable:
+                        if isinstance(arg, int):
+                            continue
+
+                        fixed_deriv_value = func_approx.derivatives[arg]
+
+                        num_deriv_value = numerical_deriv(*args_scalar)
+
+                        # This message is useful: the user can see that
+                        # tests are really performed (instead of not being
+                        # performed, silently):
+                        print("Testing derivative #%d of %s at %s" % (
+                            arg_num, funcname, args_scalar))
+
+                        if not numbers_close(fixed_deriv_value,
+                                             num_deriv_value, 1e-4):
+
+                            # It is possible that the result is NaN:
+                            if not isnan(func_approx):
+                                raise DerivativesDiffer(
+                                    "Derivative #%d of function '%s' may be"
+                                    " wrong: at args = %s,"
+                                    " value obtained = %.16f,"
+                                    " while numerical approximation = %.16f."
+                                    % (arg_num, funcname, args,
+                                       fixed_deriv_value, num_deriv_value))
+
+            except ValueError as err:  # Arguments out of range, or of wrong type
+                # Factorial(real) lands here:
+                if str(err).startswith('factorial'):
+                    integer_arg_nums = set([0])
+                continue  # We try with different arguments
+            # Some arguments might have to be integers, for instance:
+            except TypeError as err:
+                if len(integer_arg_nums) == num_args:
+                    raise Exception("Incorrect testing procedure: unable to "
+                                    "find correct argument values for %s: %s"
+                                    % (funcname, err))
+
+                # Another argument might be forced to be an integer:
+                integer_arg_nums.add(random.choice(range(num_args)))
+            else:
+                # We have found reasonable arguments, and the test passed:
+                break
+
+###############################################################################
+
+
+try:
+    import numpy
+except ImportError:
+    pass
+else:
+
+    def uarrays_close(m1, m2, precision=1e-4):
+        """
+        Returns True iff m1 and m2 are almost equal, where elements
+        can be either floats or AffineScalarFunc objects.
+
+        Two independent AffineScalarFunc objects are deemed equal if
+        both their nominal value and uncertainty are equal (up to the
+        given precision).
+
+        m1, m2 -- NumPy arrays.
+
+        precision -- precision passed through to
+        uncertainties.test_uncertainties.numbers_close().
+        """
+
+        # ! numpy.allclose() is similar to this function, but does not
+        # work on arrays that contain numbers with uncertainties, because
+        # of the isinf() function.
+
+        for (elmt1, elmt2) in zip(m1.flat, m2.flat):
+
+            # For a simpler comparison, both elements are
+            # converted to AffineScalarFunc objects:
+            elmt1 = uncert_core.to_affine_scalar(elmt1)
+            elmt2 = uncert_core.to_affine_scalar(elmt2)
+
+            if not numbers_close(elmt1.nominal_value,
+                                 elmt2.nominal_value, precision):
+                return False
+
+            if not numbers_close(elmt1.std_dev,
+                                 elmt2.std_dev, precision):
+                return False
+        
+        return True
+
diff --git a/uncertainties/unumpy/test_ulinalg.py b/tests/test_ulinalg.py
similarity index 80%
rename from uncertainties/unumpy/test_ulinalg.py
rename to tests/test_ulinalg.py
index b72a0e11..c6c1de4a 100644
--- a/uncertainties/unumpy/test_ulinalg.py
+++ b/tests/test_ulinalg.py
@@ -1,16 +1,6 @@
-"""
-Tests for uncertainties.unumpy.ulinalg.
-
-These tests can be run through the Nose testing framework.
-
-(c) 2010-2016 by Eric O. LEBIGOT (EOL) <eric.lebigot@normalesup.org>.
-"""
-
 # Some tests are already performed in test_unumpy (unumpy contains a
 # matrix inversion, for instance).  They are not repeated here.
 
-from __future__ import division
-
 try:
     import numpy
 except ImportError:
@@ -18,7 +8,7 @@
     sys.exit()  # There is no reason to test the interface to NumPy
 
 from uncertainties import unumpy, ufloat
-from uncertainties.unumpy.test_unumpy import arrays_close
+from helpers import uarrays_close
 
 def test_list_inverse():
     "Test of the inversion of a square matrix"
@@ -56,7 +46,7 @@ def test_list_inverse():
     # Internal consistency: ulinalg.inv() must coincide with the
     # unumpy.matrix inverse, for square matrices (.I is the
     # pseudo-inverse, for non-square matrices, but inv() is not).
-    assert arrays_close(unumpy.ulinalg.inv(mat), mat.I)
+    assert uarrays_close(unumpy.ulinalg.inv(mat), mat.I)
 
 
 def test_list_pseudo_inverse():
@@ -68,8 +58,8 @@ def test_list_pseudo_inverse():
 
     # Internal consistency: the inverse and the pseudo-inverse yield
     # the same result on square matrices:
-    assert arrays_close(mat.I, unumpy.ulinalg.pinv(mat), 1e-4)
-    assert arrays_close(unumpy.ulinalg.inv(mat),
+    assert uarrays_close(mat.I, unumpy.ulinalg.pinv(mat), 1e-4)
+    assert uarrays_close(unumpy.ulinalg.inv(mat),
                           # Support for the optional pinv argument is
                           # tested:
                           unumpy.ulinalg.pinv(mat, 1e-15), 1e-4)
@@ -81,5 +71,5 @@ def test_list_pseudo_inverse():
     mat2 = unumpy.matrix([[x, y], [1, 3+x], [y, 2*x]])  # "Tall" matrix
 
     # Internal consistency:
-    assert arrays_close(mat1.I, unumpy.ulinalg.pinv(mat1, 1e-10))
-    assert arrays_close(mat2.I, unumpy.ulinalg.pinv(mat2, 1e-8))
+    assert uarrays_close(mat1.I, unumpy.ulinalg.pinv(mat1, 1e-10))
+    assert uarrays_close(mat2.I, unumpy.ulinalg.pinv(mat2, 1e-8))
diff --git a/uncertainties/test_umath.py b/tests/test_umath.py
similarity index 91%
rename from uncertainties/test_umath.py
rename to tests/test_umath.py
index c8555e7c..573c6b6a 100644
--- a/uncertainties/test_umath.py
+++ b/tests/test_umath.py
@@ -1,24 +1,12 @@
-"""
-Tests of the code in uncertainties.umath.
-
-These tests can be run through the Nose testing framework.
-
-(c) 2010-2016 by Eric O. LEBIGOT (EOL).
-"""
-
-from __future__ import division
-from __future__ import absolute_import
-
-# Standard modules
 import math
+from math import isnan, isinf
 
-# Local modules:
 from uncertainties import ufloat
 import uncertainties.core as uncert_core
 import uncertainties.umath_core as umath_core
 
-from . import test_uncertainties
-
+from helpers import (power_special_cases, power_all_cases, power_wrt_ref,
+    compare_derivatives, numbers_close)
 ###############################################################################
 # Unit tests
 
@@ -30,13 +18,12 @@ def test_fixed_derivatives_math_funcs():
     """
 
     for name in umath_core.many_scalars_to_scalar_funcs:
-        # print "Checking %s..." % name
         func = getattr(umath_core, name)
         # Numerical derivatives of func: the nominal value of func() results
         # is used as the underlying function:
         numerical_derivatives = uncert_core.NumericalDerivatives(
             lambda *args: func(*args))
-        test_uncertainties.compare_derivatives(func, numerical_derivatives)
+        compare_derivatives(func, numerical_derivatives)
 
     # Functions that are not in umath_core.many_scalars_to_scalar_funcs:
 
@@ -47,11 +34,11 @@ def frac_part_modf(x):
     def int_part_modf(x):
         return umath_core.modf(x)[1]
 
-    test_uncertainties.compare_derivatives(
+    compare_derivatives(
         frac_part_modf,
         uncert_core.NumericalDerivatives(
             lambda x: frac_part_modf(x)))
-    test_uncertainties.compare_derivatives(
+    compare_derivatives(
         int_part_modf,
         uncert_core.NumericalDerivatives(
             lambda x: int_part_modf(x)))
@@ -63,11 +50,11 @@ def mantissa_frexp(x):
     def exponent_frexp(x):
         return umath_core.frexp(x)[1]
 
-    test_uncertainties.compare_derivatives(
+    compare_derivatives(
         mantissa_frexp,
         uncert_core.NumericalDerivatives(
             lambda x: mantissa_frexp(x)))
-    test_uncertainties.compare_derivatives(
+    compare_derivatives(
         exponent_frexp,
         uncert_core.NumericalDerivatives(
             lambda x: exponent_frexp(x)))
@@ -171,7 +158,7 @@ def monte_carlo_calc(n_samples):
     # or assert_array_max_ulp. This is relevant for all vectorized
     # occurrences of numbers_close.
 
-    assert numpy.vectorize(test_uncertainties.numbers_close)(
+    assert numpy.vectorize(numbers_close)(
         covariances_this_module,
         covariances_samples,
         0.06).all(), (
@@ -182,7 +169,7 @@ def monte_carlo_calc(n_samples):
         )
 
     # The nominal values must be close:
-    assert test_uncertainties.numbers_close(
+    assert numbers_close(
         nominal_value_this_module,
         nominal_value_samples,
         # The scale of the comparison depends on the standard
@@ -278,14 +265,14 @@ def test_hypot():
     # Derivatives that cannot be calculated simply return NaN, with no
     # exception being raised, normally:
     result = umath_core.hypot(x, y)
-    assert test_uncertainties.isnan(result.derivatives[x])
-    assert test_uncertainties.isnan(result.derivatives[y])
+    assert isnan(result.derivatives[x])
+    assert isnan(result.derivatives[y])
 
 def test_power_all_cases():
     '''
     Test special cases of umath_core.pow().
     '''
-    test_uncertainties.power_all_cases(umath_core.pow)
+    power_all_cases(umath_core.pow)
 
 # test_power_special_cases() is similar to
 # test_uncertainties.py:test_power_special_cases(), but with small
@@ -296,7 +283,7 @@ def test_power_special_cases():
     Checks special cases of umath_core.pow().
     '''
 
-    test_uncertainties.power_special_cases(umath_core.pow)
+    power_special_cases(umath_core.pow)
 
     # We want the same behavior for numbers with uncertainties and for
     # math.pow() at their nominal values.
@@ -340,4 +327,4 @@ def test_power_wrt_ref():
     '''
     Checks special cases of the umath_core.pow() power operator.
     '''
-    test_uncertainties.power_wrt_ref(umath_core.pow, math.pow)
+    power_wrt_ref(umath_core.pow, math.pow)
\ No newline at end of file
diff --git a/uncertainties/test_uncertainties.py b/tests/test_uncertainties.py
similarity index 82%
rename from uncertainties/test_uncertainties.py
rename to tests/test_uncertainties.py
index 79c1db90..46078747 100644
--- a/uncertainties/test_uncertainties.py
+++ b/tests/test_uncertainties.py
@@ -1,229 +1,15 @@
-# coding=utf-8
 
-"""
-Tests of the code in uncertainties/__init__.py.
-
-These tests can be run through the Nose testing framework.
-
-(c) 2010-2016 by Eric O. LEBIGOT (EOL).
-"""
-
-from __future__ import division
-from __future__ import print_function
-
-
-# Standard modules
-from builtins import str
-from builtins import zip
-from builtins import map
-from builtins import range
 import copy
 import math
-from math import isnan, isinf
-import random
 import sys
-
-# 3rd-party modules
-# import nose.tools
-
-# Local modules
+from math import isnan
 
 import uncertainties.core as uncert_core
 from uncertainties.core import ufloat, AffineScalarFunc, ufloat_fromstr
 from uncertainties import umath
+from helpers import (power_special_cases, power_all_cases, power_wrt_ref,numbers_close,
+    ufloats_close, compare_derivatives, uarrays_close)
 
-# The following information is useful for making sure that the right
-# version of Python is running the tests (for instance with the Travis
-# Continuous Integration system):
-print("Testing with Python", sys.version)
-
-###############################################################################
-
-# Utilities for unit testing
-
-def numbers_close(x, y, tolerance=1e-6):
-    """
-    Returns True if the given floats are close enough.
-
-    The given tolerance is the relative difference allowed, or the absolute
-    difference, if one of the numbers is 0.
-
-    NaN is allowed: it is considered close to itself.
-    """
-
-    # !!! Python 3.5+ has math.isclose(): maybe it could be used here.
-
-    # Instead of using a try and ZeroDivisionError, we do a test,
-    # NaN could appear silently:
-
-    if x != 0 and y != 0:
-        if isinf(x):
-            return isinf(y)
-        elif isnan(x):
-            return isnan(y)
-        else:
-            # Symmetric form of the test:
-            return 2*abs(x-y)/(abs(x)+abs(y)) < tolerance
-
-    else:  # Either x or y is zero
-        return abs(x or y) < tolerance
-
-def ufloats_close(x, y, tolerance=1e-6):
-    '''
-    Tests if two numbers with uncertainties are close, as random
-    variables: this is stronger than testing whether their nominal
-    value and standard deviation are close.
-
-    The tolerance is applied to both the nominal value and the
-    standard deviation of the difference between the numbers.
-    '''
-
-    diff = x-y
-    return (numbers_close(diff.nominal_value, 0, tolerance)
-            and numbers_close(diff.std_dev, 0, tolerance))
-
-class DerivativesDiffer(Exception):
-    pass
-
-
-def compare_derivatives(func, numerical_derivatives,
-                         num_args_list=None):
-    """
-    Checks the derivatives of a function 'func' (as returned by the
-    wrap() wrapper), by comparing them to the
-    'numerical_derivatives' functions.
-
-    Raises a DerivativesDiffer exception in case of problem.
-
-    These functions all take the number of arguments listed in
-    num_args_list.  If num_args is None, it is automatically obtained.
-
-    Tests are done on random arguments.
-    """
-
-    try:
-        funcname = func.name
-    except AttributeError:
-        funcname = func.__name__
-
-    # print "Testing", func.__name__
-
-    if not num_args_list:
-
-        # Detecting automatically the correct number of arguments is not
-        # always easy (because not all values are allowed, etc.):
-
-        num_args_table = {
-            'atanh': [1],
-            'log': [1, 2]  # Both numbers of arguments are tested
-            }
-        if funcname in num_args_table:
-            num_args_list = num_args_table[funcname]
-        else:
-
-            num_args_list = []
-
-            # We loop until we find reasonable function arguments:
-            # We get the number of arguments by trial and error:
-            for num_args in range(10):
-                try:
-                    #! Giving integer arguments is good for preventing
-                    # certain functions from failing even though num_args
-                    # is their correct number of arguments
-                    # (e.g. math.ldexp(x, i), where i must be an integer)
-                    func(*(1,)*num_args)
-                except TypeError:
-                    pass  # Not the right number of arguments
-                else:  # No error
-                    # num_args is a good number of arguments for func:
-                    num_args_list.append(num_args)
-
-            if not num_args_list:
-                raise Exception("Can't find a reasonable number of arguments"
-                                " for function '%s'." % funcname)
-
-    for num_args in num_args_list:
-
-        # Argument numbers that will have a random integer value:
-        integer_arg_nums = set()
-
-        if funcname == 'ldexp':
-            # The second argument must be an integer:
-            integer_arg_nums.add(1)
-
-        while True:
-            try:
-
-                # We include negative numbers, for more thorough tests:
-                args = []
-                for arg_num in range(num_args):
-                    if arg_num in integer_arg_nums:
-                        args.append(random.choice(range(-10, 10)))
-                    else:
-                        args.append(
-                            uncert_core.Variable(random.random()*4-2, 0))
-
-                # 'args', but as scalar values:
-                args_scalar = [uncert_core.nominal_value(v)
-                               for v in args]
-
-                func_approx = func(*args)
-
-                # Some functions yield simple Python constants, after
-                # wrapping in wrap(): no test has to be performed.
-                # Some functions also yield tuples...
-                if isinstance(func_approx, AffineScalarFunc):
-
-                    # We compare all derivatives:
-                    for (arg_num, (arg, numerical_deriv)) in (
-                        enumerate(zip(args, numerical_derivatives))):
-
-                        # Some arguments might not be differentiable:
-                        if isinstance(arg, int):
-                            continue
-
-                        fixed_deriv_value = func_approx.derivatives[arg]
-
-                        num_deriv_value = numerical_deriv(*args_scalar)
-
-                        # This message is useful: the user can see that
-                        # tests are really performed (instead of not being
-                        # performed, silently):
-                        print("Testing derivative #%d of %s at %s" % (
-                            arg_num, funcname, args_scalar))
-
-                        if not numbers_close(fixed_deriv_value,
-                                             num_deriv_value, 1e-4):
-
-                            # It is possible that the result is NaN:
-                            if not isnan(func_approx):
-                                raise DerivativesDiffer(
-                                    "Derivative #%d of function '%s' may be"
-                                    " wrong: at args = %s,"
-                                    " value obtained = %.16f,"
-                                    " while numerical approximation = %.16f."
-                                    % (arg_num, funcname, args,
-                                       fixed_deriv_value, num_deriv_value))
-
-            except ValueError as err:  # Arguments out of range, or of wrong type
-                # Factorial(real) lands here:
-                if str(err).startswith('factorial'):
-                    integer_arg_nums = set([0])
-                continue  # We try with different arguments
-            # Some arguments might have to be integers, for instance:
-            except TypeError as err:
-                if len(integer_arg_nums) == num_args:
-                    raise Exception("Incorrect testing procedure: unable to "
-                                    "find correct argument values for %s: %s"
-                                    % (funcname, err))
-
-                # Another argument might be forced to be an integer:
-                integer_arg_nums.add(random.choice(range(num_args)))
-            else:
-                # We have found reasonable arguments, and the test passed:
-                break
-
-###############################################################################
 
 def test_value_construction():
     '''
@@ -383,8 +169,6 @@ def check_op(op, num_args):
     for op in uncert_core.modified_ops_with_reflection:
         check_op(op, 2)
 
-# Additional, more complex checks, for use with the nose unit testing
-# framework.
 
 def test_copy():
     "Standard copy module integration"
@@ -1236,88 +1020,6 @@ def test_power_all_cases():
 
     power_all_cases(pow)
 
-def power_all_cases(op):
-    '''
-    Checks all cases for the value and derivatives of power-like
-    operator op (op is typically the built-in pow(), or math.pow()).
-
-    Checks only the details of special results like 0, 1 or NaN).
-
-    Different cases for the value of x**p and its derivatives are
-    tested by dividing the (x, p) plane with:
-
-    - x < 0, x = 0, x > 0
-    - p integer or not, p < 0, p = 0, p > 0
-
-    (not all combinations are distinct: for instance x > 0 gives
-    identical formulas for all p).
-    '''
-
-    zero = ufloat(0, 0.1)
-    zero2 = ufloat(0, 0.1)
-    one = ufloat(1, 0.1)
-    positive = ufloat(0.3, 0.01)
-    positive2 = ufloat(0.3, 0.01)
-    negative = ufloat(-0.3, 0.01)
-    integer = ufloat(-3, 0)
-    non_int_larger_than_one = ufloat(3.1, 0.01)
-    positive_smaller_than_one = ufloat(0.3, 0.01)
-
-    ## negative**integer
-
-    result = op(negative, integer)
-    assert not isnan(result.derivatives[negative])
-    assert isnan(result.derivatives[integer])
-
-    # Limit cases:
-    result = op(negative, one)
-    assert result.derivatives[negative] == 1
-    assert isnan(result.derivatives[one])
-
-    result = op(negative, zero)
-    assert result.derivatives[negative] == 0
-    assert isnan(result.derivatives[zero])
-
-    ## negative**non-integer
-
-    ## zero**...
-
-    result = op(zero, non_int_larger_than_one)
-    assert isnan(result.derivatives[zero])
-    assert result.derivatives[non_int_larger_than_one] == 0
-
-    # Special cases:
-    result = op(zero, one)
-    assert result.derivatives[zero] == 1
-    assert result.derivatives[one] == 0
-
-    result = op(zero, 2*one)
-    assert result.derivatives[zero] == 0
-    assert result.derivatives[one] == 0
-
-    result = op(zero, positive_smaller_than_one)
-    assert isnan(result.derivatives[zero])
-    assert result.derivatives[positive_smaller_than_one] == 0
-
-    result = op(zero, zero2)
-    assert result.derivatives[zero] == 0
-    assert isnan(result.derivatives[zero2])
-
-    ## positive**...: this is a quite regular case where the value and
-    ## the derivatives are all defined.
-
-    result = op(positive, positive2)
-    assert not isnan(result.derivatives[positive])
-    assert not isnan(result.derivatives[positive2])
-
-    result = op(positive, zero)
-    assert result.derivatives[positive] == 0
-    assert not isnan(result.derivatives[zero])
-
-    result = op(positive, negative)
-    assert not isnan(result.derivatives[positive])
-    assert not isnan(result.derivatives[negative])
-
 
 ###############################################################################
 
@@ -1361,79 +1063,12 @@ def test_power_special_cases():
     else:
         raise Exception('A proper exception should have been raised')
 
-def power_special_cases(op):
-    '''
-    Checks special cases of the uncertainty power operator op (where
-    op is typically the built-in pow or uncertainties.umath.pow).
-
-    The values x = 0, x = 1 and x = NaN are special, as are null,
-    integral and NaN values of p.
-    '''
-
-    zero = ufloat(0, 0)
-    one = ufloat(1, 0)
-    p = ufloat(0.3, 0.01)
-
-    assert op(0, p) == 0
-    assert op(zero, p) == 0
-
-    # The outcome of 1**nan and nan**0 was undefined before Python
-    # 2.6 (http://docs.python.org/library/math.html#math.pow):
-    assert op(float('nan'), zero) == 1.0
-    assert op(one, float('nan')) == 1.0
-
-    # …**0 == 1.0:
-    assert op(p, 0) == 1.0
-    assert op(zero, 0) == 1.0
-    assert op((-p), 0) == 1.0
-    # …**zero:
-    assert op((-10.3), zero) == 1.0
-    assert op(0, zero) == 1.0
-    assert op(0.3, zero) == 1.0
-    assert op((-p), zero) == 1.0
-    assert op(zero, zero) == 1.0
-    assert op(p, zero) == 1.0
-
-    # one**… == 1.0
-    assert op(one, -3) == 1.0
-    assert op(one, -3.1) == 1.0
-    assert op(one, 0) == 1.0
-    assert op(one, 3) == 1.0
-    assert op(one, 3.1) == 1.0
-
-    # … with two numbers with uncertainties:
-    assert op(one, (-p)) == 1.0
-    assert op(one, zero) == 1.0
-    assert op(one, p) == 1.0
-    # 1**… == 1.0:
-    assert op(1., (-p)) == 1.0
-    assert op(1., zero) == 1.0
-    assert op(1., p) == 1.0
-
-
 def test_power_wrt_ref():
     '''
     Checks special cases of the built-in pow() power operator.
     '''
     power_wrt_ref(pow, pow)
 
-def power_wrt_ref(op, ref_op):
-    '''
-    Checks special cases of the uncertainty power operator op (where
-    op is typically the built-in pow or uncertainties.umath.pow), by
-    comparing its results to the reference power operator ref_op
-    (which is typically the built-in pow or math.pow).
-    '''
-
-    # Negative numbers with uncertainty can be exponentiated to an
-    # integral power:
-    assert op(ufloat(-1.1, 0.1), -9).nominal_value == ref_op(-1.1, -9)
-
-    # Case of numbers with no uncertainty: should give the same result
-    # as numbers with uncertainties:
-    assert op(ufloat(-1, 0), 9) == ref_op(-1, 9)
-    assert op(ufloat(-1.1, 0), 9) == ref_op(-1.1, 9)
-
 
 ###############################################################################
 
@@ -2096,43 +1731,6 @@ def test_custom_pretty_print_and_latex():
     pass
 else:
 
-    def arrays_close(m1, m2, precision=1e-4):
-        """
-        Returns True iff m1 and m2 are almost equal, where elements
-        can be either floats or AffineScalarFunc objects.
-
-        Two independent AffineScalarFunc objects are deemed equal if
-        both their nominal value and uncertainty are equal (up to the
-        given precision).
-
-        m1, m2 -- NumPy arrays.
-
-        precision -- precision passed through to
-        uncertainties.test_uncertainties.numbers_close().
-        """
-
-        # ! numpy.allclose() is similar to this function, but does not
-        # work on arrays that contain numbers with uncertainties, because
-        # of the isinf() function.
-
-        for (elmt1, elmt2) in zip(m1.flat, m2.flat):
-
-            # For a simpler comparison, both elements are
-            # converted to AffineScalarFunc objects:
-            elmt1 = uncert_core.to_affine_scalar(elmt1)
-            elmt2 = uncert_core.to_affine_scalar(elmt2)
-
-            if not numbers_close(elmt1.nominal_value,
-                                 elmt2.nominal_value, precision):
-                return False
-
-            if not numbers_close(elmt1.std_dev,
-                                 elmt2.std_dev, precision):
-                return False
-        
-        return True
-
-
     def test_numpy_comparison():
         "Comparison with a NumPy array."
 
@@ -2198,7 +1796,7 @@ def test_correlated_values():
         covs = uncert_core.covariance_matrix([x, y, z])
 
         # Test of the diagonal covariance elements:
-        assert arrays_close(
+        assert uarrays_close(
             numpy.array([v.std_dev**2 for v in (x, y, z)]),
             numpy.array(covs).diagonal())
 
@@ -2210,15 +1808,15 @@ def test_correlated_values():
             tags = ['x', 'y', 'z'])
 
         # Even the uncertainties should be correctly reconstructed:
-        assert arrays_close(numpy.array((x, y, z)),
+        assert uarrays_close(numpy.array((x, y, z)),
                               numpy.array((x_new, y_new, z_new)))
 
         # ... and the covariances too:
-        assert arrays_close(
+        assert uarrays_close(
             numpy.array(covs),
             numpy.array(uncert_core.covariance_matrix([x_new, y_new, z_new])))
 
-        assert arrays_close(
+        assert uarrays_close(
             numpy.array([z_new]), numpy.array([-3*x_new+y_new]))
 
         ####################
@@ -2238,12 +1836,12 @@ def test_correlated_values():
             [x.nominal_value for x in [u, v, sum_value]],
             cov_matrix)
 
-        # arrays_close() is used instead of numbers_close() because
+        # uarrays_close() is used instead of numbers_close() because
         # it compares uncertainties too:
-        assert arrays_close(numpy.array([u]), numpy.array([u2]))
-        assert arrays_close(numpy.array([v]), numpy.array([v2]))
-        assert arrays_close(numpy.array([sum_value]), numpy.array([sum2]))
-        assert arrays_close(numpy.array([0]),
+        assert uarrays_close(numpy.array([u]), numpy.array([u2]))
+        assert uarrays_close(numpy.array([v]), numpy.array([v2]))
+        assert uarrays_close(numpy.array([sum_value]), numpy.array([sum2]))
+        assert uarrays_close(numpy.array([0]),
                             numpy.array([sum2-(u2+2*v2)]))
 
 
@@ -2287,7 +1885,7 @@ def test_correlated_values():
             assert numbers_close(variable.n, nom_value)
             assert numbers_close(variable.s**2, variance) 
         
-        assert arrays_close(
+        assert uarrays_close(
             cov,
             numpy.array(uncert_core.covariance_matrix(variables)))
 
@@ -2322,18 +1920,18 @@ def test_correlated_values_correlation_mat():
         x2, y2, z2 = uncert_core.correlated_values_norm(
             list(zip(nominal_values, std_devs)), corr_mat)
 
-        # arrays_close() is used instead of numbers_close() because
+        # uarrays_close() is used instead of numbers_close() because
         # it compares uncertainties too:
 
         # Test of individual variables:
-        assert arrays_close(numpy.array([x]), numpy.array([x2]))
-        assert arrays_close(numpy.array([y]), numpy.array([y2]))
-        assert arrays_close(numpy.array([z]), numpy.array([z2]))
+        assert uarrays_close(numpy.array([x]), numpy.array([x2]))
+        assert uarrays_close(numpy.array([y]), numpy.array([y2]))
+        assert uarrays_close(numpy.array([z]), numpy.array([z2]))
 
         # Partial correlation test:
-        assert arrays_close(numpy.array([0]), numpy.array([z2-(-3*x2+y2)]))
+        assert uarrays_close(numpy.array([0]), numpy.array([z2-(-3*x2+y2)]))
 
         # Test of the full covariance matrix:
-        assert arrays_close(
+        assert uarrays_close(
             numpy.array(cov_mat),
             numpy.array(uncert_core.covariance_matrix([x2, y2, z2])))
diff --git a/uncertainties/unumpy/test_unumpy.py b/tests/test_unumpy.py
similarity index 93%
rename from uncertainties/unumpy/test_unumpy.py
rename to tests/test_unumpy.py
index 5e2147a5..4958bc54 100644
--- a/uncertainties/unumpy/test_unumpy.py
+++ b/tests/test_unumpy.py
@@ -1,26 +1,14 @@
-"""
-Tests of the code in uncertainties/unumpy/__init__.py.
-
-These tests can be run through the Nose testing framework.
-
-(c) 2010-2016 by Eric O. LEBIGOT (EOL).
-"""
-
-from __future__ import division
-
-# 3rd-party modules:
 try:
     import numpy
 except ImportError:
     import sys
     sys.exit()  # There is no reason to test the interface to NumPy
 
-# Local modules:
 import uncertainties
 import uncertainties.core as uncert_core
-from uncertainties import ufloat, unumpy, test_uncertainties
+from uncertainties import ufloat, unumpy
 from uncertainties.unumpy import core
-from uncertainties.test_uncertainties import numbers_close, arrays_close
+from helpers import numbers_close, uarrays_close
 
 def test_numpy():
 
@@ -161,7 +149,7 @@ def test_inverse():
     assert numbers_close(m_double_inverse[0, 0].std_dev,
                           m[0, 0].std_dev)
 
-    assert arrays_close(m_double_inverse, m)
+    assert uarrays_close(m_double_inverse, m)
 
     # Partial test:
     assert derivatives_close(m_double_inverse[0, 0], m[0, 0])
@@ -178,7 +166,7 @@ def test_inverse():
 
     # Correlations between m and m_inverse should create a perfect
     # inversion:
-    assert arrays_close(m * m_inverse,  numpy.eye(m.shape[0]))
+    assert uarrays_close(m * m_inverse,  numpy.eye(m.shape[0]))
 
 def test_wrap_array_func():
     '''
@@ -213,7 +201,7 @@ def f(mat, *args, **kwargs):
     m_f_wrapped = f_wrapped(m, 2, factor=10)
     m_f_unc = f_unc(m, 2, factor=10)
 
-    assert arrays_close(m_f_wrapped, m_f_unc)
+    assert uarrays_close(m_f_wrapped, m_f_unc)
 
 
 def test_pseudo_inverse():
@@ -233,7 +221,7 @@ def test_pseudo_inverse():
     rcond = 1e-8  # Test of the second argument to pinv()
     m_pinv_num = pinv_num(m, rcond)
     m_pinv_package = core.pinv(m, rcond)
-    assert arrays_close(m_pinv_num, m_pinv_package)
+    assert uarrays_close(m_pinv_num, m_pinv_package)
 
     ##########
     # Example with a non-full rank rectangular matrix:
@@ -241,14 +229,14 @@ def test_pseudo_inverse():
     m = unumpy.matrix([vector, vector])
     m_pinv_num = pinv_num(m, rcond)
     m_pinv_package = core.pinv(m, rcond)
-    assert arrays_close(m_pinv_num, m_pinv_package)
+    assert uarrays_close(m_pinv_num, m_pinv_package)
 
     ##########
     # Example with a non-full-rank square matrix:
     m = unumpy.matrix([[ufloat(10, 1), 0], [3, 0]])
     m_pinv_num = pinv_num(m, rcond)
     m_pinv_package = core.pinv(m, rcond)
-    assert arrays_close(m_pinv_num, m_pinv_package)
+    assert uarrays_close(m_pinv_num, m_pinv_package)
 
 def test_broadcast_funcs():
     """