diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 4c4b63b..2d1155d 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -24,6 +24,14 @@ jobs:
       with:
         submodules: recursive
 
+    - name: Cache GFortran install
+      if: ${{ contains(matrix.os, 'windows') }}
+      id: cache
+      uses: actions/cache@v2
+      with:
+        path: ./mingw-w64
+        key: gcc-${{ matrix.gcc }}-${{ matrix.os }}
+
     - name: Install GFortran (MacOS)
       if: ${{ contains(matrix.os, 'macos') }}
       run: |
@@ -40,7 +48,7 @@ jobs:
           --slave /usr/bin/gcov gcov /usr/bin/gcov-${{ matrix.gcc }}
 
     - name: Install GFortran (Windows)
-      if: ${{ contains(matrix.os, 'windows') }}
+      if: ${{ contains(matrix.os, 'windows') && steps.cache.outputs.cache-hit != 'true' }}
       run: |
         Invoke-WebRequest -Uri ${{ env.DOWNLOAD }} -OutFile mingw-w64.zip
         Expand-Archive mingw-w64.zip
@@ -70,6 +78,7 @@ jobs:
       if: ${{ matrix.build == 'meson' }}
       run: >-
         meson setup _build
+        --libdir=lib
         --prefix=${{ contains(matrix.os, 'windows') && '$pwd\_dist' || '$PWD/_dist' }}
 
     - name: Compile project (meson)
@@ -84,6 +93,146 @@ jobs:
       if: ${{ matrix.build == 'meson' }}
       run: meson install -C _build --no-rebuild
 
+    - name: Create package (Unix)
+      if: ${{ matrix.build == 'meson' && ! contains(matrix.os, 'windows') }}
+      run: |
+        tar cvf ${{ env.OUTPUT }} _dist
+        xz -T0 ${{ env.OUTPUT }}
+        echo "MINPACK_OUTPUT=${{ env.OUTPUT }}.xz" >> $GITHUB_ENV
+      env:
+        OUTPUT: minpack-${{ matrix.os }}.tar
+
+    - name: Create package (Windows)
+      if: ${{ matrix.build == 'meson' && contains(matrix.os, 'windows') }}
+      run: |
+        tar cvf ${{ env.OUTPUT }} _dist
+        xz -T0 ${{ env.OUTPUT }}
+        echo "MINPACK_OUTPUT=${{ env.OUTPUT }}.xz" | Out-File -FilePath $env:GITHUB_ENV -Encoding utf8 -Append
+      env:
+        OUTPUT: minpack-${{ matrix.os }}.tar
+
+    - name: Upload package
+      if: ${{ matrix.build == 'meson' }}
+      uses: actions/upload-artifact@v2
+      with:
+        name: ${{ env.MINPACK_OUTPUT }}
+        path: ${{ env.MINPACK_OUTPUT }}
+
+
+  Python:
+    needs:
+      - Build
+    runs-on: ${{ matrix.os }}
+    defaults:
+      run:
+        shell: ${{ contains(matrix.os, 'windows') && 'powershell' || 'bash -l {0}' }}
+    strategy:
+      fail-fast: false
+      matrix:
+        build: [meson]
+        os: [ubuntu-latest, macos-latest]
+        gcc: [10]
+        python: ['3.7', '3.8', '3.9']
+
+        # Additional test for setuptools build
+        include:
+        - build: setuptools
+          os: ubuntu-latest
+          gcc: 10
+          python: '3.9'
+
+    env:
+      FC: gfortran
+      CC: gcc
+      MINPACK_OUTPUT: minpack-${{ matrix.os }}.tar.xz
+
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v2
+
+    - name: Cache GFortran install
+      if: ${{ contains(matrix.os, 'windows') }}
+      id: cache
+      uses: actions/cache@v2
+      with:
+        path: ./mingw-w64
+        key: gcc-${{ matrix.gcc }}-${{ matrix.os }}
+
+    - name: Install dependencies
+      uses: mamba-org/provision-with-micromamba@main
+      with:
+        environment-file: config/ci/python-env.yaml
+        extra-specs: |
+          python=${{ matrix.python }}
+
+    - name: Install GFortran (MacOS)
+      if: ${{ contains(matrix.os, 'macos') }}
+      run: |
+        brew install gcc@${{ matrix.gcc }}
+        ln -s /usr/local/bin/gfortran-${{ matrix.gcc }} /usr/local/bin/gfortran
+        ln -s /usr/local/bin/gcc-${{ matrix.gcc }} /usr/local/bin/gcc
+
+    - name: Install GFortran (Linux)
+      if: ${{ contains(matrix.os, 'ubuntu') }}
+      run: |
+        sudo update-alternatives \
+          --install /usr/bin/gcc gcc /usr/bin/gcc-${{ matrix.gcc }} 100 \
+          --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-${{ matrix.gcc }} \
+          --slave /usr/bin/gcov gcov /usr/bin/gcov-${{ matrix.gcc }}
+
+    - name: Install GFortran (Windows)
+      if: ${{ contains(matrix.os, 'windows') && steps.cache.outputs.cache-hit != 'true' }}
+      run: |
+        Invoke-WebRequest -Uri ${{ env.DOWNLOAD }} -OutFile mingw-w64.zip
+        Expand-Archive mingw-w64.zip
+        echo "$pwd\mingw-w64\mingw64\bin" | Out-File -FilePath $env:GITHUB_PATH -Encoding utf8 -Append
+      shell: pwsh
+      env:
+        DOWNLOAD: "https://github.com/brechtsanders/winlibs_mingw/releases/download/10.3.0-12.0.0-9.0.0-r2/winlibs-x86_64-posix-seh-gcc-10.3.0-mingw-w64-9.0.0-r2.zip"
+
+    - name: Download package
+      uses: actions/download-artifact@v2
+      with:
+        name: ${{ env.MINPACK_OUTPUT }}
+
+    - name: Unpack package (Unix)
+      if: ${{ ! contains(matrix.os, 'windows') }}
+      run: |
+        tar xvf ${{ env.MINPACK_OUTPUT }}
+        echo "MINPACK_PREFIX=$PWD/_dist" >> $GITHUB_ENV
+
+    - name: Unpack package (Windows)
+      if: ${{ contains(matrix.os, 'windows') }}
+      run: |
+        tar xvf ${{ env.MINPACK_OUTPUT }}
+        echo "MINPACK_OUTPUT=${{ env.OUTPUT }}.xz" | Out-File -FilePath $env:GITHUB_ENV -Encoding utf8 -Append
+
+    - name: Install Python extension module (pip)
+      if: ${{ matrix.build == 'setuptools' }}
+      run: pip3 install . -vv
+      working-directory: python
+      env:
+        PKG_CONFIG_PATH: ${{ env.PKG_CONFIG_PATH }}:${{ env.MINPACK_PREFIX }}/lib/pkgconfig
+        LD_RUNPATH_SEARCH_PATH: ${{ env.MINPACK_PREFIX }}/lib
+
+    - name: Install Python extension module (meson)
+      if: ${{ matrix.build == 'meson' }}
+      run: |
+        set -ex
+        meson setup _build --prefix=$CONDA_PREFIX --libdir=lib
+        meson compile -C _build
+        meson install -C _build
+      working-directory: python
+      env:
+        PKG_CONFIG_PATH: ${{ env.PKG_CONFIG_PATH }}:${{ env.MINPACK_PREFIX }}/lib/pkgconfig
+
+    - name: Test Python API
+      run: pytest --doctest-modules --pyargs minpack --cov=minpack -vv
+      env:
+        LD_LIBRARY_PATH: ${{ env.LD_LIBRARY_PATH }}:${{ env.MINPACK_PREFIX }}/lib
+        DYLD_LIBRARY_PATH: ${{ env.DYLD_LIBRARY_PATH }}:${{ env.MINPACK_PREFIX }}/lib
+
+
   Docs:
     runs-on: ubuntu-latest
     defaults:
diff --git a/config/ci/python-env.yaml b/config/ci/python-env.yaml
new file mode 100644
index 0000000..be9e59b
--- /dev/null
+++ b/config/ci/python-env.yaml
@@ -0,0 +1,14 @@
+name: python
+channels:
+  - conda-forge
+dependencies:
+  - python
+  - pip
+  - pkgconfig
+  - pytest
+  - pytest-cov
+  - coverage
+  - cffi
+  - numpy
+  - meson
+  - ninja
diff --git a/include/minpack.h b/include/minpack.h
index f78a16b..91f817c 100644
--- a/include/minpack.h
+++ b/include/minpack.h
@@ -19,6 +19,16 @@ typedef void (*minpack_func)(
     int* /* iflag */,
     void* /* udata */);
 
