diff --git a/.git_archival.txt b/.git_archival.txt
new file mode 100644
index 0000000..7c51009
--- /dev/null
+++ b/.git_archival.txt
@@ -0,0 +1,3 @@
+node: $Format:%H$
+node-date: $Format:%cI$
+describe-name: $Format:%(describe:tags=true,match=*[0-9]*)$
diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..7d3026f
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,2 @@
+.git_archival.txt  export-subst
+*.ipynb  linguist-documentation
diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md
new file mode 100644
index 0000000..910ad8e
--- /dev/null
+++ b/.github/CONTRIBUTING.md
@@ -0,0 +1,66 @@
+See the [Scikit-HEP Developer introduction][skhep-dev-intro] for a detailed
+description of best practices for developing Scikit-HEP packages.
+
+[skhep-dev-intro]: https://scikit-hep.org/developer/intro
+
+# Contributing
+
+## Setting up a development environment
+
+### Nox
+
+The fastest way to start with development is to use nox. If you don't have nox,
+you can use `uvx nox` to run it without installing, or `uv tool install nox`. If
+you don't have uv (or pipx, which is similar), then you can install it with
+`pip install uv` (the only case where installing an application with regular pip
+is reasonable). If you use macOS, then uv, prek, and nox are all in brew, use
+`brew install uv prek nox`.
+
+To use, run `nox`. This will lint and test using the Python version nox runs on.
+You can also run specific jobs:
+
+```console
+$ nox -l # List all the defined sessions
+$ nox -s lint  # Lint only
+$ nox -s tests  # Python tests
+$ nox -s tests -P 3.12  # Python tests on a specific Python version
+$ nox -s minimums  # Python tests with the lowest supported dependency versions
+$ nox -s cov  # Python tests with an HTML coverage report
+$ nox -s bench  # Benchmarks
+$ nox -s docs  # Build and serve the docs
+$ nox -s plots  # Regenerate the figures used in the docs
+$ nox -s build  # Make an SDist and wheel
+```
+
+Nox handles everything for you, including setting up a temporary virtual
+environment for each run.
+
+### uv
+
+For extended development, you can set up a development environment with uv:
+
+```bash
+uv sync --group test --group docs
+uv run pytest
+```
+
+## Post setup
+
+You should prepare prek, which will help you by checking that commits pass
+required checks:
+
+```bash
+uv tool install prek # or brew install prek on macOS
+prek install # Will install a pre-commit hook into the git repo
+```
+
+You can also/alternatively run `prek` (changes only) or `prek -a` to check even
+without installing the hook.
+
+## Testing
+
+Use pytest to run the unit checks:
+
+```bash
+uv run pytest
+```
diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 0000000..e39daf9
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,13 @@
+version: 2
+updates:
+  # Maintain dependencies for GitHub Actions
+  - package-ecosystem: "github-actions"
+    directory: "/"
+    schedule:
+      interval: "monthly"
+    cooldown:
+      default-days: 7
+    groups:
+      actions:
+        patterns:
+          - "*"
diff --git a/.github/release.yml b/.github/release.yml
new file mode 100644
index 0000000..5f89818
--- /dev/null
+++ b/.github/release.yml
@@ -0,0 +1,5 @@
+changelog:
+  exclude:
+    authors:
+      - dependabot[bot]
+      - pre-commit-ci[bot]
diff --git a/.github/workflows/cd.yml b/.github/workflows/cd.yml
new file mode 100644
index 0000000..bd79009
--- /dev/null
+++ b/.github/workflows/cd.yml
@@ -0,0 +1,63 @@
+name: CD
+
+on:
+  workflow_dispatch:
+  pull_request:
+  push:
+    branches:
+      - main
+  release:
+    types:
+      - published
+
+permissions: {}
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+env:
+  # Many color libraries just need this to be set to any value, but at least
+  # one distinguishes color depth, where "3" -> "256-bit color".
+  FORCE_COLOR: 3
+
+jobs:
+  dist:
+    name: Distribution build
+    runs-on: ubuntu-latest
+    permissions:
+      contents: read # to check out the repository
+
+    steps:
+      - uses: actions/checkout@v7
+        with:
+          persist-credentials: false
+          fetch-depth: 0
+
+      - uses: hynek/build-and-inspect-python-package@v3.0.1
+
+  publish:
+    needs: [dist]
+    name: Publish to PyPI
+    environment:
+      name: pypi
+      url: https://pypi.org/p/jacobi
+    permissions:
+      id-token: write # for trusted publishing to PyPI
+      attestations: write # to generate artifact attestations
+      contents: read # to check out the repository
+    runs-on: ubuntu-latest
+    if: github.event_name == 'release' && github.event.action == 'published'
+
+    steps:
+      - uses: actions/download-artifact@v8
+        with:
+          name: Packages
+          path: dist
+
+      - name: Generate artifact attestation for sdist and wheel
+        uses: actions/attest@v4
+        with:
+          subject-path: "dist/jacobi-*"
+
+      - uses: pypa/gh-action-pypi-publish@release/v1
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..02abf4f
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,117 @@
+name: CI
+
+on:
+  workflow_dispatch:
+  pull_request:
+  push:
+    branches:
+      - main
+
+permissions: {}
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+env:
+  # Many color libraries just need this variable to be set to any value.
+  # Set it to 3 to support 8-bit color graphics (256 colors per channel)
+  # for libraries that care about the value set.
+  FORCE_COLOR: 3
+
+jobs:
+  lint:
+    name: Lint
+    runs-on: ubuntu-latest
+    permissions:
+      contents: read # to check out the repository
+    steps:
+      - uses: actions/checkout@v7
+        with:
+          persist-credentials: false
+          fetch-depth: 0
+
+      - uses: actions/setup-python@v7
+        with:
+          python-version: "3.x"
+
+      - uses: astral-sh/setup-uv@v10.1.0
+
+      - uses: j178/prek-action@v3.0.0
+
+      - name: Run Pylint
+        run: uvx nox -s pylint -- --output-format=github
+
+  checks:
+    name: Check Python ${{ matrix.python-version }} on ${{ matrix.runs-on }}
+    runs-on: ${{ matrix.runs-on }}
+    permissions:
+      contents: read # to check out the repository
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
+        runs-on: [ubuntu-latest]
+
+        include:
+          - python-version: "3.10"
+            runs-on: windows-latest
+          - python-version: "3.14"
+            runs-on: windows-latest
+          - python-version: "3.10"
+            runs-on: macos-latest
+          - python-version: "3.14"
+            runs-on: macos-latest
+          - python-version: "3.14t"
+            runs-on: ubuntu-latest
+          - python-version: "pypy-3.11"
+            runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v7
+        with:
+          persist-credentials: false
+          fetch-depth: 0
+
+      - uses: actions/setup-python@v7
+        id: setup-python
+        with:
+          python-version: ${{ matrix.python-version }}
+          allow-prereleases: true
+
+      - uses: astral-sh/setup-uv@v10.1.0
+
+      - name: Install package
+        run: uv sync --group test
+        env:
+          UV_PYTHON: ${{ steps.setup-python.outputs.python-path }}
+
+      - name: Test package
+        run: uv run pytest -ra --cov --cov-report=term --durations=20
+
+  minimums:
+    name: Check minimums
+    runs-on: ubuntu-latest
+    permissions:
+      contents: read # to check out the repository
+    steps:
+      - uses: actions/checkout@v7
+        with:
+          persist-credentials: false
+          fetch-depth: 0
+
+      - uses: astral-sh/setup-uv@v10.1.0
+
+      - name: Test with the lowest supported dependency versions
+        run: uvx nox -s minimums
+
+  pass:
+    if: always()
+    needs: [lint, checks, minimums]
+    runs-on: ubuntu-slim
+    timeout-minutes: 2
+    steps:
+      - name: Decide whether the needed jobs succeeded or failed
+        uses: re-actors/alls-green@release/v1
+        with:
+          jobs: ${{ toJSON(needs) }}
diff --git a/.github/workflows/doc.yml b/.github/workflows/doc.yml
deleted file mode 100644
index 6106331..0000000
--- a/.github/workflows/doc.yml
+++ /dev/null
@@ -1,50 +0,0 @@
-name: doc
-
-on:
-  push:
-    tags:
-      - '**'
-  pull_request:
-
-concurrency:
-  group: ${{ github.workflow }}-${{ github.head_ref }}
-  cancel-in-progress: true
-
-jobs:
-  build:
-    runs-on: ubuntu-latest
-    steps:
-    - uses: actions/checkout@v3
-      with:
-        fetch-depth: 0  # needed for setuptools_scm
-    - uses: actions/setup-python@v3
-      with:
-        python-version: "3.9"
-
-    - run: sudo apt-get install pandoc
-    - run: python -m pip install --prefer-binary -v -e .[doc]
-    - run: python -m ipykernel install --user --name python3
-    - run: python doc/build.py
-
-    - uses: actions/upload-pages-artifact@v2
-      with:
-        path: doc/_build/html
-
-  deploy:
-    if: contains(github.event.ref, '/tags/')
-    needs: build
-    # Set permissions of the GITHUB_TOKEN to allow deployment to GitHub Pages
-    permissions:
-      contents: read
-      pages: write
-      id-token: write
-
-    environment:
-      name: github-pages
-      url: ${{ steps.deployment.outputs.page_url }}
-
-    runs-on: ubuntu-latest
-
-    steps:
-    - uses: actions/configure-pages@v3
-    - uses: actions/deploy-pages@v2
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
new file mode 100644
index 0000000..eca1b0c
--- /dev/null
+++ b/.github/workflows/docs.yml
@@ -0,0 +1,63 @@
+name: Docs
+
+on:
+  workflow_dispatch:
+  pull_request:
+  push:
+    branches:
+      - main
+  release:
+    types:
+      - published
+
+permissions: {}
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+env:
+  # Many color libraries just need this to be set to any value, but at least
+  # one distinguishes color depth, where "3" -> "256-bit color".
+  FORCE_COLOR: 3
+
+jobs:
+  build:
+    name: Build docs
+    runs-on: ubuntu-latest
+    permissions:
+      contents: read # to check out the repository
+    steps:
+      - uses: actions/checkout@v7
+        with:
+          persist-credentials: false
+          fetch-depth: 0
+
+      - uses: astral-sh/setup-uv@v10.1.0
+
+      - name: Build the documentation
+        run: uvx nox -s docs
+
+      - uses: actions/upload-pages-artifact@v5
+        with:
+          path: docs/_build/html
+
+  deploy:
+    name: Deploy docs to GitHub Pages
+    if: |
+      github.event_name == 'workflow_dispatch' ||
+      (github.event_name == 'release' && github.event.action == 'published')
+    needs: [build]
+    permissions:
+      contents: read # to read the repository settings
+      pages: write # to deploy to GitHub Pages
+      id-token: write # to verify the deployment originates from this workflow
+    environment:
+      name: github-pages
+      url: ${{ steps.deployment.outputs.page_url }}
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/configure-pages@v6
+
+      - id: deployment
+        uses: actions/deploy-pages@v5
diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
deleted file mode 100644
index 88a37e3..0000000
--- a/.github/workflows/publish.yml
+++ /dev/null
@@ -1,34 +0,0 @@
-name: publish
-
-on:
-  push:
-    tags:
-      - '**'
-  workflow_dispatch:
-
-jobs:
-  deploy:
-
-    runs-on: ubuntu-latest
-
-    steps:
-    - uses: actions/checkout@v3
-
-    - uses: actions/setup-python@v3
-      with:
-        python-version: '3.9'
-
-    - run: |
-        python -m pip install --upgrade pip build
-        python -m build .
-
-    - run: |
-        ls dist
-        python -m pip install $(echo dist/jacobi-*.whl)'[test]'
-        python -m pytest
-
-    - uses: pypa/gh-action-pypi-publish@release/v1
-      if: contains(github.event.ref, '/tags/')
-      with:
-        user: __token__
-        password: ${{secrets.PYPI_TOKEN}}
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
deleted file mode 100644
index 33236be..0000000
--- a/.github/workflows/test.yml
+++ /dev/null
@@ -1,34 +0,0 @@
-name: test
-
-on:
-  pull_request:
-    paths-ignore:
-      - 'doc/**'
-
-concurrency:
-  group: ${{ github.workflow }}-${{ github.head_ref }}
-  cancel-in-progress: true
-
-jobs:
-  test:
-    runs-on: ${{ matrix.os }}
-
-    strategy:
-      matrix:
-        include:
-          # version number must be string, otherwise 3.10 becomes 3.1
-          - os: windows-latest
-            python-version: "3.11"
-          - os: macos-latest
-            python-version: "3.8"
-          - os: ubuntu-latest
-            python-version: "pypy-3.8"
-      fail-fast: false
-    steps:
-    - uses: actions/checkout@v3
-    - uses: actions/setup-python@v3
-      with:
-        python-version: ${{ matrix.python-version }}
-    - run: python -m pip install --upgrade pip
-    - run: python -m pip install --prefer-binary -v -e '.[test]'
-    - run: python -m pytest
diff --git a/.github/zizmor.yml b/.github/zizmor.yml
new file mode 100644
index 0000000..ab94a2a
--- /dev/null
+++ b/.github/zizmor.yml
@@ -0,0 +1,7 @@
+# Configuration for zizmor (https://docs.zizmor.sh)
+rules:
+  unpinned-uses:
+    config:
+      # Feel free to switch to hash pinning, then this can be removed.
+      policies:
+        "*": ref-pin
diff --git a/.gitignore b/.gitignore
index 4cc1585..2cf21c7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,23 +1,44 @@
-dist
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
 *.so
