Integrating Fortran and Python with F2PY

Introduction

With future projects demanding mathematical special functions, I’ve started implementing them in Fortran for its computational efficiency. However, not everything will be built with pure Fortran, as I also plan to integrate it with the versatility of Python through F2PY. So the purpose of this post is to demonstrate how to bridge Fortran and Python with F2PY.

Methods

Tools

F2PY is a tool included with the NumPy package. To function, it requires a Fortran compiler, such as gfortran. Furthermore, for Python versions 3.12 and above, F2PY demands the installation of meson backend due to changes in NumPy’s build process.

Fortran code

The Fortran subroutine to be wrapped by F2PY contains the numerical implementation of the exponential integral:

$$ \text{Ei}(x) = \int_{-\infty}^x \frac{e^t}{t} \, dt, \quad x \in \mathbb{R}, x > 0 $$

The expressions for numerical evaluation of $\text{Ei}$ and many other special functions can be found on the book Computation of Special Functions by Zhang and Jin. The numerical implementation given below is a refit of the code found in the book.

expint.f90

subroutine expi(x, ei)
  ! Exponential integral Ei(x).
  !
  ! Parameters
  ! ----------
  ! x : real(real64)
  !   Real number ≥ 0.
  !    
  ! Returns
  ! -------
  ! ei : real(real64) 
  !   Exponential integral of x.
  !
  ! Author
  ! ------
  ! Rodrigo Castro (GitHub: rodpcastro)
  !
  ! History
  ! -------
  ! 24-04-2025 - Rodrigo Castro - Original code
  !
  ! Reference
  ! ----------
  ! Shanjie Zhang, Jianming Jin (1996). Computation of Special Functions.

  use, intrinsic :: iso_fortran_env, only: int16, real64

  implicit none
  real(real64), parameter :: gm = 0.5772156649015329d0 ! Euler's constant
  real(real64), intent(in) :: x
  real(real64), intent(out) :: ei
  real(real64) :: r
  integer(int16) :: n

  if (x == 0.0d0) then
    ei = -1.0d+300
  else if (x <= 40.0d0) then
    ei = x
    r = x
    do n = 2, 100
      r = r * x * (n-1) / n**2
      ei = ei + r
      if (abs(r) <= 1.0d-15*abs(ei)) exit
    end do
    ei = ei + gm + log(x)
  else
    ei = 1.0d0
    r = 1.0d0
    do n = 1, 20
      r = r * n / x
      ei = ei + r
    end do
    ei = exp(x) * ei / x 
  end if

end subroutine expi

The expi subroutine might not display the most correct way of writing subroutines in Fortran, but it reflects the adaptations suggested by F2PY examples for correct wrapping:

Avoid Fortran modules: In a typical Fortran setup, we might place the expi subroutine inside a module for organization. However, since F2PY automatically wraps the subroutine into a Python module, this would create a redundant nested structure. To keep the Python interface clean, the subroutine is defined outside of a Fortran module.
Handling iso_fortran_env parameters: int16 and real64 had to be used inside the subroutine. Unclear why, these parameters were the cause of parsing errors in F2PY when used outside of the subroutine (within a module).

F2PY

To ensure that F2PY correctly interprets Fortran int16 and real64 kinds, they need to be mapped to their equivalent C types. This is done by using a .f2py_c2map file placed in the same directory as the expint.f90 file:

# .f2py_c2map
dict(real=dict(real64='double'), integer=dict(int16='int'))

After that, all that remains is to build our extension module by running the following command:

f2py -c expint.f90 -m expint

shell output

The Meson build system
Version: 1.7.2
Source dir: C:\Users\Rodrigo\AppData\Local\Temp\tmpc5z677mk
Build dir: C:\Users\Rodrigo\AppData\Local\Temp\tmpc5z677mk\bbdir
Build type: native build
Project name: expint
Project version: 0.1
Fortran compiler for the host machine: gfortran (gcc 14.2.0 "GNU Fortran (Rev2, Built by MSYS2 project) 14.2.0")
Fortran linker for the host machine: gfortran ld.bfd 2.43.1
C compiler for the host machine: cc (gcc 14.2.0 "cc (Rev2, Built by MSYS2 project) 14.2.0")
C linker for the host machine: cc ld.bfd 2.43.1
Host machine cpu family: x86_64
Host machine cpu: x86_64
Program C:\Users\Rodrigo\miniconda3\python.exe found: YES (C:\Users\Rodrigo\miniconda3\python.exe)
Run-time dependency python found: YES 3.12
Library quadmath found: YES
Build targets in project: 1