+#ifdef MINPACK_CFFI
+extern "Python" void MINPACK_CALL
+func(
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    int* /* iflag */,
+    void* /* udata */);
+#endif
+
 /*
  * the purpose of hybrd is to find a zero of a system of
  * n nonlinear functions in n variables by a modification
@@ -84,6 +94,18 @@ typedef void (*minpack_fcn_hybrj)(
     int* /* iflag */,
     void* /* udata */);
 
+#ifdef MINPACK_CFFI
+extern "Python" void MINPACK_CALL
+fcn_hybrj(
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    int* /* iflag */,
+    void* /* udata */);
+#endif
+
 /*
  * the purpose of hybrj is to find a zero of a system of
  * n nonlinear functions in n variables by a modification
@@ -148,6 +170,19 @@ typedef void (*minpack_fcn_lmder)(
     int* /* iflag */,
     void* /* udata */);
 
+#ifdef MINPACK_CFFI
+extern "Python" void MINPACK_CALL
+fcn_lmder(
+    int /* m */,
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    int* /* iflag */,
+    void* /* udata */);
+#endif
+
 /*
  * the purpose of lmder is to minimize the sum of the squares of
  * m nonlinear functions in n variables by a modification of
@@ -213,6 +248,17 @@ typedef void (*minpack_func2)(
     int* /* iflag */,
     void* /* udata */);
 
+#ifdef MINPACK_CFFI
+extern "Python" void MINPACK_CALL
+func2(
+    int /* m */,
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    int* /* iflag */,
+    void* /* udata */);
+#endif
+
 /*
  * the purpose of lmdif is to minimize the sum of the squares of
  * m nonlinear functions in n variables by a modification of
@@ -279,6 +325,18 @@ typedef void (*minpack_fcn_lmstr)(
     int* /* iflag */,
     void* /* udata */);
 