+
+# Distribution / packaging
+build/
+dist/
+*.egg-info/
+*.egg
 MANIFEST
-*.py[cod]
-*.egg-info
-.mypy_cache
-.benchmarks
+.eggs/
+
+# Virtual environments
+.venv/
+venv/
+uv.lock
+
+# Unit test / coverage reports
+htmlcov/
+.coverage
+.coverage.*
+coverage.xml
+.pytest_cache/
+.benchmarks/
+.nox/
+
+# Linters and type checkers
+.mypy_cache/
+.ruff_cache/
+
+# Generated files
+src/jacobi/_version.py
+docs/index.rst
+docs/_build/
+
+# Editors and OS files
+.vscode/
+.idea/
 .DS_Store
-_build
-_generate
-build
-__pycache__
-Pipfile.lock
-.eggs
-py37
-var
-venv
-.pytest_cache
 notebook/.ipynb_checkpoints
-src/jacobi/_version.py
-doc/index.rst
-.vscode
-.coverage
diff --git a/.gitmodules b/.gitmodules
deleted file mode 100644
index 47e0f2b..0000000
--- a/.gitmodules
+++ /dev/null
@@ -1,9 +0,0 @@
-[submodule "extern/boost-math"]
-	path = extern/boost-math
-	url = https://github.com/boostorg/math.git
-[submodule "extern/pybind11"]
-	path = extern/pybind11
-	url = https://github.com/pybind/pybind11
-[submodule "extern/boost_math"]
-	path = extern/boost_math
-	url = https://github.com/boostorg/math.git
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index a149dda..8c65c6f 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,52 +1,91 @@
-# To use:
-#
-#     pre-commit run -a
-#
-# Or:
-#
-#     pre-commit install  # (runs every time you commit in git)
-#
-# To update this file:
-#
-#     pre-commit autoupdate
-#
-# See https://github.com/pre-commit/pre-commit
+ci:
+  autoupdate_commit_msg: "chore(deps): update pre-commit hooks"
+  autofix_commit_msg: "style: pre-commit fixes"
+  autoupdate_schedule: monthly
 
 repos:
-# Standard hooks
-- repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.6.0
-  hooks:
-  - id: check-case-conflict
-  - id: check-docstring-first
-  - id: check-executables-have-shebangs
-  - id: check-merge-conflict
-  - id: check-symlinks
-  - id: check-yaml
-  - id: debug-statements
-  - id: end-of-file-fixer
-  - id: mixed-line-ending
-  - id: sort-simple-yaml
-  - id: file-contents-sorter
-  - id: trailing-whitespace
-
-# Ruff linter, replacement for flake8, isort, pydocstyle
-- repo: https://github.com/astral-sh/ruff-pre-commit
-  rev: 'v0.4.5'
-  hooks:
-    - id: ruff
-      args: [--fix, --exit-non-zero-on-fix]
-
-# Python formatting
-- repo: https://github.com/psf/black
-  rev: 24.4.2
-  hooks:
-  - id: black
-
-# Python type checking
-- repo: https://github.com/pre-commit/mirrors-mypy
-  rev: 'v1.10.0'
-  hooks:
-  - id: mypy
-    args: [src]
-    pass_filenames: false
+  - repo: https://github.com/adamchainz/blacken-docs
+    rev: "1.20.0"
+    hooks:
+      - id: blacken-docs
+        additional_dependencies: [black==25.*]
+
+  - repo: https://github.com/pre-commit/pre-commit-hooks
+    rev: v6.0.0
+    hooks:
+      - id: check-added-large-files
+      - id: check-case-conflict
+      - id: check-merge-conflict
+      - id: check-symlinks
+      - id: check-yaml
+      - id: debug-statements
+      - id: end-of-file-fixer
+      - id: mixed-line-ending
+      - id: name-tests-test
+        args: ["--pytest-test-first"]
+      - id: requirements-txt-fixer
+      - id: trailing-whitespace
+
+  - repo: https://github.com/pre-commit/pygrep-hooks
+    rev: v1.10.0
+    hooks:
+      - id: python-use-type-annotations
+      - id: rst-backticks
+      - id: rst-directive-colons
+      - id: rst-inline-touching-normal
+
+  - repo: https://github.com/rbubley/mirrors-prettier
+    rev: "v3.9.6"
+    hooks:
+      - id: prettier
+        types_or: [yaml, markdown, html, css, scss, javascript, json]
+        args: [--prose-wrap=always]
+        exclude: \.ipynb$
+
+  - repo: https://github.com/astral-sh/ruff-pre-commit
+    rev: "v0.16.7"
+    hooks:
+      - id: ruff-check
+        args: ["--fix", "--show-fixes"]
+      - id: ruff-format
+
+  - repo: https://github.com/pre-commit/mirrors-mypy
+    rev: "v2.3.1"
+    hooks:
+      - id: mypy
+        files: ^src
+        args: []
+        # numpy >= 2.5 stubs need Python 3.12 syntax; mypy targets 3.10
+        additional_dependencies: ["numpy~=2.4.0"]
+
+  - repo: https://github.com/codespell-project/codespell
+    rev: "v2.4.3"
+    hooks:
+      - id: codespell
+        additional_dependencies:
+          - tomli; python_version<'3.11'
+
+  - repo: local
+    hooks:
+      - id: disallow-caps
+        name: Disallow improper capitalization
+        language: pygrep
+        entry: PyBind|Numpy|Cmake|CCache|Github|PyTest
+        exclude: .pre-commit-config.yaml
+
+  - repo: https://github.com/abravalheri/validate-pyproject
+    rev: "0.26"
+    hooks:
+      - id: validate-pyproject
+        additional_dependencies: ["validate-pyproject-schema-store[all]"]
+
+  - repo: https://github.com/python-jsonschema/check-jsonschema
+    rev: "0.38.0"
+    hooks:
+      - id: check-dependabot
+      - id: check-github-workflows
+
+  - repo: https://github.com/zizmorcore/zizmor-pre-commit
+    rev: "v1.30.1"
+    hooks:
+      - id: zizmor
diff --git a/README.rst b/README.rst
index 0719e19..1428313 100644
--- a/README.rst
+++ b/README.rst
@@ -6,6 +6,8 @@
 
 .. image:: https://img.shields.io/pypi/v/jacobi
    :target: https://pypi.org/project/jacobi
+.. image:: https://github.com/hdembinski/jacobi/actions/workflows/ci.yml/badge.svg
+   :target: https://github.com/hdembinski/jacobi/actions/workflows/ci.yml
 .. image:: https://img.shields.io/badge/github-docs-success
    :target: https://hdembinski.github.io/jacobi
 .. image:: https://img.shields.io/badge/github-source-blue
@@ -87,6 +89,7 @@ Examples
       r[2] = np.polyval([1, 2, 3], x[0])
       return r  # x and r have different lengths
 
+
   # fn accepts a parameter vector x, which has an associated covariance matrix xcov
   x = [1.0, 2.0]
   xcov = [[1.1, 0.1], [0.1, 2.3]]
diff --git a/tests/bench.py b/bench/test_bench_jacobi.py
similarity index 93%
rename from tests/bench.py
rename to bench/test_bench_jacobi.py
index ebe2be9..156619c 100644
--- a/tests/bench.py
+++ b/bench/test_bench_jacobi.py
@@ -1,6 +1,9 @@
-from jacobi import jacobi
-import pytest
+from __future__ import annotations
+
 import numpy as np
+import pytest
+
+from jacobi import jacobi
 
 fn = ["0", "x ** 2", "x ** -1", "exp(x)", "sin(x)", "ones_like(x)", "mean(x ** 2)"]
 