Found ninja-1.11.1.git.kitware.jobserver-1 at C:\Users\Rodrigo\miniconda3\Scripts\ninja.EXE
ninja: Entering directory `C:/Users/Rodrigo/AppData/Local/Temp/tmpc5z677mk/bbdir'
[1/6] Compiling C object expint.cp312-win_amd64.pyd.p/expintmodule.c.obj
[2/6] Scanning modules
[3/6] Compiling Fortran object expint.cp312-win_amd64.pyd.p/expint-f2pywrappers.f.obj
[4/6] Compiling Fortran object expint.cp312-win_amd64.pyd.p/expint.f90.obj
[5/6] Compiling C object expint.cp312-win_amd64.pyd.p/fdef2c95a4dd72a4b86c832312e3647fd2047c6d_.._.._f2py_src_fortranobject.c.obj
[6/6] Linking target expint.cp312-win_amd64.pyd
INFO: autodetecting backend as ninja
INFO: calculating backend command to run: C:\Users\Rodrigo\miniconda3\Scripts\ninja.EXE -C C:/Users/Rodrigo/AppData/Local/Temp/tmpc5z677mk/bbdir
Cannot use distutils backend with Python>=3.12, using meson backend instead.
Using meson backend
Will pass --lower to f2py
See https://numpy.org/doc/stable/f2py/buildtools/meson.html
Reading f2cmap from '.f2py_f2cmap' ...
	Mapping "real(kind=real64)" to "double"
	Mapping "integer(kind=int16)" to "int"
Successfully applied user defined f2cmap changes
Reading fortran codes...
	Reading file 'expint.f90' (format:free)
{'before': '', 'this': 'use', 'after': ', intrinsic :: iso_fortran_env, only: int16, real64 '}
Line #27 in expint.f90:"  use, intrinsic :: iso_fortran_env, only: int16, real64 "
	analyzeline: Could not crack the use statement.
Post-processing...
	Block: expint
			Block: expi
In: :expint:expint.f90:expi
param_eval: got "eval() arg 1 must be a string, bytes or code object" on 8
In: :expint:expint.f90:expi
param_eval: got "eval() arg 1 must be a string, bytes or code object" on 8
Applying post-processing hooks...
  character_backward_compatibility_hook
Post-processing (stage 2)...
Building modules...
    Building module "expint"...
    Generating possibly empty wrappers"
    Maybe empty "expint-f2pywrappers.f"
        Constructing wrapper function "expi"...
          ei = expi(x)
    Wrote C/API module "expint" to file ".\expintmodule.c"

The command produces a shared library (expint.cp312-win_amd64.pyd). The expi procedure can now be called from a Python module named expint.

Results

To demonstrate that the expi procedure works as expected, its outputs are compared against mpmath, a Python library for arbitrary precision arithmetic. The mpmath.ei function is used with 16 digits of precision and the relative error between expi and mpmath.ei is computed and displayed through the following Python code:

expint_test.py

from expint import expi
import mpmath
import matplotlib.pyplot as plt
import numpy as np

mpmath.dps = 16 # mpmath precision

xv = np.logspace(-10, 2, 1000, dtype=np.float64)
ei_expint = np.empty(xv.shape, dtype=np.float64)
ei_mpmath = np.empty(xv.shape, dtype=np.float64)

for i, x in enumerate(xv):
    ei_expint[i] = expi(x)
    ei_mpmath[i] = mpmath.ei(x)

ei_diff = np.abs(ei_expint - ei_mpmath)/np.abs(ei_mpmath)

fig, ax = plt.subplots(figsize=(7,4))
ax.loglog(xv, ei_diff, 'k')
ax.set_xticks(np.logspace(-10, 2, 5))
ax.set_yticks(np.logspace(-16, -13, 4))
ax.set_xlabel('x')
ax.set_ylabel(r'$ \frac{|Ei_{expint} - Ei_{mpmath}|}{|Ei_{mpmath}|} $')
ax.set_title('Relative error')
plt.savefig('expi_error.svg')

The resulting plot shows that the relative difference between expi and mpmath.ei is very small (≤ $10^{-13}$) for $10^{-10} ≤ x ≤ 10^{2}$.

References

NumPy Developers. 2024. F2PY. Distributed as part of NumPy. https://numpy.org/doc/stable/f2py/index.html. (2025).

Introduction#

Methods#

Tools#

Fortran code#

F2PY#

Results#

References#

Appendices#