+#ifdef MINPACK_CFFI
+extern "Python" void MINPACK_CALL
+fcn_lmstr(
+    int /* m */,
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    double* /* fjrow */,
+    int* /* iflag */,
+    void* /* udata */);
+#endif
+
 /*
  * the purpose of lmstr is to minimize the sum of the squares of
  * m nonlinear functions in n variables by a modification of
diff --git a/meson.build b/meson.build
index 52a600d..94eb4d0 100644
--- a/meson.build
+++ b/meson.build
@@ -8,6 +8,7 @@ project(
     'buildtype=debugoptimized',
   ],
 )
+has_cc = add_languages('c', required: get_option('python'), native: false)
 
 minpack_lib = library(
   meson.project_name(),
@@ -59,3 +60,7 @@ subdir('examples')
 
 # add the testsuite
 subdir('test')
+
+if get_option('python')
+  subdir('python/minpack')
+endif
diff --git a/meson_options.txt b/meson_options.txt
new file mode 100644
index 0000000..df0f2b1
--- /dev/null
+++ b/meson_options.txt
@@ -0,0 +1,12 @@
+option(
+  'python',
+  type: 'boolean',
+  value: false,
+  description: 'Build Python extension module',
+)
+option(
+  'python_version',
+  type: 'string',
+  value: 'python3',
+  description: 'Python version to link against.',
+)
diff --git a/python/.gitignore b/python/.gitignore
new file mode 100644
index 0000000..9842ebd
--- /dev/null
+++ b/python/.gitignore
@@ -0,0 +1,172 @@
+# Prerequisites
+*.d
+
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app
+
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+cover/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+.pybuilder/
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+#   For a library or package, you might want to ignore these files since the code is
+#   intended to run in multiple environments; otherwise, check them in:
+# .python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pytype static type analyzer
+.pytype/
+
+# Cython debug symbols
+cython_debug/
+
+# Directories
+/build*/
+/_*/
diff --git a/python/README.rst b/python/README.rst
new file mode 100644
index 0000000..dc9e38c
--- /dev/null
+++ b/python/README.rst
@@ -0,0 +1,85 @@
+Minpack Python bindings
+=======================
+
+Python bindings for Minpack.
+
+
+Building the extension module
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+To perform a build some version of ``minpack`` has to be available on your system and preferably findable by ``pkg-config``.
+Try to find a ``minpack`` installation you build against first with
+
+.. code:: sh
+
+   pkg-config --modversion minpack
+
+Adjust the ``PKG_CONFIG_PATH`` environment variable to include the correct directories to find the installation if necessary.
+
+
+Using pip
+^^^^^^^^^
+
+This project support installation with pip as an easy way to build the Python API.
+
+- C compiler to build the C-API and compile the extension module (the compiler name should be exported in the ``CC`` environment variable)
+- Python 3.6 or newer
+- The following Python packages are required additionally
+
+  - `cffi <https://cffi.readthedocs.io/>`_
+  - `numpy <https://numpy.org/>`_
+  - `pkgconfig <https://pypi.org/project/pkgconfig/>`_ (setup only)
+
+Make sure to have your C compiler set to the ``CC`` environment variable
+
+.. code:: sh
+
+   export CC=gcc
+
+Install the project with pip
+
+.. code:: sh
+
+   pip install .
+
+If you already have a ``minpack`` installation, *e.g.* from conda-forge, you can build the Python extension module directly without cloning this repository
+
+.. code:: sh
+
+   pip install "https://github.com/fortran-lang/minpack/archive/refs/heads/main.zip#egg=minpack&subdirectory=python"
+
+
+
+Using meson
+^^^^^^^^^^^
+
+This directory contains a separate meson build file to allow the out-of-tree build of the CFFI extension module.
+The out-of-tree build requires
+
+- C compiler to build the C-API and compile the extension module
+- `meson <https://mesonbuild.com>`_ version 0.53 or newer
+- a build-system backend, *i.e.* `ninja <https://ninja-build.org>`_ version 1.7 or newer
+- Python 3.6 or newer with the `CFFI <https://cffi.readthedocs.io/>`_ package installed
+
+Setup a build with
+
+.. code:: sh
+
+   meson setup _build -Dpython_version=$(which python3)
+
+The Python version can be used to select a different Python version, it defaults to ``'python3'``.
+Python 2 is not supported with this project, the Python version key is meant to select between several local Python 3 versions.
+
+Compile the project with
+
+.. code:: sh
+
+   meson compile -C _build
+
+The extension module is now available in ``_build/minpack/_libminpack.*.so``.
+You can install as usual with
+
+.. code:: sh
+
+   meson configure _build --prefix=/path/to/install
+   meson install -C _build
diff --git a/python/ffi-builder.py b/python/ffi-builder.py
new file mode 100644
index 0000000..d221f31
--- /dev/null
+++ b/python/ffi-builder.py
@@ -0,0 +1,83 @@
+"""
+FFI builder module for minpack for usage from meson and from setup.py.
+
+Since meson has the full knowledge about the build, it will handle
+the generation of the C definitions in the meson.build file rather
+than in the FFI builder. This allows to correctly keep track of
+dependencies and updates in the build process.
+
+For setup.py we have to do the preprocessing ourselves here, this
+requires us to use the C compiler to preprocess the header file
+of minpack because the CFFI C parser cannot handle certain C
+preprocessor constructs. Also, we cannot rely on an external build
+system fixing dependencies for us and therefore have to find those
+ourselves using pkg-config.
+"""
+
+import os
+import cffi
+
+library = "minpack"
+include_header = '#include "minpack.h"'
+prefix_var = "MINPACK_PREFIX"
+if prefix_var not in os.environ:
+    prefix_var = "CONDA_PREFIX"
+
+if __name__ == "__main__":
+    import sys
+
+    kwargs = dict(libraries=[library])
+
+    header_file = sys.argv[1]
+    module_name = sys.argv[2]
+
+    with open(header_file) as f:
+        cdefs = f.read()
+else:
+    import subprocess
+
+    try:
+        import pkgconfig
+
+        if not pkgconfig.exists(library):
+            raise ModuleNotFoundError("Unable to find pkg-config package 'minpack'")
+        if pkgconfig.installed(library, "< 2.0.0"):
+            raise Exception(
+                "Installed 'minpack' version is too old, 2.0.0 or newer is required"
+            )
+
+        kwargs = pkgconfig.parse(library)
+        cflags = pkgconfig.cflags(library).split()
+
+    except ModuleNotFoundError:
+        kwargs = dict(libraries=[library])
+        cflags = []
+        if prefix_var in os.environ:
+            prefix = os.environ[prefix_var]
+            kwargs.update(
+                include_dirs=[os.path.join(prefix, "include")],
+                library_dirs=[os.path.join(prefix, "lib")],
+                runtime_library_dirs=[os.path.join(prefix, "lib")],
+            )
+            cflags.append("-I" + os.path.join(prefix, "include"))
+
+    cc = os.environ["CC"] if "CC" in os.environ else "cc"
+
+    module_name = "minpack._libminpack"
+
+    p = subprocess.Popen(
+        [cc, *cflags, "-DMINPACK_CFFI=1", "-E", "-"],
+        stdin=subprocess.PIPE,
+        stdout=subprocess.PIPE,
+        stderr=subprocess.PIPE,
+    )
+    out, err = p.communicate(include_header.encode())
+
+    cdefs = out.decode()
+
+ffibuilder = cffi.FFI()
+ffibuilder.set_source(module_name, include_header, **kwargs)
+ffibuilder.cdef(cdefs)
+
+if __name__ == "__main__":
+    ffibuilder.distutils_extension(".")
diff --git a/python/meson.build b/python/meson.build
new file mode 100644
index 0000000..91d3eaa
--- /dev/null
+++ b/python/meson.build
@@ -0,0 +1,14 @@
+project(
+  'minpack',
+  'c',
+  meson_version: '>=0.53',
+  default_options: [
+    'buildtype=debugoptimized',
+  ],
+)
+install = true
+
+minpack_dep = dependency('minpack', version: '>=2.0.0')
+minpack_header = files('../include/minpack.h')
+
+subdir('minpack')
diff --git a/python/meson_options.txt b/python/meson_options.txt
new file mode 100644
index 0000000..7c91d21
--- /dev/null
+++ b/python/meson_options.txt
@@ -0,0 +1,6 @@
+option(
+  'python_version',
+  type: 'string',
+  value: 'python3',
+  description: 'Python version to link against.',
+)
diff --git a/python/minpack/__init__.py b/python/minpack/__init__.py
new file mode 100644
index 0000000..127d5a8
--- /dev/null
+++ b/python/minpack/__init__.py
@@ -0,0 +1,18 @@
+# -*- coding: utf-8 -*-
+"""
+Python bindings for the Minpack library.
+"""
+
+from .library import (
+    hybrd,
+    hybrj,
+    lmdif,
+    lmder,
+    lmstr,
+    hybrd1,
+    hybrj1,
+    lmdif1,
+    lmder1,
+    lmstr1,
+    chkder,
+)
diff --git a/python/minpack/exception.py b/python/minpack/exception.py
new file mode 100644
index 0000000..ad0d653
--- /dev/null
+++ b/python/minpack/exception.py
@@ -0,0 +1,101 @@
+# -*- coding: utf-8 -*-
+"""
+Possible exceptions for the `minpack` package mapping the info codes produced
+by the library to exceptions.
+"""
+
+from typing import Optional, Type
+
+
+class MinpackError(Exception):
+    """
+    Exception raised when Minpack returns an error.
+    """
+
+    pass
+
+
+class MinpackInputError(MinpackError):
+    """
+    Exception raised when Minpack input is invalid.
+    """
+
+    pass
+
+
+class MinpackMaxIterations(MinpackError):
+    """
+    The maximum number of calls to the objective function is reached.
+    """
+
+    pass
+
+
+class MinpackFunctionTolerance(MinpackError):
+    """
+    `ftol` is too small. No further reduction in the sum of squares is possible.
+    """
+
+    pass
+
+
+class MinpackSolutionTolerance(MinpackError):
+    """
+    `xtol` is too small. No further improvement in the approximate
+    solution x is possible.
+    """
+
+    pass
+
+
+class MinpackJacobianTolerance(MinpackError):
+    """
+    `gtol` is too small. `fvec` is orthogonal to the columns of
+    the Jacobian to machine precision.
+    """
+
+    pass
+
+
+class MinpackSlowProgress(MinpackError):
+    """
+    Iteration is not making good progress, as measured by the improvement
+    from the last ten iterations.
+    """
+
+    pass
+
+
+class MinpackSlowProgressJacobian(MinpackError):
+    """
+    Iteration is not making good progress, as measured by the improvement
+    from the last five jacobian evaluations.
+    """
+
+
+def info_hy(info: int) -> Optional[Type[MinpackError]]:
+    """
+    Get possible errors for `hybrd` and `hybrj` drivers.
+    """
+
+    return {
+        0: MinpackInputError,
+        2: MinpackMaxIterations,
+        3: MinpackFunctionTolerance,
+        4: MinpackSlowProgressJacobian,
+        5: MinpackSlowProgress,
+    }.get(info)
+
+
+def info_lm(info: int) -> Optional[Type[MinpackError]]:
+    """
+    Get possible errors for `lmdif`, `lmder`, and `lmstr` drivers.
+    """
+
+    return {
+        0: MinpackInputError,
+        5: MinpackMaxIterations,
+        6: MinpackFunctionTolerance,
+        7: MinpackSolutionTolerance,
+        8: MinpackJacobianTolerance,
+    }.get(info)
diff --git a/python/minpack/library.py b/python/minpack/library.py
new file mode 100644
index 0000000..72cd762
--- /dev/null
+++ b/python/minpack/library.py
@@ -0,0 +1,979 @@
+# -*- coding: utf-8 -*-
+"""
+Low-level Python bindings to the Minpack library.
+
+This module forwards the CFFI generated bindings to the Minpack library and provides
+a Pythonic interface to the C API.
+"""
+
+import numpy as np
+import functools
+import math
+from typing import Optional
+
+from .typing import (
+    CallableHybrd,
+    CallableHybrj,
+    CallableLmder,
+    CallableLmdif,
+    CallableLmstr,
+)
+from ._libminpack import ffi, lib
+from .exception import info_hy, info_lm
+
+
+class UserData:
+    """
+    Carry Python callable object through callback and propagate exceptions without
+    disturbing foreign runtime.
+    """
+
+    def __init__(self, fcn):
+        self.fcn = fcn
+        self.exception = None
+
+
+@ffi.def_extern()
+def func(n, x, fvec, iflag, data) -> None:
+    """
+    Entry point for callback from minpack_hybrd and minpack_hybrd1 library functions.
+    Restores type information for NDArray objects and calls the user-provided callback.
+    """
+    if iflag[0] <= 0:
+        return
+    handle: UserData = ffi.from_handle(data)
+    try:
+        handle.fcn(
+            np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real),
+            np.frombuffer(ffi.buffer(fvec, n * real.itemsize), dtype=real),
+        )
+    except BaseException as e:
+        iflag[0] = -1
+        handle.exception = e
+
+
+@ffi.def_extern()
+def fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag, data) -> None:
+    """
+    Entry point for callback from minpack_hybrj and minpack_hybrj1 library functions.
+    Restores type information for NDArray objects and calls the user-provided callback.
+    """
+    if iflag[0] <= 0:
+        return
+    handle: UserData = ffi.from_handle(data)
+    try:
+        fjac = np.frombuffer(ffi.buffer(fjac, ldfjac * n * real.itemsize), dtype=real)
+        handle.fcn(
+            np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real),
+            np.frombuffer(ffi.buffer(fvec, n * real.itemsize), dtype=real),
+            np.reshape(fjac, (n, ldfjac)),
+            iflag[0] == 2,
+        )
+    except BaseException as e:
+        iflag[0] = -1
+        handle.exception = e
+
+
+@ffi.def_extern()
+def fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag, data) -> None:
+    """
+    Entry point for callback from minpack_lmder and minpack_lmder1 library functions.
+    Restores type information for NDArray objects and calls the user-provided callback.
+    """
+    if iflag[0] <= 0:
+        return
+    handle: UserData = ffi.from_handle(data)
+    try:
+        fjac = np.frombuffer(ffi.buffer(fjac, ldfjac * n * real.itemsize), dtype=real)
+        handle.fcn(
+            np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real),
+            np.frombuffer(ffi.buffer(fvec, m * real.itemsize), dtype=real),
+            np.reshape(fjac, (n, ldfjac)),
+            iflag[0] == 2,
+        )
+    except BaseException as e:
+        iflag[0] = -1
+        handle.exception = e
+
+
+@ffi.def_extern()
+def func2(m, n, x, fvec, iflag, data) -> None:
+    """
+    Entry point for callback from minpack_lmdif and minpack_lmdif1 library functions.
+    Restores type information for NDArray objects and calls the user-provided callback.
+    """
+    if iflag[0] <= 0:
+        return
+    handle: UserData = ffi.from_handle(data)
+    try:
+        handle.fcn(
+            np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real),
+            np.frombuffer(ffi.buffer(fvec, m * real.itemsize), dtype=real),
+        )
+    except BaseException as e:
+        iflag[0] = -1
+        handle.exception = e
+
+
+@ffi.def_extern()
+def fcn_lmstr(m, n, x, fvec, fjrow, iflag, data) -> None:
+    """
+    Entry point for callback from minpack_lmstr and minpack_lmstr1 library functions.
+    Restores type information for NDArray objects and calls the user-provided callback.
+    """
+    if iflag[0] <= 0:
+        return
+    handle: UserData = ffi.from_handle(data)
+    try:
+        handle.fcn(
+            np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real),
+            np.frombuffer(ffi.buffer(fvec, m * real.itemsize), dtype=real),
+            np.frombuffer(ffi.buffer(fjrow, n * real.itemsize), dtype=real),
+            iflag[0] - 2 if iflag[0] > 1 else None,
+        )
+    except BaseException as e:
+        iflag[0] = -1
+        handle.exception = e
+
+
+def extern_python(dec):
+    """
+    Meta-decorator to attach a CFFI extern "Python" callback to a decorator
+    handling the Python side of the callback.
+    """
+
+    def layer(*args, **kwargs):
+        def wrapper(func):
+            return dec(func, *args, **kwargs)
+
+        return wrapper
+
+    return layer
+
+
+@extern_python
+def cffi_callback(func, callback):
+    """
+    Attach Python callback to a library function with extern "Python" callback.
+
+    This decorator wraps the user-provided Python callback in a `UserData` object
+    to carry it through the foreign runtime. It also propagates exceptions from
+    the user-provided Python callback back through the foreign runtime and re-raises.
+    """
+
+    @functools.wraps(func)
+    def entry_point(fcn, *args):
+        data = UserData(fcn)
+        handle = ffi.new_handle(data)
+        func(callback, *args, handle)
+        if data.exception is not None:
+            raise data.exception
+
+    return entry_point
+
+
+real = np.dtype("f8")
+minpack_hybrd1 = cffi_callback(lib.func)(lib.minpack_hybrd1)
+minpack_hybrd = cffi_callback(lib.func)(lib.minpack_hybrd)
+minpack_hybrj1 = cffi_callback(lib.fcn_hybrj)(lib.minpack_hybrj1)
+minpack_hybrj = cffi_callback(lib.fcn_hybrj)(lib.minpack_hybrj)
+minpack_lmder1 = cffi_callback(lib.fcn_lmder)(lib.minpack_lmder1)
+minpack_lmder = cffi_callback(lib.fcn_lmder)(lib.minpack_lmder)
+minpack_lmdif1 = cffi_callback(lib.func2)(lib.minpack_lmdif1)
+minpack_lmdif = cffi_callback(lib.func2)(lib.minpack_lmdif)
+minpack_lmstr1 = cffi_callback(lib.fcn_lmstr)(lib.minpack_lmstr1)
+minpack_lmstr = cffi_callback(lib.fcn_lmstr)(lib.minpack_lmstr)
+minpack_chkder = lib.minpack_chkder
+
+
+def hybrd1(
+    fcn: CallableHybrd,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    tol: float = math.sqrt(np.finfo(real).eps),
+) -> int:
+    """
+    Find a zero of a system of n nonlinear functions in n variables
+    by a modification of the Powell hybrid method.
+    This is done by using the more general nonlinear equation solver `hybrd`.
+    The user must provide a subroutine which calculates the functions.
+    The Jacobian is then calculated by a forward-difference approximation.
+
+    Parameters
+    ----------
+    func : callable ``f(x, fvec)``
+        A function that takes at least one (possibly vector) argument,
+        and returns a value of the same length.
+    x : ndarray
+        The starting estimate for the roots of ``func(x) = 0``.
+    fvec: ndarray
+        Function evaluated at the output
+    tol : float
+        The calculation will terminate if the relative error between two
+        consecutive iterates is at most `tol`.
+
+    Returns
+    -------
+    info : int
+        Set to 1 if a solution was found.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSlowProgressJacobian
+        When the Jacobian is not changing.
+    MinpackSlowProgress
+        When the function is not changing.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from minpack import hybrd
+    >>>
+    >>> def fcn(x, fvec) -> None:
+    ...     fvec[0] = x[0] * np.cos(x[1]) - 4
+    ...     fvec[1] = x[1] * x[0] - x[1] - 5
+    ...
+    >>> x = np.array(2 * [1.0])
+    >>> fvec = np.zeros(2, dtype=np.float64)
+    >>> hybrd1(fcn, x, fvec)
+    1
+    >>> x
+    array([6.50409711, 0.90841421])
+    >>> np.isclose(fvec, [0.0, 0.0])  # fvec should be almost 0.0.
+    array([ True,  True])
+    """
+    n = x.size
+    lwa = n * (3 * n + 13) // 2
+    wa = np.zeros(lwa, dtype=real)
+    info = ffi.new("int *")
+
+    minpack_hybrd1(
+        fcn,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        tol,
+        info,
+        ffi.cast("double*", wa.ctypes.data),
+        lwa,
+    )
+    ex = info_hy(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def hybrd(
+    fcn: CallableHybrd,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    xtol: float = math.sqrt(np.finfo(real).eps),
+    *,
+    maxfev: Optional[int] = None,
+    ml: Optional[int] = None,
+    mu: Optional[int] = None,
+    epsfcn: float = 0.0,
+    diag: Optional[np.ndarray] = None,
+    mode: int = 2,
+    factor: float = 100.0,
+    nprint: int = 0,
+    fjac: Optional[np.ndarray] = None,
+    r: Optional[np.ndarray] = None,
+    qtf: Optional[np.ndarray] = None,
+) -> int:
+    """
+    Find a zero of a system of n nonlinear functions in n variables
+    by a modification of the Powell hybrid method.
+    The user must provide a subroutine which calculates the functions.
+    The Jacobian is then calculated by a forward-difference approximation.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSlowProgressJacobian
+        When the Jacobian is not changing.
+    MinpackSlowProgress
+        When the function is not changing.
+    """
+    n = x.size
+    info = ffi.new("int *")
+    nfev = ffi.new("int *")
+    if maxfev is None:
+        maxfev = 200 * (n + 1)
+    if ml is None:
+        ml = n - 1
+    if mu is None:
+        mu = n - 1
+    if diag is None:
+        diag = np.ones(n, dtype=real)
+    if fjac is None:
+        fjac = np.zeros((n, n), dtype=real)
+    if r is None:
+        r = np.zeros(n * (n + 1) // 2, dtype=real)
+    if qtf is None:
+        qtf = np.zeros(n, dtype=real)
+    wa1 = np.zeros(n, dtype=real)
+    wa2 = np.zeros(n, dtype=real)
+    wa3 = np.zeros(n, dtype=real)
+    wa4 = np.zeros(n, dtype=real)
+
+    minpack_hybrd(
+        fcn,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        xtol,
+        maxfev,
+        ml,
+        mu,
+        epsfcn,
+        ffi.cast("double*", diag.ctypes.data),
+        mode,
+        factor,
+        nprint,
+        info,
+        nfev,
+        ffi.cast("double*", fjac.ctypes.data),
+        n,
+        ffi.cast("double*", r.ctypes.data),
+        r.size,
+        ffi.cast("double*", qtf.ctypes.data),
+        ffi.cast("double*", wa1.ctypes.data),
+        ffi.cast("double*", wa2.ctypes.data),
+        ffi.cast("double*", wa3.ctypes.data),
+        ffi.cast("double*", wa4.ctypes.data),
+    )
+    ex = info_hy(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def hybrj1(
+    fcn: CallableHybrj,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    fjac: np.ndarray,
+    tol: float = math.sqrt(np.finfo(real).eps),
+) -> int:
+    """
+    Find a zero of a system of n nonlinear functions in n variables
+    by a modification of the Powell hybrid method.
+    This is done by using the more general nonlinear equation solver `hybrj`.
+    The user must provide a subroutine which calculates the functions
+    and the Jacobian.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSlowProgressJacobian
+        When the Jacobian is not changing.
+    MinpackSlowProgress
+        When the function is not changing.
+    """
+    n = x.size
+    lwa = (n * (n + 13)) // 2
+    wa = np.zeros(lwa, dtype=real)
+    info = ffi.new("int *")
+
+    minpack_hybrj1(
+        fcn,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        ffi.cast("double*", fjac.ctypes.data),
+        n,
+        tol,
+        info,
+        ffi.cast("double*", wa.ctypes.data),
+        lwa,
+    )
+    ex = info_hy(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def hybrj(
+    fcn: CallableHybrj,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    fjac: np.ndarray,
+    xtol: float = math.sqrt(np.finfo(real).eps),
+    *,
+    maxfev: Optional[int] = None,
+    diag: Optional[np.ndarray] = None,
+    mode: int = 2,
+    factor: float = 100.0,
+    nprint: int = 0,
+    r: Optional[np.ndarray] = None,
+    qtf: Optional[np.ndarray] = None,
+) -> int:
+    """
+    Find a zero of a system of n nonlinear functions in n variables
+    by a modification of the Powell hybrid method.
+    The user must provide a subroutine which calculates the functions
+    and the Jacobian.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSlowProgressJacobian
+        When the Jacobian is not changing.
+    MinpackSlowProgress
+        When the function is not changing.
+    """
+    n = x.size
+    info = ffi.new("int *")
+    nfev = ffi.new("int *")
+    njev = ffi.new("int *")
+    if maxfev is None:
+        maxfev = 200 * (n + 1)
+    if diag is None:
+        diag = np.ones(n, dtype=real)
+    if fjac is None:
+        fjac = np.zeros((n, n), dtype=real)
+    if r is None:
+        r = np.zeros(n * (n + 1) // 2, dtype=real)
+    if qtf is None:
+        qtf = np.zeros(n, dtype=real)
+    wa1 = np.zeros(n, dtype=real)
+    wa2 = np.zeros(n, dtype=real)
+    wa3 = np.zeros(n, dtype=real)
+    wa4 = np.zeros(n, dtype=real)
+
+    minpack_hybrj(
+        fcn,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        ffi.cast("double*", fjac.ctypes.data),
+        n,
+        xtol,
+        maxfev,
+        ffi.cast("double*", diag.ctypes.data),
+        mode,
+        factor,
+        nprint,
+        info,
+        nfev,
+        njev,
+        ffi.cast("double*", r.ctypes.data),
+        r.size,
+        ffi.cast("double*", qtf.ctypes.data),
+        ffi.cast("double*", wa1.ctypes.data),
+        ffi.cast("double*", wa2.ctypes.data),
+        ffi.cast("double*", wa3.ctypes.data),
+        ffi.cast("double*", wa4.ctypes.data),
+    )
+    ex = info_hy(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def lmder1(
+    fcn: CallableLmder,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    fjac: np.ndarray,
+    tol: float = math.sqrt(np.finfo(real).eps),
+) -> int:
+    """
+    Minimize the sum of the squares of m nonlinear functions in n variables
+    by a modification of the Levenberg-Marquardt algorithm.
+    This is done by using the more general least-squares solver `lmder`.
+    The user must provide a subroutine which calculates the functions and the Jacobian.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSolutionTolerance
+        When no further improvement in the approximate solution is possible.
+    MinpackJacobianTolerance
+        The solution is orthogonal to the jacobian.
+    """
+    n = x.size
+    m = fvec.size
+    lwa = 5 * n + m
+    wa = np.zeros(lwa, dtype=real)
+    ipvt = np.zeros(n, dtype=np.int32)
+    info = ffi.new("int *")
+
+    minpack_lmder1(
+        fcn,
+        m,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        ffi.cast("double*", fjac.ctypes.data),
+        m,
+        tol,
+        info,
+        ffi.cast("int*", ipvt.ctypes.data),
+        ffi.cast("double*", wa.ctypes.data),
+        lwa,
+    )
+    ex = info_lm(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def lmder(
+    fcn: CallableLmder,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    fjac: np.ndarray,
+    ftol: float = math.sqrt(np.finfo(real).eps),
+    xtol: float = math.sqrt(np.finfo(real).eps),
+    *,
+    gtol: float = 0.0,
+    maxfev: Optional[int] = None,
+    diag: Optional[np.ndarray] = None,
+    mode: int = 1,
+    factor=100.0,
+    nprint=0,
+    ipvt: Optional[np.ndarray] = None,
+    qtf: Optional[np.ndarray] = None,
+) -> int:
+    """
+    Minimize the sum of the squares of m nonlinear functions in n variables
+    by a modification of the Levenberg-Marquardt algorithm.
+    The user must provide a subroutine which calculates the functions and the Jacobian.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSolutionTolerance
+        When no further improvement in the approximate solution is possible.
+    MinpackJacobianTolerance
+        The solution is orthogonal to the jacobian.
+    """
+    n = x.size
+    m = fvec.size
+    info = ffi.new("int *")
+    nfev = ffi.new("int *")
+    njev = ffi.new("int *")
+    if diag is None:
+        diag = np.ones(n, dtype=real)
+    if maxfev is None:
+        maxfev = 100 * (n + 1)
+    if ipvt is None:
+        ipvt = np.zeros(n, dtype=np.int32)
+    if qtf is None:
+        qtf = np.zeros(n, dtype=real)
+    wa1 = np.zeros(n, dtype=real)
+    wa2 = np.zeros(n, dtype=real)
+    wa3 = np.zeros(n, dtype=real)
+    wa4 = np.zeros(m, dtype=real)
+
+    minpack_lmder(
+        fcn,
+        m,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        ffi.cast("double*", fjac.ctypes.data),
+        m,
+        ftol,
+        xtol,
+        gtol,
+        maxfev,
+        ffi.cast("double*", diag.ctypes.data),
+        mode,
+        factor,
+        nprint,
+        info,
+        nfev,
+        njev,
+        ffi.cast("int*", ipvt.ctypes.data),
+        ffi.cast("double*", qtf.ctypes.data),
+        ffi.cast("double*", wa1.ctypes.data),
+        ffi.cast("double*", wa2.ctypes.data),
+        ffi.cast("double*", wa3.ctypes.data),
+        ffi.cast("double*", wa4.ctypes.data),
+    )
+    ex = info_lm(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def lmdif1(
+    fcn: CallableLmdif,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    tol: float = math.sqrt(np.finfo(real).eps),
+) -> int:
+    """
+    Minimize the sum of the squares of m nonlinear functions in n variables
+    by a modification of the Levenberg-Marquardt algorithm.
+    This is done by using the more general least-squares solver `lmdif`.
+    The user must provide a subroutine which calculates the functions.
+    The jacobian is then calculated by a forward-difference approximation.
+
+    Parameters
+    ----------
+    func : callable ``f(x, fvec)``
+        Should take at least one (possibly length n vector) argument and
+        compute m floating point numbers in fvec. It must not return NaNs or
+        fitting might fail. m must be greater than or equal to n.
+    x : ndarray
+        The starting estimate for the minimization.
+    fvec : ndarray
+        The function evaluated at the output.
+    tol : float, optional
+        Relative error desired in the sum of squares and the approximate solution.
+
+    Returns
+    -------
+    info : int
+        An integer flag. If it is equal to 1, 2, 3 or 4, the solution was
+        found. Otherwise, the solution was not found.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSolutionTolerance
+        When no further improvement in the approximate solution is possible.
+    MinpackJacobianTolerance
+        The solution is orthogonal to the jacobian.
+
+    Example
+    -------
+    >>> import numpy as np
+    >>> from minpack import lmdif1
+    >>>
+    >>> def func(x, fvec):
+    ...     fvec[:] = 2*(x-3)**2+1
+    ...
+    >>> x = np.array(0.0)
+    >>> fvec = np.zeros(1, dtype=np.float64)
+    >>> lmdif1(func, x, fvec)
+    1
+    >>> x
+    array(2.99999999)
+    """
+    n = x.size
+    m = fvec.size
+    lwa = m * n + 5 * n + m
+    wa = np.zeros(lwa, dtype=real)
+    ipvt = np.zeros(n, dtype=np.int32)
+    info = ffi.new("int *")
+
+    minpack_lmdif1(
+        fcn,
+        m,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        tol,
+        info,
+        ffi.cast("int*", ipvt.ctypes.data),
+        ffi.cast("double*", wa.ctypes.data),
+        lwa,
+    )
+    ex = info_lm(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def lmdif(
+    fcn: CallableLmdif,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    ftol: float = math.sqrt(np.finfo(real).eps),
+    xtol: float = math.sqrt(np.finfo(real).eps),
+    *,
+    gtol: float = 0.0,
+    maxfev: Optional[int] = None,
+    epsfcn: float = 0.0,
+    diag: Optional[np.ndarray] = None,
+    mode=1,
+    factor: float = 100.0,
+    nprint: int = 0,
+    fjac: Optional[np.ndarray] = None,
+    ipvt: Optional[np.ndarray] = None,
+    qtf: Optional[np.ndarray] = None,
+) -> int:
+    """
+    Minimize the sum of the squares of m nonlinear functions in n variables
+    by a modification of the Levenberg-Marquardt algorithm.
+    The user must provide a subroutine which calculates the functions.
+    The jacobian is then calculated by a forward-difference approximation.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSolutionTolerance
+        When no further improvement in the approximate solution is possible.
+    MinpackJacobianTolerance
+        The solution is orthogonal to the jacobian.
+    """
+    n = x.size
+    m = fvec.size
+    info = ffi.new("int *")
+    nfev = ffi.new("int *")
+    if maxfev is None:
+        maxfev = 200 * (n + 1)
+    if diag is None:
+        diag = np.ones(n, dtype=real)
+    if fjac is None:
+        fjac = np.zeros((n, m), dtype=real)
+    if ipvt is None:
+        ipvt = np.zeros(n, dtype=np.int32)
+    if qtf is None:
+        qtf = np.zeros(n, dtype=real)
+    wa1 = np.zeros(n, dtype=real)
+    wa2 = np.zeros(n, dtype=real)
+    wa3 = np.zeros(n, dtype=real)
+    wa4 = np.zeros(m, dtype=real)
+
+    minpack_lmdif(
+        fcn,
+        m,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        ftol,
+        xtol,
+        gtol,
+        maxfev,
+        epsfcn,
+        ffi.cast("double*", diag.ctypes.data),
+        mode,
+        factor,
+        nprint,
+        info,
+        nfev,
+        ffi.cast("double*", fjac.ctypes.data),
+        m,
+        ffi.cast("int*", ipvt.ctypes.data),
+        ffi.cast("double*", qtf.ctypes.data),
+        ffi.cast("double*", wa1.ctypes.data),
+        ffi.cast("double*", wa2.ctypes.data),
+        ffi.cast("double*", wa3.ctypes.data),
+        ffi.cast("double*", wa4.ctypes.data),
+    )
+    ex = info_lm(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def lmstr1(
+    fcn: CallableLmstr,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    fjac: np.ndarray,
+    tol: float = math.sqrt(np.finfo(real).eps),
+) -> int:
+    """
+    Minimize the sum of the squares of m nonlinear functions in n variables by
+    a modification of the Levenberg-Marquardt algorithm which uses minimal storage.
+    This is done by using the more general least-squares solver `lmstr`.
+    The user must provide a subroutine which calculates the functions and
+    the rows of the Jacobian.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSolutionTolerance
+        When no further improvement in the approximate solution is possible.
+    MinpackJacobianTolerance
+        The solution is orthogonal to the jacobian.
+    """
+    n = x.size
+    m = fvec.size
+    lwa = m * n + 5 * n + m
+    wa = np.zeros(lwa, dtype=real)
+    ipvt = np.zeros(n, dtype=np.int32)
+    info = ffi.new("int *")
+
+    minpack_lmstr1(
+        fcn,
+        m,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        ffi.cast("double*", fjac.ctypes.data),
+        n,
+        tol,
+        info,
+        ffi.cast("int*", ipvt.ctypes.data),
+        ffi.cast("double*", wa.ctypes.data),
+        lwa,
+    )
+    ex = info_lm(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def lmstr(
+    fcn: CallableLmstr,
+    x: np.ndarray,
+    fvec: np.ndarray,
+    fjac: np.ndarray,
+    ftol: float = math.sqrt(np.finfo(real).eps),
+    xtol: float = math.sqrt(np.finfo(real).eps),
+    *,
+    gtol: float = 0.0,
+    maxfev: Optional[int] = None,
+    diag: Optional[np.ndarray] = None,
+    mode: int = 1,
+    factor: float = 100.0,
+    nprint=0,
+    ipvt: Optional[np.ndarray] = None,
+    qtf: Optional[np.ndarray] = None,
+) -> int:
+    """
+    Minimize the sum of the squares of m nonlinear functions in n variables by
+    a modification of the Levenberg-Marquardt algorithm which uses minimal storage.
+    The user must provide a subroutine which calculates the functions and
+    the rows of the Jacobian.
+
+    Raises
+    ------
+    MinpackInputError
+        In case of invalid input parameters.
+    MinpackMaxIterations
+        When the maximum number of iterations is exceeded.
+    MinpackFunctionTolerance
+        When the function tolerance cannot be satisfied.
+    MinpackSolutionTolerance
+        When no further improvement in the approximate solution is possible.
+    MinpackJacobianTolerance
+        The solution is orthogonal to the jacobian.
+    """
+    n = x.size
+    m = fvec.size
+    info = ffi.new("int *")
+    nfev = ffi.new("int *")
+    njev = ffi.new("int *")
+    if maxfev is None:
+        maxfev = 100 * (n + 1)
+    if diag is None:
+        diag = np.ones(n, dtype=real)
+    if ipvt is None:
+        ipvt = np.zeros(n, dtype=np.int32)
+    if qtf is None:
+        qtf = np.zeros(n, dtype=real)
+    wa1 = np.zeros(n, dtype=real)
+    wa2 = np.zeros(n, dtype=real)
+    wa3 = np.zeros(n, dtype=real)
+    wa4 = np.zeros(m, dtype=real)
+
+    minpack_lmstr(
+        fcn,
+        m,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        ffi.cast("double*", fjac.ctypes.data),
+        n,
+        ftol,
+        xtol,
+        gtol,
+        maxfev,
+        ffi.cast("double*", diag.ctypes.data),
+        mode,
+        factor,
+        nprint,
+        info,
+        nfev,
+        njev,
+        ffi.cast("int*", ipvt.ctypes.data),
+        ffi.cast("double*", qtf.ctypes.data),
+        ffi.cast("double*", wa1.ctypes.data),
+        ffi.cast("double*", wa2.ctypes.data),
+        ffi.cast("double*", wa3.ctypes.data),
+        ffi.cast("double*", wa4.ctypes.data),
+    )
+    ex = info_lm(info[0])
+    if ex is not None:
+        raise ex(ex.__doc__)
+    return info[0]
+
+
+def chkder(x, fvec, fjac, xp, fvecp, check, error):
+    """
+    This subroutine checks the gradients of m nonlinear functions
+    in n variables, evaluated at a point x, for consistency with
+    the functions themselves.
+
+    The subroutine does not perform reliably if cancellation or
+    rounding errors cause a severe loss of significance in the
+    evaluation of a function. Therefore, none of the components
+    of x should be unusually small (in particular, zero) or any
+    other value which may cause loss of significance.
+    """
+    if not fvec.size == fjac.shape[-1] == fvecp.size == error.size:
+        raise ValueError("fvec, fjac, fvecp, error must have the same size")
+    if not x.size == fjac.shape[0] == xp.size:
+        raise ValueError("x, fjac, xp must have the same size")
+
+    m = fvec.size
+    n = x.size
+    ldfjac = fjac.shape[-1]
+
+    minpack_chkder(
+        m,
+        n,
+        ffi.cast("double*", x.ctypes.data),
+        ffi.cast("double*", fvec.ctypes.data),
+        ffi.cast("double*", fjac.ctypes.data),
+        ldfjac,
+        ffi.cast("double*", xp.ctypes.data),
+        ffi.cast("double*", fvecp.ctypes.data),
+        2 if check else 1,
+        ffi.cast("double*", error.ctypes.data),
+    )
diff --git a/python/minpack/meson.build b/python/minpack/meson.build
new file mode 100644
index 0000000..812ba00
--- /dev/null
+++ b/python/minpack/meson.build
@@ -0,0 +1,45 @@
+cc = meson.get_compiler('c')
+ext_module = '_libminpack'
+
+pymod = import('python')
+python = pymod.find_installation(
+  get_option('python_version'),
+  modules: [
+    'cffi',
+  ],
+)
+python_dep = python.dependency(required: true)
+
+# Python's CFFI is horrible in working with preprocessor statements,
+# therefore, we have to preprocess the header before passing it to the ffibuilder
+minpack_pp = configure_file(
+  command: [cc, '-DMINPACK_CFFI=1', '-E', '@INPUT@'],
+  input: minpack_header[0],
+  output: '@0@.h'.format(ext_module),
+  capture: true,
+)
+
+# This is the actual out-of-line API processing of the ffibuilder
+minpack_cffi_srcs = configure_file(
+  command: [python, files('..'/'ffi-builder.py'), '@INPUT@', '@BASENAME@'],
+  input: minpack_pp,
+  output: '@BASENAME@.c',
+)
+
+# Actual generation of the Python extension
+minpack_pyext = python.extension_module(
+  ext_module,
+  minpack_cffi_srcs,
+  dependencies: [minpack_dep, python_dep],
+  install: true,
+  subdir: 'minpack',
+)
+
+python.install_sources(
+  '__init__.py',
+  'exception.py',
+  'library.py',
+  'typing.py',
+  'test_library.py',
+  subdir: 'minpack',
+)
diff --git a/python/minpack/test_library.py b/python/minpack/test_library.py
new file mode 100644
index 0000000..6d729ee
--- /dev/null
+++ b/python/minpack/test_library.py
@@ -0,0 +1,272 @@
+# -*- coding: utf-8 -*-
+
+import pytest
+import minpack.library
+import numpy as np
+from math import sqrt
+
+
+@pytest.mark.parametrize("driver", [minpack.library.hybrd1, minpack.library.hybrd])
+def test_hybrd(driver):
+    def fcn(x, fvec) -> None:
+        for k in range(x.size):
+            tmp = (3.0 - 2.0 * x[k]) * x[k]
+            tmp1 = x[k - 1] if k > 0 else 0.0
+            tmp2 = x[k + 1] if k < len(x) - 1 else 0.0
+            fvec[k] = tmp - tmp1 - 2.0 * tmp2 + 1.0
+
+    x = np.array(9 * [-1.0])
+    fvec = np.zeros(9, dtype=np.float64)
+    tol = sqrt(np.finfo(np.float64).eps)
+
+    assert driver(fcn, x, fvec, tol) == 1
+
+    assert pytest.approx(x, abs=10 * tol) == [
+        -0.5706545,
+        -0.6816283,
+        -0.7017325,
+        -0.7042129,
+        -0.7013690,
+        -0.6918656,
+        -0.6657920,
+        -0.5960342,
+        -0.4164121,
+    ]
+
+
+@pytest.mark.parametrize("driver", [minpack.library.hybrd1, minpack.library.hybrd])
+def test_hybrd_exception(driver):
+    class DummyException(Exception):
+        ...
+
+    def fcn(x, fvec) -> None:
+        raise DummyException()
+
+    x = np.array(9 * [-1.0])
+    fvec = np.zeros(9, dtype=np.float64)
+
+    with pytest.raises(DummyException):
+        driver(fcn, x, fvec)
+
+
+@pytest.mark.parametrize("driver", [minpack.library.hybrj1, minpack.library.hybrj])
+def test_hybrj(driver):
+    def fcn(x, fvec, fjac, jacobian: bool) -> None:
+
+        if jacobian:
+            for k in range(x.size):
+                for j in range(x.size):
+                    fjac[k, j] = 0.0
+                fjac[k, k] = 3.0 - 4.0 * x[k]
+                if k > 0:
+                    fjac[k, k - 1] = -1.0
+                if k < x.size - 1:
+                    fjac[k, k + 1] = -2.0
+        else:
+            for k in range(x.size):
+                tmp = (3.0 - 2.0 * x[k]) * x[k]
+                tmp1 = x[k - 1] if k > 0 else 0.0
+                tmp2 = x[k + 1] if k < len(x) - 1 else 0.0
+                fvec[k] = tmp - tmp1 - 2.0 * tmp2 + 1.0
+
+    x = np.array(9 * [-1.0])
+    fvec = np.zeros(9, dtype=np.float64)
+    fjac = np.zeros((9, 9), dtype=np.float64)
+    tol = sqrt(np.finfo(np.float64).eps)
+
+    assert driver(fcn, x, fvec, fjac, tol) == 1
+
+    assert pytest.approx(x, abs=10 * tol) == [
+        -0.5706545,
+        -0.6816283,
+        -0.7017325,
+        -0.7042129,
+        -0.7013690,
+        -0.6918656,
+        -0.6657920,
+        -0.5960342,
+        -0.4164121,
+    ]
+
+
+@pytest.mark.parametrize("driver", [minpack.library.hybrj1, minpack.library.hybrj])
+def test_hybrj_exception(driver):
+    class DummyException(Exception):
+        ...
+
+    def fcn(x, fvec, fjac, jacobian) -> None:
+        raise DummyException()
+
+    x = np.array(9 * [-1.0])
+    fvec = np.zeros(9, dtype=np.float64)
+    fjac = np.zeros((9, 9), dtype=np.float64)
+
+    with pytest.raises(DummyException):
+        driver(fcn, x, fvec, fjac)
+
+
+@pytest.mark.parametrize("driver", [minpack.library.lmder1, minpack.library.lmder])
+def test_lmder(driver):
+    y = np.array(
+        [
+            1.4e-1,
+            1.8e-1,
+            2.2e-1,
+            2.5e-1,
+            2.9e-1,
+            3.2e-1,
+            3.5e-1,
+            3.9e-1,
+            3.7e-1,
+            5.8e-1,
+            7.3e-1,
+            9.6e-1,
+            1.34e0,
+            2.1e0,
+            4.39e0,
+        ]
+    )
+
+    def fcn(x, fvec, fjac, jacobian: bool) -> None:
+        if jacobian:
+            for i in range(fvec.size):
+                tmp1, tmp2 = i + 1, 16 - i - 1
+                tmp3 = tmp2 if i >= 8 else tmp1
+                tmp4 = (x[1] * tmp2 + x[2] * tmp3) ** 2
+                fjac[0, i] = -1.0
+                fjac[1, i] = tmp1 * tmp2 / tmp4
+                fjac[2, i] = tmp1 * tmp3 / tmp4
+        else:
+            for i in range(fvec.size):
+                tmp1, tmp2 = i + 1, 16 - i - 1
+                tmp3 = tmp2 if i >= 8 else tmp1
+                fvec[i] = y[i] - (x[0] + tmp1 / (x[1] * tmp2 + x[2] * tmp3))
+
+    x = np.array([1.0, 1.0, 1.0])
+    fvec = np.zeros(15, dtype=np.float64)
+    fjac = np.zeros((3, 15), dtype=np.float64)
+    tol = sqrt(np.finfo(np.float64).eps)
+
+    xp = np.zeros(3, dtype=np.float64)
+    fvecp = np.zeros(15, dtype=np.float64)
+    err = np.zeros(15, dtype=np.float64)
+    minpack.library.chkder(x, fvec, fjac, xp, fvecp, False, err)
+    fcn(x, fvec, fjac, False)
+    fcn(x, fvec, fjac, True)
+    fcn(xp, fvecp, fjac, False)
+    minpack.library.chkder(x, fvec, fjac, xp, fvecp, True, err)
+
+    assert pytest.approx(err) == 15 * [1.0]
+
+    assert driver(fcn, x, fvec, fjac, tol) == 1
+
+    assert pytest.approx(x, abs=100 * tol) == [0.8241058e-1, 0.1133037e1, 0.2343695e1]
+
+
+@pytest.mark.parametrize("driver", [minpack.library.lmder1, minpack.library.lmder])
+def test_lmder_exception(driver):
+    class DummyException(Exception):
+        ...
+
+    def fcn(x, fvec, fjac, jacobian: bool) -> None:
+        raise DummyException()
+
+    x = np.array([1.0, 1.0, 1.0])
+    fvec = np.zeros(15, dtype=np.float64)
+    fjac = np.zeros((3, 15), dtype=np.float64)
+    tol = sqrt(np.finfo(np.float64).eps)
+
+    with pytest.raises(DummyException):
+        driver(fcn, x, fvec, fjac, tol)
+
+
+@pytest.mark.parametrize("driver", [minpack.library.lmdif1, minpack.library.lmdif])
+def test_lmdif(driver):
+    y = np.array(
+        [
+            1.4e-1,
+            1.8e-1,
+            2.2e-1,
+            2.5e-1,
+            2.9e-1,
+            3.2e-1,
+            3.5e-1,
+            3.9e-1,
+            3.7e-1,
+            5.8e-1,
+            7.3e-1,
+            9.6e-1,
+            1.34e0,
+            2.1e0,
+            4.39e0,
+        ]
+    )
+
+    def fcn(x, fvec) -> None:
+        for i in range(fvec.size):
+            tmp1, tmp2 = i + 1, 16 - i - 1
+            tmp3 = tmp2 if i >= 8 else tmp1
+            fvec[i] = y[i] - (x[0] + tmp1 / (x[1] * tmp2 + x[2] * tmp3))
+
+    x = np.array([1.0, 1.0, 1.0])
+    fvec = np.zeros(15, dtype=np.float64)
+    fjac = np.zeros((3, 15), dtype=np.float64)
+    tol = sqrt(np.finfo(np.float64).eps)
+
+    assert driver(fcn, x, fvec, tol) == 1
+
+    assert pytest.approx(x, abs=100 * tol) == [0.8241058e-1, 0.1133037e1, 0.2343695e1]
+
+
+@pytest.mark.parametrize("driver", [minpack.library.lmdif1, minpack.library.lmdif])
+def test_lmdif_exception(driver):
+    class DummyException(Exception):
+        ...
+
+    def fcn(x, fvec) -> None:
+        raise DummyException()
+
+    x = np.array([1.0, 1.0, 1.0])
+    fvec = np.zeros(15, dtype=np.float64)
+    fjac = np.zeros((3, 15), dtype=np.float64)
+    tol = sqrt(np.finfo(np.float64).eps)
+
+    with pytest.raises(DummyException):
+        driver(fcn, x, fvec, tol)
+
+
+@pytest.mark.parametrize("driver", [minpack.library.lmstr1, minpack.library.lmstr])
+def test_lmstr(driver):
+    def fcn(x, fvec, fjrow, row) -> None:
+        if row is None:
+            fvec[0] = 10.0 * (x[1] - x[0] ** 2)
+            fvec[1] = 1.0 - x[0]
+        else:
+            fjrow[0] = -20.0 * x[0] if row == 0 else -1.0
+            fjrow[1] = 10.0 if row == 0 else 0.0
+
+    x = np.array([-1.2, 1.0])
+    fvec = np.zeros(2, dtype=np.float64)
+    fjac = np.zeros((2, 2), dtype=np.float64)
+    tol = sqrt(np.finfo(np.float64).eps)
+
+    assert driver(fcn, x, fvec, fjac, tol) == 4
+
+    assert pytest.approx(x, abs=100 * tol) == 2 * [1.0]
+
+
+@pytest.mark.parametrize("driver", [minpack.library.lmstr1, minpack.library.lmstr])
+def test_lmstr_exception(driver):
+    class DummyException(Exception):
+        ...
+
+    def fcn(x, fvec, fjac, row) -> None:
+        raise DummyException()
+
+    x = np.array([-1.2, 1.0])
+    fvec = np.zeros(2, dtype=np.float64)
+    fjac = np.zeros((2, 2), dtype=np.float64)
+    tol = sqrt(np.finfo(np.float64).eps)
+
+    with pytest.raises(DummyException):
+        driver(fcn, x, fvec, fjac, tol)
diff --git a/python/minpack/typing.py b/python/minpack/typing.py
new file mode 100644
index 0000000..8fc5663
--- /dev/null
+++ b/python/minpack/typing.py
@@ -0,0 +1,13 @@
+# -*- coding: utf-8 -*-
+"""
+Callback signatures for Minpack drivers.
+"""
+
+import numpy as np
+from typing import Optional, Callable
+
+CallableHybrd = Callable[[np.ndarray, np.ndarray], None]
+CallableHybrj = Callable[[np.ndarray, np.ndarray, np.ndarray, bool], None]
+CallableLmder = Callable[[np.ndarray, np.ndarray, np.ndarray, bool], None]
+CallableLmdif = Callable[[np.ndarray, np.ndarray], None]
+CallableLmstr = Callable[[np.ndarray, np.ndarray, np.ndarray, Optional[int]], None]
diff --git a/python/setup.cfg b/python/setup.cfg
new file mode 100644
index 0000000..e47d585
--- /dev/null
+++ b/python/setup.cfg
@@ -0,0 +1,38 @@
+[metadata]
+name = minpack
+version = 2.0.0
+desciption = Python bindings for MINPACK
+long_desciption = file: README.rst
+long_description_content_type = text/x-rst
+url = https://github.com/fortran-lang/minpack
+license = MIT
+license_files =
+    ../LICENSE.txt
+classifiers =
+    Development Status :: 5 - Production
+    Intended Audience :: Science/Research
+    Programming Language :: Fortran
+    Programming Language :: Python :: 3 :: Only
+    Programming Language :: Python :: 3
+    Programming Language :: Python :: 3.6
+    Programming Language :: Python :: 3.7
+    Programming Language :: Python :: 3.8
+    Programming Language :: Python :: 3.9
+    Programming Language :: Python :: 3.10
+
+[options]
+packages = find:
+install_requires =
+    cffi
+    numpy
+tests_require =
+    pytest
+    pytest-cov
+python_requires = >=3.6
+
+[coverage:run]
+omit =
+    */test_*.py
+
+[aliases]
+test = pytest
diff --git a/python/setup.py b/python/setup.py
new file mode 100644
index 0000000..f48c3eb
--- /dev/null
+++ b/python/setup.py
@@ -0,0 +1,7 @@
+from setuptools import setup
+
+setup(
+    cffi_modules=["ffi-builder.py:ffibuilder"],
+    package_data={"minpack": ["_libminpack*.so"]},
+)
+
diff --git a/src/minpack.f90 b/src/minpack.f90
index 562d886..65fe5b1 100644
--- a/src/minpack.f90
+++ b/src/minpack.f90
@@ -1576,7 +1576,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, &
         real(wp), intent(inout) :: Wa1(n) !! work array of length n.
         real(wp), intent(inout) :: Wa2(n) !! work array of length n.
         real(wp), intent(inout) :: Wa3(n) !! work array of length n.
-        real(wp), intent(inout) :: Wa4(m) !! work array of length n.
+        real(wp), intent(inout) :: Wa4(m) !! work array of length m.
 
         integer :: i, iflag, iter, j, l
         real(wp) :: actred, delta, dirder, fnorm, fnorm1, gnorm, par, &
diff --git a/test/api/tester.c b/test/api/tester.c
index 03fb536..d609490 100644
--- a/test/api/tester.c
+++ b/test/api/tester.c
@@ -210,18 +210,18 @@ test_lmder1 (void)
 
     minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 1, err);
     info = 1;
-    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y);
+    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, (void*)y);
     info = 2;
-    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y);
+    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, (void*)y);
     info = 1;
-    trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, y);
+    trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, (void*)y);
     minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 2, err);
 
     for (int i = 0; i < 15; i++) {
         if (!check(err[i], 1.0, tol, "Unexpected derivatives")) return 1;
     }
 
-    minpack_lmder1(trial_lmder_fcn, m, n, x, fvec, fjac, m, tol, &info, ipvt, wa, 30, y);
+    minpack_lmder1(trial_lmder_fcn, m, n, x, fvec, fjac, m, tol, &info, ipvt, wa, 30, (void*)y);
     if (!check(info, 1, "Unexpected info value")) return 1;
     if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1;
     if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1;
@@ -247,11 +247,11 @@ test_lmder (void)
 
     minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 1, err);
     info = 1;
-    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y);
+    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, (void*)y);
     info = 2;
