diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 42e1048482743acf9ad5b8bb3addae9032dbb84d..a66afc54dd8e542de94010b66626795a30726b48 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -119,7 +119,7 @@ deploy2zenodo:
 
 pypi:
   stage: deploy
-  needs: 
+  needs:
     - job: "build"
       artifacts: true
   script:
diff --git a/pyproject.toml b/pyproject.toml
index 00e1d0509bcf8065954533285e315d25223d25ba..ac9197b047386272e0ee445d2648caeb57fbdb9a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,7 @@ dynamic = [
 description = "Utility package of the Quantum Technology group of RWTH Aachen"
 readme = "README.md"
 license = "GPL-3.0"
-requires-python = ">= 3.8"
+requires-python = ">= 3.9"
 authors = [
     { name = "Quantum Technology Group, RWTH Aachen University" },
 ]
@@ -20,10 +20,10 @@ classifiers = [
     "Intended Audience :: Science/Research",
     "License :: OSI Approved :: GNU General Public License v3 (GPLv3)",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
     "Topic :: Scientific/Engineering :: Physics",
     "Topic :: Utilities",
 ]
@@ -104,6 +104,7 @@ features = [
 ]
 dependencies = [
   "pytest",
+  "numpy>=2.0"
 ]
 
 [tool.hatch.envs.test.scripts]
diff --git a/qutil/functools.py b/qutil/functools.py
index e5568ebf8c700793418cc7a7be8a97ef95b99f62..cee0f80b6d9040ac611d1436fbd7259c330a979b 100644
--- a/qutil/functools.py
+++ b/qutil/functools.py
@@ -43,9 +43,16 @@ class FunctionChain(Generic[_C]):
     >>> x = np.arange(12).reshape(3, 4)
     >>> axis = 1
     >>> f_chain = FunctionChain(adder, multiplier, n_args=2)
-    >>> f_chain(x, axis)
-    (5016, -1)
+    >>> result, axis = f_chain(x, axis)
+    >>> print(result)
+    5016
     """
+    __array_interface__ = {
+        'shape': (),
+        'typestr': '|O',
+        'version': 3
+    }
+    """Describes to NumPy how to convert this object into an array."""
 
     def __init__(self, *functions: _C, n_args: int = 1, inspect_kwargs: bool = False):
         self.functions = functions
@@ -62,7 +69,8 @@ class FunctionChain(Generic[_C]):
         return len(self.functions)
 
     def __iter__(self) -> Iterator[_C]:
-        yield from self.functions
+        raise TypeError('The __iter__ method has been removed. Iterate over '
+                        'FunctionChain.functions instead')
 
     def __repr__(self) -> str:
         return (
diff --git a/qutil/itertools.py b/qutil/itertools.py
index 6b18ccd9a94da967606f70026d1ca46dd2f26aea..97d32cedad271c1981740fd7f45d8407a1e1007c 100644
--- a/qutil/itertools.py
+++ b/qutil/itertools.py
@@ -451,7 +451,7 @@ def absmin(iterable_or_value: SupportsAbs[_T] | Iterable[SupportsAbs[_T]],
     1
     >>> import numpy as np
     >>> x = np.random.default_rng().normal(20)
-    >>> absmin(x) == np.abs(x).min()
+    >>> print(absmin(x) == np.abs(x).min())
     True
     """
     iterable_or_value = (iterable_or_value, *others) if others else iterable_or_value
@@ -477,7 +477,7 @@ def absmax(iterable_or_value: SupportsAbs[_T] | Iterable[SupportsAbs[_T]],
     7
     >>> import numpy as np
     >>> x = np.random.default_rng().normal(20)
-    >>> absmax(x) == np.abs(x).max()
+    >>> print(absmax(x) == np.abs(x).max())
     True
     """
     iterable_or_value = (iterable_or_value, *others) if others else iterable_or_value
diff --git a/qutil/linalg.py b/qutil/linalg.py
index 6fcf8bcb245c1b3bcfcb3ce40a7be71a40d0b4f7..ef910ba3aca61dd534f7f98aa70740e461d76390 100644
--- a/qutil/linalg.py
+++ b/qutil/linalg.py
@@ -431,8 +431,11 @@ def check_phase_eq(psi: Sequence,
     >>> import qutil.qi
     >>> psi = qutil.qi.paulis[1]
     >>> phi = qutil.qi.paulis[1]*np.exp(1j*1.2345)
-    >>> check_phase_eq(psi, phi)
-    (True, 1.2345000000000603)
+    >>> eq, phase = check_phase_eq(psi, phi)
+    >>> print(eq)
+    True
+    >>> print(phase)
+    1.2345000000000603
     """
     if eps is None:
         # Tolerance the floating point eps times the # of flops for the matrix
@@ -499,9 +502,9 @@ def dot_HS(U: ndarray,
     --------
     >>> import qutil.qi
     >>> U, V = qutil.qi.paulis[1:3]
-    >>> dot_HS(U, V)
+    >>> print(dot_HS(U, V))
     0.0
-    >>> dot_HS(U, U)
+    >>> print(dot_HS(U, U))
     2
     """
     if eps is None:
diff --git a/qutil/pandas_tools.py b/qutil/pandas_tools.py
index 58414d369408ed69a200817d940998a48d3b20a5..2a1f71a7cfb1ebafe29a422ebfa1e8988a40bd2e 100644
--- a/qutil/pandas_tools.py
+++ b/qutil/pandas_tools.py
@@ -42,17 +42,17 @@ def consecutive_groupby(df: pd.DataFrame,
       4    a     4.0     43)]
 
     >>> list(df.groupby(['name', 'magic'])) # doctest: +NORMALIZE_WHITESPACE
-    [(('a', 42),
+    [(('a', np.int64(42)),
         name  number  magic
       0    a     0.0     42
       1    a     1.0     42),
-     (('a', 43),
+     (('a', np.int64(43)),
         name  number  magic
       4    a     4.0     43),
-     (('b', 42),
+     (('b', np.int64(42)),
         name  number  magic
       2    b     2.0     42),
-     (('b', 43),
+     (('b', np.int64(43)),
         name  number  magic
       3    b     3.0     43)]
 
diff --git a/qutil/signal_processing/fourier_space.py b/qutil/signal_processing/fourier_space.py
index 87810df3e460e7a729ed55d0ef0632a3b6333c17..cc330261b5703d4a53dbf9b965be7d17ae887314 100644
--- a/qutil/signal_processing/fourier_space.py
+++ b/qutil/signal_processing/fourier_space.py
@@ -95,7 +95,7 @@ def rms(x, f, /, out=None, *, axis: Optional[int] = None, where=True, dtype=None
     >>> x = 2*np.sqrt(2)*np.sin(2*np.pi*10*t)
     >>> xf = np.fft.fft(x)  # nb rfft would need to be scaled by factor √2
     >>> r, _ = rms(xf, ...)  # f not actually needed
-    >>> r  # doctest: +ELLIPSIS
+    >>> print(r)  # doctest: +ELLIPSIS
     1.9990007493755...
     >>> np.allclose(r, np.sqrt(x.mean()**2 + x.var()))
     True
diff --git a/qutil/signal_processing/real_space.py b/qutil/signal_processing/real_space.py
index 65e82e4d41a9920cc3bf663606b0100a5f201a90..2643122a51d9152c611fa93ef95200bb272fd774 100644
--- a/qutil/signal_processing/real_space.py
+++ b/qutil/signal_processing/real_space.py
@@ -70,7 +70,7 @@ def rms(x, /, out=None, *, axis: Optional[int] = None, where=True, dtype=None, k
     >>> t = np.linspace(0, 2*np.pi, 1001)
     >>> x = 2*np.sqrt(2)*np.sin(t)
     >>> r = rms(x)
-    >>> r  # doctest: +ELLIPSIS
+    >>> print(r)  # doctest: +ELLIPSIS
     1.9990007493755...
     >>> np.allclose(r, np.sqrt(x.mean()**2 + x.var()))
     True