diff --git a/bench/test_diagonal.py b/bench/test_diagonal.py
index 9edb0b3..2972eba 100644
--- a/bench/test_diagonal.py
+++ b/bench/test_diagonal.py
@@ -1,9 +1,12 @@
-from jacobi import jacobi, propagate
+from __future__ import annotations
+
 import numpy as np
 import pytest
 
+from jacobi import jacobi, propagate
+
 
-@pytest.mark.parametrize("diagonal", (False, True))
+@pytest.mark.parametrize("diagonal", [False, True])
 def test_jacobi(diagonal, benchmark):
     x = np.linspace(0, 1, 100)
 
@@ -16,7 +19,7 @@ def run():
     benchmark(run)
 
 
-@pytest.mark.parametrize("diagonal", (False, True))
+@pytest.mark.parametrize("diagonal", [False, True])
 def test_propagate(diagonal, benchmark):
     x = np.linspace(0, 1, 100)
 
@@ -29,7 +32,7 @@ def run():
     benchmark(run)
 
 
-@pytest.mark.parametrize("diagonal", (False, True))
+@pytest.mark.parametrize("diagonal", [False, True])
 def test_propagate_2dcov(diagonal, benchmark):
     x = np.linspace(0, 1, 100)
 
diff --git a/doc/build.py b/doc/build.py
deleted file mode 100755
index fd5e8d5..0000000
--- a/doc/build.py
+++ /dev/null
@@ -1,8 +0,0 @@
-#!/usr/bin/env python3
-import subprocess as subp
-from pathlib import Path
-import sys
-
-cmd = "sphinx-build -W -b html -d _build/doctrees . _build/html"
-
-sys.exit(subp.call(cmd.split(), cwd=Path(__file__).parent))
diff --git a/doc/conf.py b/doc/conf.py
deleted file mode 100644
index 948f131..0000000
--- a/doc/conf.py
+++ /dev/null
@@ -1,69 +0,0 @@
-# -*- coding: utf-8 -*-
-#
-# python_example documentation build configuration file, created by
-# sphinx-quickstart on Fri Feb 26 00:29:33 2016.
-#
-# This file is execfile()d with the current directory set to its
-# containing dir.
-#
-# Note that not all possible configuration values are present in this
-# autogenerated file.
-#
-# All configuration values have a default; values that are commented out
-# serve to show the default.
-
-from pathlib import Path
-
-with open("index.rst.in") as f:
-    stub = f.read()
-
-with open(Path("..") / "README.rst") as f:
-    readme = f.read().replace(
-        "https://hdembinski.github.io/jacobi/_images/", "_static/"
-    )
-
-with open("index.rst", "w") as f:
-    f.write(stub + "\n" + readme)
-
-extensions = [
-    "sphinx.ext.autodoc",
-    "sphinx.ext.intersphinx",
-    "sphinx.ext.autosummary",
-    "sphinx.ext.napoleon",
-    "sphinx.ext.doctest",
-]
-
-autosummary_generate = True
-
-# General information about the project.
-project = "Jacobi"
-copyright = "2021, Hans Dembinski"
-author = "Hans Dembinski"
-
-# The version info for the project you're documenting, acts as replacement for
-# |version| and |release|, also used in various other places throughout the
-# built documents.
-#
-with open(Path() / ".." / "src" / "jacobi" / "_version.py") as f:
-    d = {}
-    exec(f.read(), d)
-
-# The full version, including alpha/beta/rc tags.
-release = d["version"]
-
-# The short X.Y version.
-version = ".".join(release.split(".")[:2])
-
-language = "en"
-
-exclude_patterns = ["_build"]
-
-pygments_style = "sphinx"
-
-# -- Options for HTML output ----------------------------------------------
-
-# The theme to use for HTML and HTML Help pages.  See the documentation for
-# a list of builtin themes.
-html_theme = "sphinx_rtd_theme"
-
-html_static_path = ["_static"]
diff --git a/doc/_static/example.svg b/docs/_static/example.svg
similarity index 100%
rename from doc/_static/example.svg
rename to docs/_static/example.svg
diff --git a/doc/_static/logo.svg b/docs/_static/logo.svg
similarity index 100%
rename from doc/_static/logo.svg
rename to docs/_static/logo.svg
diff --git a/doc/_static/precision.svg b/docs/_static/precision.svg
similarity index 100%
rename from doc/_static/precision.svg
rename to docs/_static/precision.svg
diff --git a/doc/_static/speed.svg b/docs/_static/speed.svg
similarity index 100%
rename from doc/_static/speed.svg
rename to docs/_static/speed.svg
diff --git a/docs/conf.py b/docs/conf.py
new file mode 100644
index 0000000..a0d462c
--- /dev/null
+++ b/docs/conf.py
@@ -0,0 +1,64 @@
+"""Sphinx configuration."""
+
+from __future__ import annotations
+
+import importlib.metadata
+from pathlib import Path
+
+DIR = Path(__file__).parent.resolve()
+
+# index.rst is generated from index.rst.in and the README, so that the
+# front page of the documentation and the README stay in sync.
+readme = (DIR.parent / "README.rst").read_text(encoding="utf-8")
+readme = readme.replace("https://hdembinski.github.io/jacobi/_images/", "_static/")
+stub = (DIR / "index.rst.in").read_text(encoding="utf-8")
+index = DIR / "index.rst"
+content = stub + "\n" + readme
+# Only write when needed, otherwise sphinx-autobuild rebuilds in a loop.
+if not index.is_file() or index.read_text(encoding="utf-8") != content:
+    index.write_text(content, encoding="utf-8")
+
+project = "jacobi"
+copyright = "2020, Hans Dembinski"
+author = "Hans Dembinski"
+version = release = importlib.metadata.version("jacobi")
+
+extensions = [
+    "sphinx.ext.autodoc",
+    "sphinx.ext.doctest",
+    "sphinx.ext.intersphinx",
+    "sphinx.ext.mathjax",
+    "sphinx.ext.napoleon",
+    "sphinx_copybutton",
+]
+
+exclude_patterns = [
+    "_build",
+    "**.ipynb_checkpoints",
+    "Thumbs.db",
+    ".DS_Store",
+    ".env",
+    ".venv",
+]
+
+# Types are documented in the numpy-style docstrings.
+autodoc_typehints = "none"
+
+napoleon_google_docstring = False
+napoleon_numpy_docstring = True
+# Render parameter and return types verbatim instead of cross-referencing them.
+napoleon_use_param = False
+napoleon_use_rtype = False
+
+intersphinx_mapping = {
+    "python": ("https://docs.python.org/3", None),
+    "numpy": ("https://numpy.org/doc/stable/", None),
+}
+
+html_theme = "furo"
+html_static_path = ["_static"]
+html_theme_options = {
+    "source_repository": "https://github.com/hdembinski/jacobi",
+    "source_branch": "main",
+    "source_directory": "docs/",
+}
diff --git a/doc/index.rst.in b/docs/index.rst.in
similarity index 100%
rename from doc/index.rst.in
rename to docs/index.rst.in
diff --git a/doc/plot/example.py b/docs/plot/example.py
similarity index 90%
rename from doc/plot/example.py
rename to docs/plot/example.py
index 5d92339..b91aaff 100644
--- a/doc/plot/example.py
+++ b/docs/plot/example.py
@@ -1,5 +1,8 @@
-from matplotlib import pyplot as plt
+from __future__ import annotations
+
 import numpy as np
+from matplotlib import pyplot as plt
+
 from jacobi import jacobi
 
 
@@ -28,4 +31,4 @@ def f(x):
 )
 plt.legend()
 
-plt.savefig("doc/_static/example.svg")
+plt.savefig("docs/_static/example.svg")
diff --git a/doc/plot/precision.py b/docs/plot/precision.py
similarity index 92%
rename from doc/plot/precision.py
rename to docs/plot/precision.py
index c1df2d5..2d96151 100644
--- a/doc/plot/precision.py
+++ b/docs/plot/precision.py
@@ -1,8 +1,11 @@
-from matplotlib import pyplot as plt
+from __future__ import annotations
+
 import numpy as np
-from jacobi import jacobi
+from matplotlib import pyplot as plt
 from numdifftools import Derivative
 
+from jacobi import jacobi
+
 
 # function of one variable with auxiliary argument; returns a vector
 def f(x):
@@ -29,4 +32,4 @@ def fp(x):
 plt.ylabel("$|f'_\\mathrm{num}(x) - f'(x)|$")
 plt.axhline(np.finfo(float).resolution, color="k", ls="--")
 
-plt.savefig("doc/_static/precision.svg")
+plt.savefig("docs/_static/precision.svg")
diff --git a/doc/plot/speed.py b/docs/plot/speed.py
similarity index 78%
rename from doc/plot/speed.py
rename to docs/plot/speed.py
index 820a06d..d0972b6 100644
--- a/doc/plot/speed.py
+++ b/docs/plot/speed.py
@@ -1,8 +1,12 @@
-from jacobi import jacobi
-import numpy as np
-from numdifftools import Derivative
+from __future__ import annotations
+
 from timeit import timeit
+
+import numpy as np
 from matplotlib import pyplot as plt