-    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y);
+    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, (void*)y);
     info = 1;
-    trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, y);
+    trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, (void*)y);
     minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 2, err);
 
     for (int i = 0; i < 15; i++) {
@@ -259,7 +259,7 @@ test_lmder (void)
     }
 
     minpack_lmder(trial_lmder_fcn, m, n, x, fvec, fjac, m, tol, tol, 0.0, 2000, diag, 1,
-            100.0, 0, &info, &nfev, &njev, ipvt, qtf, wa1, wa2, wa3, wa4, y);
+            100.0, 0, &info, &nfev, &njev, ipvt, qtf, wa1, wa2, wa3, wa4, (void*)y);
     if (!check(info, 1, "Unexpected info value")) return 1;
     if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1;
     if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1;
@@ -298,7 +298,7 @@ test_lmdif1 (void)
     int lwa = m*n + 5*n + m;
     double wa[lwa];
 
-    minpack_lmdif1(trial_lmdif_fcn, 15, 3, x, fvec, tol, &info, ipvt, wa, lwa, y);
+    minpack_lmdif1(trial_lmdif_fcn, 15, 3, x, fvec, tol, &info, ipvt, wa, lwa, (void*)y);
     if (!check(info, 1, "Unexpected info value")) return 1;
     if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1;
     if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1;
@@ -321,7 +321,7 @@ test_lmdif (void)
     double fjac[m*n], diag[n], qtf[n], wa1[n], wa2[n], wa3[n], wa4[m];
 
     minpack_lmdif(trial_lmdif_fcn, 15, 3, x, fvec, tol, tol, 0.0, 2000, 0.0, diag, 1,
-            100.0, 0, &info, &nfev, fjac, 15, ipvt, qtf, wa1, wa2, wa3, wa4, y);
+            100.0, 0, &info, &nfev, fjac, 15, ipvt, qtf, wa1, wa2, wa3, wa4, (void*)y);
     if (!check(info, 1, "Unexpected info value")) return 1;
     if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1;
     if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1;
diff --git a/test/meson.build b/test/meson.build
index 6279d4b..1f744e6 100644
--- a/test/meson.build
+++ b/test/meson.build
@@ -19,7 +19,7 @@ foreach t : tests
   )
 endforeach
 
-if add_languages('c', required: false, native: false)
+if has_cc
   test(
     'c-api',
     executable(