+from numdifftools import Derivative
+
+from jacobi import jacobi
 
 fn = [
     "x ** 2",
@@ -25,10 +29,13 @@
         print(i, ni)
         x = np.linspace(0.1, 10, ni)
         number = 500 // ni + 10
-        r = timeit(lambda: jacobi(f, x, diagonal=True), number=number) / number
+        r = timeit(lambda f=f, x=x: jacobi(f, x, diagonal=True), number=number) / number
         t["jacobi"].append(r)
         number = 500 // ni + 1
-        r = timeit(lambda: Derivative(lambda p: f(p + x))(0), number=number) / number
+        r = (
+            timeit(lambda f=f, x=x: Derivative(lambda p: f(p + x))(0), number=number)
+            / number
+        )
         t["numdifftools"].append(r)
     plt.sca(ax.flat[i])
     for k, v in t.items():
@@ -36,7 +43,7 @@
         plt.plot(n, v, ls=ls, label=k)
     r = np.divide(t["numdifftools"], t["jacobi"])
     plt.plot(n, r, ls=":", color="k", lw=3, label="time ratio")
-    for ni, ri in zip(n, r):
+    for ni, ri in zip(n, r, strict=True):
         plt.text(
             ni, ri * 0.5, f"{ri:.0f}" if i != 3 else f"{ri:.1f}", va="top", ha="center"
         )
@@ -46,4 +53,4 @@
 fig.supxlabel("n")
 fig.supylabel("t/sec")
 
-plt.savefig("doc/_static/speed.svg")
+plt.savefig("docs/_static/speed.svg")
diff --git a/doc/reference.rst b/docs/reference.rst
similarity index 100%
rename from doc/reference.rst
rename to docs/reference.rst
diff --git a/noxfile.py b/noxfile.py
new file mode 100755
index 0000000..f91b26e
--- /dev/null
+++ b/noxfile.py
@@ -0,0 +1,138 @@
+#!/usr/bin/env -S uv run --script
+
+# /// script
+# dependencies = ["nox>=2025.2.9"]
+# ///
+
+"""Nox runner."""
+
+from __future__ import annotations
+
+import argparse
+import shutil
+from pathlib import Path
+
+import nox
+
+DIR = Path(__file__).parent.resolve()
+PROJECT = nox.project.load_toml()
+MINIMUM_PYTHON = PROJECT["project"]["requires-python"].removeprefix(">=")
+
+nox.needs_version = ">=2025.2.9"
+nox.options.default_venv_backend = "uv|virtualenv"
+
+
+@nox.session
+def lint(session: nox.Session) -> None:
+    """Run the linter."""
+    session.install("prek")
+    session.run(
+        "prek", "run", "--all-files", "--show-diff-on-failure", *session.posargs
+    )
+
+
+@nox.session
+def pylint(session: nox.Session) -> None:
+    """Run Pylint."""
+    # This needs to be installed into the package environment, and is slower
+    # than a pre-commit check
+    session.install("-e.", "pylint>=3.2")
+    session.run("pylint", "jacobi", *session.posargs)
+
+
+@nox.session
+def tests(session: nox.Session) -> None:
+    """Run the unit and regular tests."""
+    test_deps = nox.project.dependency_groups(PROJECT, "test")
+    session.install("-e.", *test_deps)
+    session.run("pytest", *session.posargs)
+
+
+@nox.session(python=MINIMUM_PYTHON, venv_backend="uv", default=False)
+def minimums(session: nox.Session) -> None:
+    """Run the tests with the minimum versions of the direct dependencies."""
+    test_deps = nox.project.dependency_groups(PROJECT, "test")
+    session.install("-e.", *test_deps, "--resolution=lowest-direct")
+    session.run("uv", "pip", "list")
+    session.run("pytest", *session.posargs)
+
+
+@nox.session(default=False)
+def cov(session: nox.Session) -> None:
+    """Run the tests with coverage and write an HTML report to "htmlcov"."""
+    test_deps = nox.project.dependency_groups(PROJECT, "test")
+    session.install("-e.", *test_deps)
+    session.run(
+        "pytest",
+        "--cov",
+        "--cov-report=term-missing",
+        "--cov-report=html",
+        *session.posargs,
+    )
+
+
+@nox.session(default=False)
+def bench(session: nox.Session) -> None:
+    """Run the benchmarks."""
+    bench_deps = nox.project.dependency_groups(PROJECT, "bench")
+    session.install("-e.", *bench_deps)
+    session.run("pytest", "bench", *session.posargs)
+
+
+@nox.session(reuse_venv=True, default=False)
+def docs(session: nox.Session) -> None:
+    """
+    Build the docs.
+
+    Serves the docs when run interactively; pass --non-interactive to nox to
+    avoid serving. The first positional argument is the target directory.
+    """
+    doc_deps = nox.project.dependency_groups(PROJECT, "docs")
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "-b", dest="builder", default="html", help="Build target (default: html)"
+    )
+    parser.add_argument("output", nargs="?", help="Output directory")
+    args, posargs = parser.parse_known_args(session.posargs)
+    serve = args.builder == "html" and session.interactive
+
+    session.install("-e.", *doc_deps, *(["sphinx-autobuild"] if serve else []))
+
+    shared_args = (
+        "-n",  # nitpicky mode
+        "-T",  # full tracebacks
+        "-W",  # warnings are errors
+        f"-b={args.builder}",
+        "docs",
+        args.output or f"docs/_build/{args.builder}",
+        *posargs,
+    )
+
+    if serve:
+        session.run("sphinx-autobuild", "--open-browser", *shared_args)
+    else:
+        session.run("sphinx-build", "--keep-going", *shared_args)
+
+
+@nox.session(reuse_venv=True, default=False)
+def plots(session: nox.Session) -> None:
+    """Regenerate the figures in docs/_static."""
+    plot_deps = nox.project.dependency_groups(PROJECT, "plot")
+    session.install("-e.", *plot_deps)
+    for script in sorted((DIR / "docs" / "plot").glob("*.py")):
+        session.run("python", str(script))
+
+
+@nox.session(default=False)
+def build(session: nox.Session) -> None:
+    """Build an SDist and wheel."""
+    build_path = DIR.joinpath("build")
+    if build_path.exists():
+        shutil.rmtree(build_path)
+
+    session.install("build")
+    session.run("python", "-m", "build")
+
+
+if __name__ == "__main__":
+    nox.main()
diff --git a/pyproject.toml b/pyproject.toml
index fd164af..6e99646 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,66 +1,191 @@
 [build-system]
-requires = ["setuptools>=42", "setuptools_scm[toml]>=3.4"]
-build-backend = "setuptools.build_meta"
+requires = ["hatchling>=1.27", "hatch-vcs"]
+build-backend = "hatchling.build"
 
 [project]
 name = "jacobi"
-requires-python = ">=3.8"
 description = "Fast numerical derivatives for analytic functions with arbitrary round-off error and error propagation"
-authors = [{ name = "Hans Dembinski" }, { email = "hans.dembinski@gmail.com" }]
-dynamic = ["version"]
-dependencies = ["numpy"]
 readme = "README.rst"
+authors = [{ name = "Hans Dembinski", email = "hans.dembinski@gmail.com" }]
+maintainers = [
+  { name = "Hans Dembinski", email = "hans.dembinski@gmail.com" },
+  { name = "The Scikit-HEP admins", email = "scikit-hep-admins@googlegroups.com" },
+]
+license = "MIT"
+license-files = ["LICENSE"]
+requires-python = ">=3.10"
+keywords = [
+  "numerical differentiation",
+  "derivative",
+  "jacobian",
+  "error propagation",
+]
 classifiers = [
-    "Development Status :: 4 - Beta",
-    "License :: OSI Approved :: MIT License",
-    "Programming Language :: Python :: 3",
-    "Operating System :: OS Independent",
-    "Operating System :: POSIX :: Linux",
-    "Operating System :: MacOS",
-    "Operating System :: Microsoft :: Windows",
-    "Topic :: Scientific/Engineering",
+  "Development Status :: 4 - Beta",
+  "Intended Audience :: Developers",
+  "Intended Audience :: Science/Research",
+  "Operating System :: OS Independent",
+  "Programming Language :: Python",
+  "Programming Language :: Python :: 3",
+  "Programming Language :: Python :: 3 :: Only",
+  "Programming Language :: Python :: 3.10",
+  "Programming Language :: Python :: 3.11",
+  "Programming Language :: Python :: 3.12",
+  "Programming Language :: Python :: 3.13",
+  "Programming Language :: Python :: 3.14",
+  "Programming Language :: Python :: Implementation :: CPython",
+  "Programming Language :: Python :: Implementation :: PyPy",
+  "Topic :: Scientific/Engineering",
+  "Topic :: Scientific/Engineering :: Mathematics",
+  "Typing :: Typed",
 ]
-license = { file = "LICENSE" }
+dynamic = ["version"]
+dependencies = ["numpy>=1.24"]
 
 [project.urls]
-repository = "https://github.com/hdembinski/jacobi"
-documentation = "https://hdembinski.github.io/jacobi/"
+Homepage = "https://github.com/hdembinski/jacobi"
+Documentation = "https://hdembinski.github.io/jacobi/"
+Repository = "https://github.com/hdembinski/jacobi"
+"Bug Tracker" = "https://github.com/hdembinski/jacobi/issues"
+Changelog = "https://github.com/hdembinski/jacobi/releases"
 
-[project.optional-dependencies]
-test = ["pytest", "pytest-benchmark"]
-doc = ["sphinx", "sphinx-rtd-theme", "ipykernel"]
-plot = ["numdifftools", "matplotlib"]
+[dependency-groups]
+test = [
+  "pytest>=9",
+  "pytest-cov>=7",
+]
+bench = [
+  { include-group = "test" },
+  "pytest-benchmark>=4",
+]
+plot = [
+  "matplotlib>=3.5",
+  "numdifftools>=0.9.41",
+]
+docs = [
+  "sphinx>=7",
+  "furo>=2023.8.17",
+  "sphinx-copybutton",
+]
+dev = [
+  { include-group = "test" },
+]
 
-[tool.setuptools.packages.find]
-where = ["src"]
+[tool.hatch]
+version.source = "vcs"
+build.hooks.vcs.version-file = "src/jacobi/_version.py"
 
-[tool.setuptools_scm]
-write_to = "src/jacobi/_version.py"
+[tool.hatch.envs.default]
+dependency-groups = ["dev"]
 
-[tool.mypy]
-ignore_missing_imports = true
-allow_redefinition = true
-no_implicit_optional = false
+[tool.hatch.envs.hatch-test]
+dependency-groups = ["test"]
+env-vars.PYTHONWARNDEFAULTENCODING = "1"
 
-[tool.pytest.ini_options]
-minversion = "6.0"
-addopts = "--doctest-modules --strict-config --strict-markers -q -ra --ff"
+[tool.hatch.build.targets.sdist]
+# The exploratory notebooks are large and not part of the package.
+exclude = ["/notebook"]
+
+[tool.pytest]
+minversion = "9.0"
+addopts = ["-ra", "--showlocals", "--doctest-modules"]
+strict = true
+filterwarnings = ["error"]
+log_level = "INFO"
 testpaths = ["src/jacobi", "tests"]
-filterwarnings = [
-    "error::numpy.VisibleDeprecationWarning",
-    "error::DeprecationWarning",
+
+[tool.coverage]
+run.source = ["jacobi"]
+report.exclude_also = [
+  'if typing.TYPE_CHECKING:',
+  'if TYPE_CHECKING:',
 ]
 
+[tool.mypy]
+files = ["src"]
+python_version = "3.10"
+warn_unused_configs = true
+strict = true
+enable_error_code = ["ignore-without-code", "redundant-expr", "truthy-bool"]
+warn_unreachable = true
+
 [tool.ruff]
-select = ["E", "F", "D"]
-extend-ignore = ["D203", "D212"]
+show-fixes = true
+# The exploratory notebooks are not maintained as code.
+extend-exclude = ["notebook"]
+
+[tool.ruff.lint]
+extend-select = [
+  "ARG",  # flake8-unused-arguments
+  "B",    # flake8-bugbear
+  "BLE",  # flake8-blind-except
+  "C4",   # flake8-comprehensions
+  "D",    # pydocstyle
+  "DTZ",  # flake8-datetimez
+  "E",    # pycodestyle errors
+  "EM",   # flake8-errmsg
+  "EXE",  # flake8-executable
+  "FA",   # flake8-future-annotations
+  "FLY",  # flynt
+  "FURB", # refurb
+  "G",    # flake8-logging-format
+  "I",    # isort
+  "ICN",  # flake8-import-conventions
+  "ISC",  # flake8-implicit-str-concat
+  "LOG",  # flake8-logging
+  "NPY",  # NumPy specific rules
+  "PERF", # perflint
+  "PGH",  # pygrep-hooks
+  "PIE",  # flake8-pie
+  "PL",   # pylint
+  "PT",   # flake8-pytest-style
+  "PTH",  # flake8-use-pathlib
+  "PYI",  # flake8-pyi
+  "Q",    # flake8-quotes
+  "RET",  # flake8-return
+  "RSE",  # flake8-raise
+  "RUF",  # Ruff-specific
+  "SIM",  # flake8-simplify
+  "SLOT", # flake8-slots
+  "T10",  # flake8-debugger
+  "T20",  # flake8-print
+  "TC",   # flake8-type-checking
+  "TRY",  # tryceratops
+  "UP",   # pyupgrade
+  "W",    # pycodestyle warnings
+  "YTT",  # flake8-2020
+]
+ignore = [
+  "ISC001",  # Conflicts with the formatter
+  "PLR09",   # Too many <...>
+  "PLR2004", # Magic value used in comparison
+  "PT011",   # pytest.raises(ValueError) is too broad
+]
+isort.required-imports = ["from __future__ import annotations"]
+flake8-unused-arguments.ignore-variadic-names = true
 
-[tool.ruff.pydocstyle]
+[tool.ruff.lint.pydocstyle]
 convention = "numpy"
 
-[tool.ruff.per-file-ignores]
-"test_*.py" = ["B", "D"]
-"tests/bench.py" = ["D"]
-".ci/*.py" = ["D"]
-"bench/*.py" = ["D"]
-"doc/*.py" = ["D"]
+[tool.ruff.lint.per-file-ignores]
+"tests/**" = ["ARG", "D", "PLR", "T20"]
+"bench/**" = ["ARG", "D", "PLR", "T20"]
+"docs/**" = ["D", "T20"]
+"noxfile.py" = ["T20"]
+
+[tool.pylint]
+py-version = "3.10"
+ignore-paths = [".*/_version.py"]
+reports.output-format = "colorized"
+similarities.ignore-imports = "yes"
+messages_control.disable = [
+  "design",
+  "fixme",
+  "line-too-long",
+  "missing-module-docstring",
+  "missing-function-docstring",
+  "wrong-import-position",
+]
+
+[tool.codespell]
+skip = "*.svg,*.ipynb,*.lock"
diff --git a/src/jacobi/__init__.py b/src/jacobi/__init__.py
index 5a31d72..7a1b16b 100644
--- a/src/jacobi/__init__.py
+++ b/src/jacobi/__init__.py
@@ -1,7 +1,9 @@
 """Fast numerical derivatives for analytic functions with arbitrary round-off error."""
 
-from ._propagate import propagate
+from __future__ import annotations
+
 from ._jacobi import jacobi
-from ._version import version as __version__  # noqa
+from ._propagate import propagate
+from ._version import version as __version__
 
-__all__ = ["jacobi", "propagate"]
+__all__ = ["__version__", "jacobi", "propagate"]
diff --git a/src/jacobi/_jacobi.py b/src/jacobi/_jacobi.py
index 09972f0..e8b631d 100644
--- a/src/jacobi/_jacobi.py
+++ b/src/jacobi/_jacobi.py
@@ -1,21 +1,30 @@
+"""Numerical computation of the Jacobi matrix with error estimates."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Any
+
 import numpy as np
-import typing as _tp
-from ._typing import Indexable as _Indexable
+
+if TYPE_CHECKING:
+    from collections.abc import Callable
+
+    from numpy.typing import ArrayLike, NDArray
 
 
 def jacobi(
-    fn: _tp.Callable,
-    x: _tp.Union[float, _Indexable[float]],
-    *args,
+    fn: Callable[..., ArrayLike],
+    x: ArrayLike,
+    *args: Any,
     diagonal: bool = False,
-    method: int = None,
-    mask: np.ndarray = None,
+    method: int | None = None,
+    mask: ArrayLike | None = None,
     rtol: float = 0,
     maxiter: int = 10,
     maxgrad: int = 3,
-    step: _tp.Tuple[float, float] = None,
-    diagnostic: dict = None,
-) -> _tp.Tuple[np.ndarray, np.ndarray]:
+    step: tuple[float, float] | None = None,
+    diagnostic: dict[str, Any] | None = None,
+) -> tuple[NDArray[Any], NDArray[Any]]:
     """
     Return first derivative and its error estimate.
 
@@ -33,11 +42,11 @@ def jacobi(
     *args : tuple
         Additional arguments passed to the function.
     diagonal : boolean, optional
-        If it is known that the Jacobian computed from the function contains has
-        off-diagonal entries that are all zero, the calculation can be speed up
-        significantly. Set this to true to only compute the diagonal entries of the
-        Jacobi matrix, which are returned as a 1D array. This is faster and uses much
-        less memory if the vector x is very large. Default is False.
+        If it is known that the Jacobian computed from the function has off-diagonal
+        entries that are all zero, the calculation can be sped up significantly. Set
+        this to true to only compute the diagonal entries of the Jacobi matrix, which
+        are returned as a 1D array. This is faster and uses much less memory if the
+        vector x is very large. Default is False.
     method : {-1, 0, 1} or None, optional
         Whether to compute central (0), forward (1) or backward derivatives (-1).
         The default (None) uses auto-detection.
@@ -51,7 +60,7 @@ def jacobi(
     maxiter : int, optional
         Maximum number of iterations of the algorithm.
     maxgrad : int, optional
-        Maximum grad of the extrapolation polynomial.
+        Maximum degree of the extrapolation polynomial.
     step : tuple of float or None, optional
         Factors that reduce the step size in each iteration relative to the previous
         step.
@@ -65,6 +74,40 @@ def jacobi(
     array, array
         Derivative and its error estimate.
     """
+    if maxiter <= 0:
+        msg = "maxiter must be > 0"
+        raise ValueError(msg)
+    if maxgrad < 0:
+        msg = "maxgrad must be >= 0"
+        raise ValueError(msg)
+    if step is not None:
+        if not 0 < step[0] < 0.5:
+            msg = "step[0] must be between 0 and 0.5"
+            raise ValueError(msg)
+        if not 0 < step[1] < 1:
+            msg = "step[1] must be between 0 and 1"
+            raise ValueError(msg)
+    if method is not None and method not in (-1, 0, 1):
+        msg = "method must be -1, 0, 1"
+        raise ValueError(msg)
+
+    xa = np.asarray(x, dtype=float)
+    ma: NDArray[Any] | None = None
+    if mask is not None:
+        ma = np.asarray(mask)
+        if ma.dtype != bool:
+            msg = "mask must be a boolean array"
+            raise ValueError(msg)
+        if ma.shape != xa.shape:
+            msg = "mask shape must match x shape"
+            raise ValueError(msg)
+
+    if xa.size == 0 or (ma is not None and not ma.any()):
+        # nothing to compute, but the output shape is needed for the result
+        _, fval = _wrap_function_if_needed(fn, fn(xa, *args))
+        shape = np.shape(fval) if diagonal else np.shape(fval) + xa.shape
+        return np.zeros(shape), np.zeros(shape)
+
     if diagonal:
         # TODO maybe solve this without introducing a wrapper function
         j, je = jacobi(
@@ -79,77 +122,57 @@ def jacobi(
             step=step,
             diagnostic=diagnostic,
         )
-        if mask is not None:
-            j[~mask] = 0.0
-            je[~mask] = 0.0
+        if ma is not None:
+            j[~ma] = 0.0
+            je[~ma] = 0.0
         return j, je
 
-    if maxiter <= 0:
-        raise ValueError("maxiter must be > 0")
-    if maxgrad < 0:
-        raise ValueError("maxgrad must be >= 0")
-    if step is not None:
-        if not (0 < step[0] < 0.5):
-            raise ValueError("step[0] must be between 0 and 0.5")
-        if not (0 < step[1] < 1):
-            raise ValueError("step[1] must be between 0 and 1")
-    if method is not None and method not in (-1, 0, 1):
-        raise ValueError("method must be -1, 0, 1")
-
-    x = np.asarray(x, dtype=float)
-    if mask is not None:
-        mask = np.asarray(mask)
-        if mask.dtype != bool:
-            raise ValueError("mask must be a boolean array")
-        if mask.shape != x.shape:
-            raise ValueError("mask shape must match x shape")
-
     if diagnostic is not None:
-        diagnostic["method"] = np.zeros(x.size, dtype=np.int8)
-        diagnostic["iteration"] = np.zeros(x.size, dtype=np.uint8)
-        diagnostic["residual"] = [[] for _ in range(x.size)]
-
-    f0 = None
-    jac = None
-    err = None
-    it = np.nditer(x, flags=["c_index", "multi_index"])
+        diagnostic["method"] = np.zeros(xa.size, dtype=np.int8)
+        diagnostic["iteration"] = np.zeros(xa.size, dtype=np.uint8)
+        diagnostic["residual"] = [[] for _ in range(xa.size)]
+
+    f0: Any = None
+    jac: NDArray[Any] | None = None
+    err: NDArray[Any] | None = None
+    it = np.nditer(xa, flags=["c_index", "multi_index"])
     while not it.finished:
         k = it.index
-        kx = it.multi_index if it.has_multi_index else ...
-        if mask is not None and not mask[kx]:
+        kx = it.multi_index
+        if ma is not None and not ma[kx]:
             it.iternext()
             continue
         xk = it[0]
         # if step is None, use optimal step sizes for central derivatives
         h = _steps(xk, step or (0.25, 0.5), maxiter)
         # if method is None, auto-detect for each x[k]
-        fn, md, f0, r = _first(method, f0, fn, x, kx, h[0], args)
+        fn, md, f0, r = _first(method, f0, fn, xa, kx, h[0], args)
         # f0 is not guaranteed to be set here and can be still None
 
         if md != 0 and step is None:
             # need different step sizes for forward derivatives to avoid overlap
             h = _steps(xk, (0.25, 0.125), maxiter)
 
-        r = np.asarray(r, dtype=float)
-        re = np.full_like(r, np.inf)
-        todo = np.ones_like(r, dtype=bool)
-        fd = [np.reshape(r.copy(), -1)]
+        res = np.asarray(r, dtype=float)
+        res_err = np.full_like(res, np.inf)
+        todo = np.ones_like(res, dtype=bool)
+        fd = [np.reshape(res.copy(), -1)]
 
-        if jac is None:  # first iteration
-            jac = np.zeros(r.shape + x.shape, dtype=r.dtype)
-            err = np.zeros(r.shape + x.shape, dtype=r.dtype)
+        if jac is None or err is None:  # first iteration
+            jac = np.zeros(res.shape + xa.shape, dtype=res.dtype)
+            err = np.zeros(res.shape + xa.shape, dtype=res.dtype)
             if diagnostic is not None:
-                diagnostic["call"] = np.zeros((r.size, x.size), dtype=np.uint8)
+                diagnostic["call"] = np.zeros((res.size, xa.size), dtype=np.uint8)
 
         if diagnostic is not None:
             diagnostic["method"][k] = md
             diagnostic["call"][:, k] = 2 if md == 0 else 3
 
         for i in range(1, len(h)):
-            fdi = _derive(md, f0, fn, x, kx, h[i], args)
+            fdi = np.asarray(_derive(md, f0, fn, xa, kx, h[i], args))
             fd.append(np.reshape(fdi, -1) if i == 1 else fdi[todo])
             if diagnostic is not None:
-                diagnostic["call"][todo, k] += 2
+                diagnostic["call"][todo.reshape(-1), k] += 2
                 diagnostic["iteration"][k] += 1
 
             # polynomial fit with one extra degree of freedom;
@@ -167,11 +190,11 @@ def jacobi(
             rei = c[-1, -1] ** 0.5
 
             # update estimates that have smaller estimated error
-            sub_todo = rei < re[todo]
+            sub_todo = rei < res_err[todo]
             todo1 = todo.copy()
             todo[todo1] = sub_todo
-            r[todo] = ri[sub_todo]
-            re[todo] = rei[sub_todo]
+            res[todo] = ri[sub_todo]
+            res_err[todo] = rei[sub_todo]
 
             # do not improve estimates further which meet the tolerance
             if rtol > 0:
@@ -179,7 +202,7 @@ def jacobi(
                 todo[todo1] = sub_todo
 
             if diagnostic is not None:
-                re2 = re.copy()
+                re2 = res_err.copy()
                 re2[todo1] = rei
                 diagnostic["residual"][k].append(re2)
 
@@ -187,32 +210,36 @@ def jacobi(
                 break
 
             # shrink previous vectors of estimates
-            fd = [fdi[sub_todo] for fdi in fd]
-
-        if jac.ndim == 0:
-            jac[...] = r
-            err[...] = re
-        elif jac.ndim == 1:
-            jac[kx] = r
-            err[kx] = re
-        else:
-            jac[(...,) + kx] = r
-            err[(...,) + kx] = re
+            fd = [v[sub_todo] for v in fd]
+
+        idx: tuple[Any, ...] = (..., *kx)
+        jac[idx] = res
+        err[idx] = res_err
 
         it.iternext()
 
+    assert jac is not None
+    assert err is not None
     return jac, err
 
 
-def _steps(p, step, maxiter):
+def _steps(p: Any, step: tuple[float, float], maxiter: int) -> NDArray[Any]:
     h0, factor = step
     h = p * h0
-    if not h != 0:  # also works if p is NaN
+    if h == 0:  # if p is NaN, h stays NaN
         h = h0
-    return h * factor ** np.arange(maxiter)
+    return np.asarray(h * factor ** np.arange(maxiter))
 
 
-def _derive(mode, f0, f, x, i, h, args):
+def _derive(
+    mode: int,
+    f0: Any,
+    f: Callable[..., Any],
+    x: NDArray[Any],
+    i: tuple[int, ...],
+    h: Any,
+    args: tuple[Any, ...],
+) -> Any:
     x1 = x.copy()
     x2 = x.copy()
     if mode == 0:
@@ -227,13 +254,21 @@ def _derive(mode, f0, f, x, i, h, args):
     return (-3 * f0 + 4 * f1 - f2) * (0.5 / h)
 
 
-def _first(method, f0, fn, x, i, h, args):
+def _first(
+    method: int | None,
+    f0: Any,
+    fn: Callable[..., Any],
+    x: NDArray[Any],
+    i: tuple[int, ...],
+    h: Any,
+    args: tuple[Any, ...],
+) -> tuple[Callable[..., Any], int, Any, Any]:
     # This is the first derivative that we calculate.
     # This function is special because we collect a lot of diagnostic
     # information about the function for the remainder of the iterations.
     norm = 0.5 / h
-    f1 = None
-    f2 = None
+    f1: Any = None
+    f2: Any = None
     if method is None or method == 0:
         x1 = x.copy()
         x2 = x.copy()
@@ -269,14 +304,21 @@ def _first(method, f0, fn, x, i, h, args):
     return fn, method, f0, (-3 * f0 + 4 * f1 - f2) * norm
 
 
-def _wrap_function_if_needed(fn, fval):
+def _wrap_function_if_needed(
+    fn: Callable[..., Any], fval: Any
+) -> tuple[Callable[..., Any], Any]:
     if not isinstance(fval, float):
         try:
-            fval = np.asarray(fval, dtype=float)
+            fval_a = np.asarray(fval, dtype=float)
         except ValueError as e:
-            raise ValueError(
+            msg = (
                 "function return value cannot be converted into "
                 "1D numpy array of floats"
-            ) from e
-        return lambda *args: np.asarray(fn(*args)), fval
+            )
+            raise ValueError(msg) from e
+        if isinstance(fval, np.ndarray):
+            # fn already returns arrays and needs no wrapper; wrapping it again
+            # for every element of x exceeds the recursion limit for large x
+            return fn, fval_a
+        return lambda *args: np.asarray(fn(*args)), fval_a
     return fn, fval
diff --git a/src/jacobi/_propagate.py b/src/jacobi/_propagate.py
index 1092a55..7014513 100644
--- a/src/jacobi/_propagate.py
+++ b/src/jacobi/_propagate.py
@@ -1,21 +1,30 @@
-from typing import Callable, Union, Tuple, List
-from ._typing import Indexable
-from ._jacobi import jacobi
+"""Numerical error propagation based on the Jacobi matrix."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Any
+
 import numpy as np
 
+from ._jacobi import jacobi
+
+if TYPE_CHECKING:
+    from collections.abc import Callable
+
+    from numpy.typing import ArrayLike, NDArray
 
 __all__ = ["propagate"]
 
 
 def propagate(
-    fn: Callable,
-    x: Union[float, Indexable[float]],
-    cov: Union[float, Indexable[float], Indexable[Indexable[float]]],
-    *args,
-    **kwargs,
-) -> Tuple[np.ndarray, np.ndarray]:
+    fn: Callable[..., ArrayLike],
+    x: ArrayLike,
+    cov: ArrayLike,
+    *args: ArrayLike,
+    **kwargs: Any,
+) -> tuple[NDArray[Any], NDArray[Any]]:
     """
-    Numerically propagates the covariance of function inputs to function outputs.
+    Numerically propagate the covariance of function inputs to function outputs.
 
     The function computes C' = J C J^T, where C is the covariance matrix of the input,
     C' the matrix of the output, and J is the Jacobi matrix of first derivatives of the
@@ -23,23 +32,23 @@ def propagate(
 
     Parameters
     ----------
-    fn: callable
+    fn : callable
         Function with the signature `fn(x, [y, ...])`, where `x` is a number or a
         sequence of numbers, likewise if other arguments are present they must have the
         same format. The function must return a number or a sequence of numbers (ideally
         as a numpy array). The length of `x` can differ from the output sequence. The
         function should accept more than one argument only if there are no
         correlations between these arguments. See example below for use cases.
-    x: float or array-like with shape (N,)
+    x : float or array-like with shape (N,)
         Input vector. An array-like is converted before passing it to the callable.
-    cov: float or array-like with shape (N,) or shape(N, N)
+    cov : float or array-like with shape (N,) or shape (N, N)
         Covariance matrix of input vector. If the array is one-dimensional, it is
         interpreted as the diagonal of a covariance matrix with zero off-diagonal
         elements.
-    *args:
-        If the function accepts several arguments that are mutually independent, these
+    *args : tuple
+        If the function accepts several arguments that are mutually independent, it
         is possible to pass those values and covariance matrices pairwise, see examples.
-    **kwargs:
+    **kwargs : dict
         Extra arguments are passed to :func:`jacobi`.
 
     Returns
@@ -47,14 +56,14 @@ def propagate(
     y, ycov
         y is the result of fn(x).
         ycov is the propagated covariance matrix.
-        If ycov is a matrix, unless y is a number. In that case, ycov is also
+        ycov is a matrix, unless y is a number. In that case, ycov is also
         reduced to a number.
 
     Notes
     -----
     For callables `fn` which perform only element-wise computation, the jacobian is
     a diagonal matrix. This special case is detected and the computation optimised,
-    although can further speed up the computation by passing the argument
+    although one can further speed up the computation by passing the argument
     `diagonal=True`.
 
     In this special case, error propagation works correctly even if the output of `fn`
@@ -65,10 +74,9 @@ def propagate(
     General error propagation maps input vectors to output vectors.
 
     >>> def fn(x):
-    ...     return x ** 2 + 1
+    ...     return x**2 + 1
     >>> x = [1, 2]
-    >>> xcov = [[3, 1],
-    ...         [1, 4]]
+    >>> xcov = [[3, 1], [1, 4]]
     >>> y, ycov = propagate(fn, x, xcov)
 
     In the previous example, the function ``y = fn(x)`` treats all x values
@@ -82,7 +90,7 @@ def propagate(
     accepts several arguments, their uncertainties are treated as uncorrelated.
 
     >>> def fn(x, y):
-    ...    return x + y
+    ...     return x + y
     >>> x = 1
     >>> y = 2
     >>> xcov = 2
@@ -105,9 +113,10 @@ def propagate(
     """
     if args:
         if len(args) % 2 != 0:
-            raise ValueError("number of extra positional arguments must be even")
+            msg = "number of extra positional arguments must be even"
+            raise ValueError(msg)
 
-        args_a = [np.asarray(_) for _ in ((x, cov) + args)]
+        args_a = [np.asarray(arg) for arg in (x, cov, *args)]
         x_parts = args_a[::2]
         cov_parts = args_a[1::2]
         y_a = np.asarray(fn(*x_parts))
@@ -118,26 +127,34 @@ def propagate(
     try:
         y_a = np.asarray(fn(x_a))
     except ValueError as e:
-        raise ValueError(
-            "function return value cannot be converted into numpy array"
-        ) from e
+        msg = "function return value cannot be converted into numpy array"
+        raise ValueError(msg) from e
 
     # TODO lift this limitation
     if x_a.ndim > 1:
-        raise ValueError("x must have dimension 0 or 1")
+        msg = "x must have dimension 0 or 1"
+        raise ValueError(msg)
 
     # TODO lift this limitation
     if y_a.ndim > 1:
-        raise ValueError("function return value must have dimension 0 or 1")
+        msg = "function return value must have dimension 0 or 1"
+        raise ValueError(msg)
 
     # TODO lift this limitation
     if cov_a.ndim > 2:
-        raise ValueError("cov must have dimension 0, 1, or 2")
+        msg = "cov must have dimension 0, 1, or 2"
+        raise ValueError(msg)
 
     return _propagate(fn, y_a, x_a, cov_a, **kwargs)
 
 
-def _propagate(fn: Callable, y: np.ndarray, x: np.ndarray, xcov: np.ndarray, **kwargs):
+def _propagate(
+    fn: Callable[..., Any],
+    y: NDArray[Any],
+    x: NDArray[Any],
+    xcov: NDArray[Any],
+    **kwargs: Any,
+) -> tuple[NDArray[Any], NDArray[Any]]:
     _check_x_xcov_compatibility(x, xcov)
 
     jac = jacobi(fn, x, **kwargs)[0]
@@ -151,7 +168,7 @@ def _propagate(fn: Callable, y: np.ndarray, x: np.ndarray, xcov: np.ndarray, **k
         # contain NaN values.
         jac = _try_reduce_jacobian(jac)
     elif not diagonal:
-        jac.shape = (y.size, x.size)
+        jac = jac.reshape(y.size, x.size)
 
     ycov = _jac_cov_product(jac, xcov)
 
@@ -162,15 +179,15 @@ def _propagate(fn: Callable, y: np.ndarray, x: np.ndarray, xcov: np.ndarray, **k
 
 
 def _propagate_independent(
-    fn: Callable,
-    y: np.ndarray,
-    x_parts: List[np.ndarray],
-    xcov_parts: List[np.ndarray],
-    **kwargs,
-):
-    ycov: Union[float, np.ndarray] = 0
-
-    mask = kwargs.get("mask", None)
+    fn: Callable[..., Any],
+    y: NDArray[Any],
+    x_parts: list[NDArray[Any]],
+    xcov_parts: list[NDArray[Any]],
+    **kwargs: Any,
+) -> tuple[NDArray[Any], NDArray[Any]]:
+    ycov: NDArray[Any] | None = None
+
+    mask = kwargs.get("mask")
     mask_parts = []
     if mask is None:
         kwargs2 = kwargs
@@ -183,54 +200,69 @@ def _propagate_independent(
             elif np.shape(x) == np.shape(mask):
                 mask_parts.append(mask)
             else:
-                raise ValueError("mask shapes do not match arguments")
+                msg = "mask shapes do not match arguments"
+                raise ValueError(msg)
 
     for i, x in enumerate(x_parts):
-
-        def wrapped(x):
-            args = x_parts[:i] + [x] + x_parts[i + 1 :]
-            return fn(*args)
-
+        wrapped = _fix_other_arguments(fn, x_parts, i)
         xcov = xcov_parts[i]
 
         if mask_parts:
             kwargs2["mask"] = mask_parts[i]
 
         yc = _propagate(wrapped, y, x, xcov, **kwargs2)[1]
-        if np.ndim(ycov) == 2 and yc.ndim == 1:
-            for i, yci in enumerate(yc):
-                ycov[i, i] += yci  # type:ignore
+        if ycov is None:
+            ycov = np.array(yc)
+        elif ycov.ndim == 2 and yc.ndim == 1:
+            ycov[np.diag_indices_from(ycov)] += yc
+        elif ycov.ndim == 1 and yc.ndim == 2:
+            ycov = np.diag(ycov) + yc
         else:
-            ycov += yc
+            ycov = ycov + yc
+
+    if ycov is None:
+        msg = "at least one pair of value and covariance must be passed"
+        raise ValueError(msg)
 
     return y, ycov
 
 
-def _jac_cov_product(jac: np.ndarray, xcov: np.ndarray):
+def _fix_other_arguments(
+    fn: Callable[..., Any], x_parts: list[NDArray[Any]], i: int
+) -> Callable[[NDArray[Any]], Any]:
+    def wrapped(x: NDArray[Any]) -> Any:
+        return fn(*x_parts[:i], x, *x_parts[i + 1 :])
+
+    return wrapped
+
+
+def _jac_cov_product(jac: NDArray[Any], xcov: NDArray[Any]) -> NDArray[Any]:
     # if jac or xcov are 1D, they represent diagonal matrices
     if xcov.ndim == 2:
         if jac.ndim == 2:
-            return np.einsum("ij,kl,jl", jac, jac, xcov)
+            return np.asarray(np.einsum("ij,kl,jl", jac, jac, xcov))
         if jac.ndim == 1:
-            return np.einsum("i,j,ij -> ij", jac, jac, xcov)
-        return jac**2 * xcov
+            return np.asarray(np.einsum("i,j,ij -> ij", jac, jac, xcov))
+        return np.asarray(jac**2 * xcov)
     assert xcov.ndim < 2
     if jac.ndim == 2:
         if xcov.ndim == 1:
-            return np.einsum("ij,kj,j", jac, jac, xcov)
+            return np.asarray(np.einsum("ij,kj,j", jac, jac, xcov))
         assert xcov.ndim == 0  # xcov.ndim == 2 is already covered above
-        return np.einsum("ij,kj", jac, jac) * xcov
-    assert jac.ndim < 2 and xcov.ndim < 2
-    return xcov * jac**2
+        return np.asarray(np.einsum("ij,kj", jac, jac) * xcov)
+    assert jac.ndim < 2
+    assert xcov.ndim < 2
+    return np.asarray(xcov * jac**2)
 
 
-def _check_x_xcov_compatibility(x: np.ndarray, xcov: np.ndarray):
+def _check_x_xcov_compatibility(x: NDArray[Any], xcov: NDArray[Any]) -> None:
     if xcov.ndim > 0 and len(xcov) != (len(x) if x.ndim == 1 else 1):
         # this works for 1D and 2D xcov
-        raise ValueError("x and cov have incompatible shapes")
+        msg = "x and cov have incompatible shapes"
+        raise ValueError(msg)
 
 
-def _try_reduce_jacobian(jac: np.ndarray):
+def _try_reduce_jacobian(jac: NDArray[Any]) -> NDArray[Any]:
     if jac.ndim != 2 or jac.shape[0] != jac.shape[1]:
         return jac
     # if jacobian contains only off-diagonal elements
@@ -243,7 +275,7 @@ def _try_reduce_jacobian(jac: np.ndarray):
     return jac
 
 
-def _nodiag_view(a: np.ndarray):
+def _nodiag_view(a: NDArray[Any]) -> NDArray[Any]:
     # https://stackoverflow.com/a/43761941/ @Divakar
     m = a.shape[0]
     p, q = a.strides
diff --git a/src/jacobi/_typing.py b/src/jacobi/_typing.py
deleted file mode 100644
index 4beba56..0000000
--- a/src/jacobi/_typing.py
+++ /dev/null
@@ -1,13 +0,0 @@
-from typing import Sized, Generic, TypeVar
-import abc
-
-T = TypeVar("T")
-
-
-class Indexable(Sized, Generic[T]):
-    """Indexable type for mypy."""
-
-    @abc.abstractmethod
-    def __getitem__(self, idx: int) -> T:
-        """Get item at index idx."""
-        ...  # pragma: no cover
diff --git a/src/jacobi/_version.pyi b/src/jacobi/_version.pyi
new file mode 100644
index 0000000..55b0e01
--- /dev/null
+++ b/src/jacobi/_version.pyi
@@ -0,0 +1,15 @@
+__all__ = [
+    "__commit_id__",
+    "__version__",
+    "__version_tuple__",
+    "commit_id",
+    "version",
+    "version_tuple",
+]
+
+version: str
+__version__: str
+__version_tuple__: tuple[int | str, ...]
+version_tuple: tuple[int | str, ...]
+commit_id: str | None
+__commit_id__: str | None
diff --git a/src/jacobi/py.typed b/src/jacobi/py.typed
new file mode 100644
index 0000000..e69de29
diff --git a/tests/test_jacobi.py b/tests/test_jacobi.py
index 52d44b9..535daf3 100644
--- a/tests/test_jacobi.py
+++ b/tests/test_jacobi.py
@@ -1,6 +1,9 @@
-from numpy.testing import assert_allclose, assert_equal
+from __future__ import annotations
+
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose, assert_equal
+
 from jacobi import jacobi
 
 
@@ -64,9 +67,9 @@ def f7(x):
     "fn",
     [
         (lambda x: 1.0, lambda x: 0.0),
-        (lambda x: np.exp(x), lambda x: np.diagflat(np.exp(x))),
+        (np.exp, lambda x: np.diagflat(np.exp(x))),
         (lambda x: x**2, lambda x: np.diagflat(2 * x)),
-        (lambda x: np.ones_like(x), lambda x: np.diagflat(np.zeros_like(x))),
+        (np.ones_like, lambda x: np.diagflat(np.zeros_like(x))),
         (lambda x: x**-1, lambda x: np.diagflat(-(x**-2))),
         (lambda x: (x + 1) ** 0.5, lambda x: np.diagflat(0.5 * (x + 1) ** -0.5)),
         (f5, fd5),
@@ -138,7 +141,7 @@ def test_mask():
     assert_allclose(jac, jac_ref)
 
 
-@pytest.mark.parametrize("method", (-1, 0, 1))
+@pytest.mark.parametrize("method", [-1, 0, 1])
 def test_method(method):
     d = {}
     fp, fpe = jacobi(lambda x: x, 0, method=method, diagnostic=d)
@@ -173,7 +176,7 @@ def test_diagonal(x):
     assert_allclose(jac, fd1(np.atleast_1d(x), 3))
 
 
-@pytest.mark.parametrize("maxgrad", (0, 2, 5))
+@pytest.mark.parametrize("maxgrad", [0, 2, 5])
 def test_maxgrad(maxgrad):
     x = [1, 2, 3]
 
@@ -193,13 +196,13 @@ def test_rtol():
     assert_allclose(jac, jac_ref, rtol=0.03)
 
 
-@pytest.mark.parametrize("maxiter", (0, -1))
+@pytest.mark.parametrize("maxiter", [0, -1])
 def test_bad_maxiter(maxiter):
     with pytest.raises(ValueError):
         jacobi(f1, 1, 3, maxiter=maxiter)
 
 
-@pytest.mark.parametrize("maxgrad", (-1, -0.1))
+@pytest.mark.parametrize("maxgrad", [-1, -0.1])
 def test_bad_maxgrad(maxgrad):
     with pytest.raises(ValueError):
         jacobi(f1, 1, 3, maxgrad=maxgrad)
@@ -217,7 +220,7 @@ def test_bad_step_1(step):
         jacobi(f1, 1, 3, step=(0.25, step))
 
 
-@pytest.mark.parametrize("method", (-2, -3, 2, 3))
+@pytest.mark.parametrize("method", [-2, -3, 2, 3])
 def test_bad_method(method):
     with pytest.raises(ValueError, match="method"):
         jacobi(f1, 1, 3, method=method)
@@ -240,9 +243,8 @@ def test_on_nan():
     )
 
 
-@pytest.mark.parametrize("method", (None, -1, 0, 1))
-@pytest.mark.parametrize("fn", (lambda x: (x[0], 2 * x[1], x[1]), lambda x: x[0]))
-def test_non_array_arguments_and_return_value(method, fn):
+@pytest.mark.parametrize("method", [None, -1, 0, 1])
+def test_non_array_arguments_and_return_value(method):
     def fn(x):
         return [x[0], 2 * x[1], x[1] ** 2]
 
@@ -250,10 +252,69 @@ def fn(x):
     assert_allclose(j, [(1, 0), (0, 2), (0, 4)])
 
 
-@pytest.mark.parametrize("method", (None, -1, 0, 1))
+@pytest.mark.parametrize("method", [None, -1, 0, 1])
 @pytest.mark.parametrize(
-    "fn", (lambda x: [1, [1, 2]], lambda x: "s", lambda x: ("a", "b"))
+    "fn",
+    [lambda x: [1, [1, 2]], lambda x: "s", lambda x: ("a", "b")],
+    ids=["ragged", "str", "tuple_of_str"],
 )
 def test_bad_return_value_2(method, fn):
     with pytest.raises(ValueError, match="function return value cannot be converted"):
         jacobi(fn, (1, 2), method=method)
+
+
+@pytest.mark.parametrize("diagonal", [False, True])
+def test_empty_x(diagonal):
+    jac, jace = jacobi(np.exp, [], diagonal=diagonal)
+    assert jac.shape == ((0,) if diagonal else (0, 0))
+    assert jace.shape == jac.shape
+
+
+@pytest.mark.parametrize("diagonal", [False, True])
+def test_mask_selects_nothing(diagonal):
+    jac, jace = jacobi(np.exp, [1.0, 2.0], diagonal=diagonal, mask=[False, False])
+    assert_equal(jac, np.zeros((2,) if diagonal else (2, 2)))
+    assert_equal(jace, np.zeros_like(jac))
+
+
+@pytest.mark.parametrize("diagonal", [False, True])
+@pytest.mark.parametrize("x", [[], [1.0, 2.0]], ids=["empty", "nonempty"])
+@pytest.mark.parametrize(
+    "kwargs",
+    [{"maxiter": 0}, {"maxgrad": -1}, {"step": (0, 0.5)}, {"method": 5}],
+    ids=["maxiter", "maxgrad", "step", "method"],
+)
+def test_bad_arguments_on_all_paths(diagonal, x, kwargs):
+    # arguments are validated before any shortcut for empty or masked input
+    with pytest.raises(ValueError):
+        jacobi(np.exp, x, diagonal=diagonal, mask=np.zeros(len(x), bool), **kwargs)
+
+
+@pytest.mark.parametrize("diagonal", [False, True])
+def test_bad_mask(diagonal):
+    x = [1.0, 2.0, 3.0]
+    with pytest.raises(ValueError, match="boolean"):
+        jacobi(np.exp, x, diagonal=diagonal, mask=[0, 1, 1])
+    with pytest.raises(ValueError, match="shape"):
+        jacobi(np.exp, x, diagonal=diagonal, mask=[True, False])
+
+
+def test_large_x():
+    # regression test: fn used to be wrapped once per element of x, which hit
+    # the recursion limit for about 1000 elements
+    x = np.linspace(1, 2, 1000)
+    jac, jace = jacobi(np.exp, x)
+    assert_allclose(jac, np.diag(np.exp(x)))
+    assert_allclose(jace, 0, atol=1e-10)
+
+
+@pytest.mark.parametrize(
+    ("fn", "call_shape"),
+    [(lambda x: np.sum(x**2), (1, 2)), (lambda x: np.outer(x, x), (4, 2))],
+    ids=["scalar_output", "matrix_output"],
+)
+def test_diagnostic_output_ndim(fn, call_shape):
+    d = {}
+    jacobi(fn, [1.0, 2.0], diagnostic=d)
+    assert d["call"].shape == call_shape
+    assert np.all(d["call"] > 0)
diff --git a/tests/test_propagate.py b/tests/test_propagate.py
index 79336a6..5fed294 100644
--- a/tests/test_propagate.py
+++ b/tests/test_propagate.py
@@ -1,8 +1,10 @@
+from __future__ import annotations
+
 import numpy as np
-from numpy.testing import assert_allclose
-from jacobi import propagate, jacobi
 import pytest
-from numpy.testing import assert_equal
+from numpy.testing import assert_allclose, assert_equal
+
+from jacobi import jacobi, propagate
 
 
 def test_00():
@@ -77,9 +79,9 @@ def fn(x):
         assert_allclose(ycov, np.linalg.multi_dot([jac, xcov2, jac.T]))
 
 
-@pytest.mark.parametrize("ndim", (1, 2))
-@pytest.mark.parametrize("diagonal", (False, True))
-@pytest.mark.parametrize("len", (1, 2))
+@pytest.mark.parametrize("ndim", [1, 2])
+@pytest.mark.parametrize("diagonal", [False, True])
+@pytest.mark.parametrize("len", [1, 2])
 def test_cov_1d_2d(ndim, diagonal, len):
     def fn(x):
         return 2 * x
@@ -109,6 +111,8 @@ def fn1(x, y):
     z, zcov = propagate(fn1, x, xcov, y, ycov)
     assert z.ndim == 0
     assert zcov.ndim == 0
+    # a scalar result must behave like a float for callers
+    assert isinstance(zcov, float)
 
     def fn2(r):
         return fn1(r[0], r[1])
@@ -172,7 +176,7 @@ def fprime(x):
     x = 2
     xcov = 3
 
-    y, ycov = propagate(fn, x, xcov, diagonal=True)
+    _y, ycov = propagate(fn, x, xcov, diagonal=True)
 
     assert ycov.ndim == 0
     assert_allclose(ycov, fprime(x) ** 2 * xcov)
@@ -287,7 +291,7 @@ def f(a, b):
         return a * b
 
     mask = [False, True]
-    c, c_var = propagate(f, a, a_var, b, b_var, mask=mask)
+    _c, c_var = propagate(f, a, a_var, b, b_var, mask=mask)
 
     assert c_var[0] == 0
     assert c_var[1] > 0
@@ -303,15 +307,14 @@ def f(a, b):
         return np.outer(a, b).ravel()
 
     mask = [[False, True], [True, False, True]]
-    c1, c1_var = propagate(f, a, a_var, b, b_var, mask=mask)
-    c2, c2_var = propagate(f, a, a_var, b, b_var)
+    _c1, c1_var = propagate(f, a, a_var, b, b_var, mask=mask)
+    _c2, c2_var = propagate(f, a, a_var, b, b_var)
 
     assert np.sum(np.diag(c2_var) > np.diag(c1_var)) > 0
 
 
-@pytest.mark.parametrize("method", (None, -1, 0, 1))
-@pytest.mark.parametrize("fn", (lambda x: (x[0], 2 * x[1], x[1]), lambda x: x[0]))
-def test_non_array_arguments_and_return_value(method, fn):
+@pytest.mark.parametrize("method", [None, -1, 0, 1])
+def test_non_array_arguments_and_return_value(method):
     def fn(x):
         return [x[0], 2 * x[1], x[1] ** 3]
 
@@ -326,10 +329,61 @@ def fn(x):
     assert_allclose(ycov, ycov_ref)
 
 
-@pytest.mark.parametrize("method", (None, -1, 0, 1))
+@pytest.mark.parametrize("method", [None, -1, 0, 1])
 @pytest.mark.parametrize(
-    "fn", (lambda x: [1, [1, 2]], lambda x: "s", lambda x: ("a", "b"))
+    "fn",
+    [lambda x: [1, [1, 2]], lambda x: "s", lambda x: ("a", "b")],
+    ids=["ragged", "str", "tuple_of_str"],
 )
 def test_bad_return_value_2(method, fn):
     with pytest.raises(ValueError, match="function return value cannot be converted"):
         propagate(fn, (1, 2), ((1, 0), (0, 1)), method=method)
+
+
+def test_two_arguments_mixed_cov_dimensions():
+    # a diagonal (1D) covariance followed by a full (2D) one and vice versa
+    # must give the same result as the equivalent single-argument call
+    def fn1(x, y):
+        return x * y
+
+    x = [1.0, 2.0]
+    xcov = [0.1, 0.2]
+    y = [3.0, 4.0]
+    ycov = [[0.3, 0.05], [0.05, 0.4]]
+
+    z1, zcov1 = propagate(fn1, x, xcov, y, ycov)
+    z2, zcov2 = propagate(lambda y, x: fn1(x, y), y, ycov, x, xcov)
+
+    def fn2(r):
+        return fn1(r[:2], r[2:])
+
+    r = [*x, *y]
+    rcov = np.zeros((4, 4))
+    rcov[:2, :2] = np.diag(xcov)
+    rcov[2:, 2:] = ycov
+    z_ref, zcov_ref = propagate(fn2, r, rcov)
+
+    assert_allclose(z1, z_ref)
+    assert_allclose(z2, z_ref)
+    assert_allclose(zcov1, zcov_ref)
+    assert_allclose(zcov2, zcov_ref)
+
+
+@pytest.mark.parametrize("diagonal", [False, True])
+def test_fully_masked_argument(diagonal):
+    # argument a is exactly known, only b contributes to the uncertainty
+    def f(a, b):
+        return a * b
+
+    a = np.array([1.0, 2.0])
+    a_var = np.array([0.1, 0.1])
+    b = np.array([3.0, 4.0])
+    b_var = np.array([0.2, 0.3])
+
+    mask = [[False, False], [True, True]]
+    c, c_var = propagate(f, a, a_var, b, b_var, mask=mask, diagonal=diagonal)
+
+    assert_allclose(c, a * b)
+    if c_var.ndim == 2:
+        c_var = np.diag(c_var)
+    assert_allclose(c_var, a**2 * b_var)