NumPy core math library#
The numpy core math library (npymath
) is a first step in this direction. This
library contains most mathrelated C99 functionality, which can be used on
platforms where C99 is not well supported. The core math functions have the
same API as the C99 ones, except for the npy_*
prefix.
The available functions are defined in <numpy/npy_math.h>
 please refer to
this header when in doubt.
Note
An effort is underway to make npymath
smaller (since C99 compatibility
of compilers has improved over time) and more easily vendorable or usable as
a headeronly dependency. That will avoid problems with shipping a static
library built with a compiler which may not match the compiler used by a
downstream package or end user. See
gh20880 for details.
Floating point classification#

NPY_NAN#
This macro is defined to a NaN (Not a Number), and is guaranteed to have the signbit unset (‘positive’ NaN). The corresponding single and extension precision macro are available with the suffix F and L.

NPY_INFINITY#
This macro is defined to a positive inf. The corresponding single and extension precision macro are available with the suffix F and L.

NPY_PZERO#
This macro is defined to positive zero. The corresponding single and extension precision macro are available with the suffix F and L.

NPY_NZERO#
This macro is defined to negative zero (that is with the sign bit set). The corresponding single and extension precision macro are available with the suffix F and L.

npy_isnan(x)#
This is an alias for C99 isnan: works for single, double and extended precision, and return a non 0 value if x is a NaN.

npy_isfinite(x)#
This is an alias for C99 isfinite: works for single, double and extended precision, and return a non 0 value if x is neither a NaN nor an infinity.

npy_isinf(x)#
This is an alias for C99 isinf: works for single, double and extended precision, and return a non 0 value if x is infinite (positive and negative).

npy_signbit(x)#
This is an alias for C99 signbit: works for single, double and extended precision, and return a non 0 value if x has the signbit set (that is the number is negative).

npy_copysign(x, y)#
This is an alias for C99 copysign: return x with the same sign as y. Works for any value, including inf and nan. Single and extended precisions are available with suffix f and l.
Useful math constants#
The following math constants are available in npy_math.h
. Single
and extended precision are also available by adding the f
and
l
suffixes respectively.

NPY_E#
Base of natural logarithm (\(e\))

NPY_LOG2E#
Logarithm to base 2 of the Euler constant (\(\frac{\ln(e)}{\ln(2)}\))

NPY_LOG10E#
Logarithm to base 10 of the Euler constant (\(\frac{\ln(e)}{\ln(10)}\))

NPY_LOGE2#
Natural logarithm of 2 (\(\ln(2)\))

NPY_LOGE10#
Natural logarithm of 10 (\(\ln(10)\))

NPY_PI#
Pi (\(\pi\))

NPY_PI_2#
Pi divided by 2 (\(\frac{\pi}{2}\))

NPY_PI_4#
Pi divided by 4 (\(\frac{\pi}{4}\))

NPY_1_PI#
Reciprocal of pi (\(\frac{1}{\pi}\))

NPY_2_PI#
Two times the reciprocal of pi (\(\frac{2}{\pi}\))

NPY_EULER#
 The Euler constant
\(\lim_{n\rightarrow\infty}({\sum_{k=1}^n{\frac{1}{k}}\ln n})\)
Lowlevel floating point manipulation#
Those can be useful for precise floating point comparison.

double npy_nextafter(double x, double y)#
This is an alias to C99 nextafter: return next representable floating point value from x in the direction of y. Single and extended precisions are available with suffix f and l.

double npy_spacing(double x)#
This is a function equivalent to Fortran intrinsic. Return distance between x and next representable floating point value from x, e.g. spacing(1) == eps. spacing of nan and +/ inf return nan. Single and extended precisions are available with suffix f and l.

void npy_set_floatstatus_divbyzero()#
Set the divide by zero floating point exception

void npy_set_floatstatus_overflow()#
Set the overflow floating point exception

void npy_set_floatstatus_underflow()#
Set the underflow floating point exception

void npy_set_floatstatus_invalid()#
Set the invalid floating point exception

int npy_get_floatstatus()#
Get floating point status. Returns a bitmask with following possible flags:
NPY_FPE_DIVIDEBYZERO
NPY_FPE_OVERFLOW
NPY_FPE_UNDERFLOW
NPY_FPE_INVALID
Note that
npy_get_floatstatus_barrier
is preferable as it prevents aggressive compiler optimizations reordering the call relative to the code setting the status, which could lead to incorrect results.

int npy_get_floatstatus_barrier(char*)#
Get floating point status. A pointer to a local variable is passed in to prevent aggressive compiler optimizations from reordering this function call relative to the code setting the status, which could lead to incorrect results.
Returns a bitmask with following possible flags:
NPY_FPE_DIVIDEBYZERO
NPY_FPE_OVERFLOW
NPY_FPE_UNDERFLOW
NPY_FPE_INVALID

int npy_clear_floatstatus()#
Clears the floating point status. Returns the previous status mask.
Note that
npy_clear_floatstatus_barrier
is preferable as it prevents aggressive compiler optimizations reordering the call relative to the code setting the status, which could lead to incorrect results.

int npy_clear_floatstatus_barrier(char*)#
Clears the floating point status. A pointer to a local variable is passed in to prevent aggressive compiler optimizations from reordering this function call. Returns the previous status mask.
Support for complex numbers#
C99like complex functions have been added. Those can be used if you wish to implement portable C extensions. Since NumPy 2.0 we use C99 complex types as the underlying type:
typedef double _Complex npy_cdouble;
typedef float _Complex npy_cfloat;
typedef long double _Complex npy_clongdouble;
MSVC does not support the _Complex
type itself, but has added support for
the C99 complex.h
header by providing its own implementation. Thus, under
MSVC, the equivalent MSVC types will be used:
typedef _Dcomplex npy_cdouble;
typedef _Fcomplex npy_cfloat;
typedef _Lcomplex npy_clongdouble;
Because MSVC still does not support C99 syntax for initializing a complex number, you need to restrict to C90compatible syntax, e.g.:
/* a = 1 + 2i \*/
npy_complex a = npy_cpack(1, 2);
npy_complex b;
b = npy_log(a);
A few utilities have also been added in
numpy/npy_math.h
, in order to retrieve or set the real or the imaginary
part of a complex number:
npy_cdouble c;
npy_csetreal(&c, 1.0);
npy_csetimag(&c, 0.0);
printf("%d + %di\n", npy_creal(c), npy_cimag(c));
Changed in version 2.0.0: The underlying C types for all of numpy’s complex types have been changed to use C99 complex types. Up until now the following was being used to represent complex types:
typedef struct { double real, imag; } npy_cdouble;
typedef struct { float real, imag; } npy_cfloat;
typedef struct {npy_longdouble real, imag;} npy_clongdouble;
Using the struct
representation ensured that complex numbers could be used
on all platforms, even the ones without support for builtin complex types. It
also meant that a static library had to be shipped together with NumPy to
provide a C99 compatibility layer for downstream packages to use. In recent
years however, support for native complex types has been improved immensely,
with MSVC adding builtin support for the complex.h
header in 2019.
To ease crossversion compatibility, macros that use the new set APIs have been added.
#define NPY_CSETREAL(z, r) npy_csetreal(z, r)
#define NPY_CSETIMAG(z, i) npy_csetimag(z, i)
A compatibility layer is also provided in numpy/npy_2_complexcompat.h
. It
checks whether the macros exist, and falls back to the 1.x syntax in case they
don’t.
#include <numpy/npy_math.h>
#ifndef NPY_CSETREALF
#define NPY_CSETREALF(c, r) (c)>real = (r)
#endif
#ifndef NPY_CSETIMAGF
#define NPY_CSETIMAGF(c, i) (c)>imag = (i)
#endif
We suggest all downstream packages that need this functionality to copypaste
the compatibility layer code into their own sources and use that, so that
they can continue to support both NumPy 1.x and 2.x without issues. Note also
that the complex.h
header is included in numpy/npy_common.h
, which
makes complex
a reserved keyword.
Linking against the core math library in an extension#
To use the core math library that NumPy ships as a static library in your own
Python extension, you need to add the npymath
compile and link options to your
extension. The exact steps to take will depend on the build system you are using.
The generic steps to take are:
Add the numpy include directory (= the value of
np.get_include()
) to your include directories,The
npymath
static library resides in thelib
directory right next to numpy’s include directory (i.e.,pathlib.Path(np.get_include()) / '..' / 'lib'
). Add that to your library search directories,Link with
libnpymath
andlibm
.
Note
Keep in mind that when you are cross compiling, you must use the numpy
for the platform you are building for, not the native one for the build
machine. Otherwise you pick up a static library built for the wrong
architecture.
When you build with numpy.distutils
(deprecated), then use this in your setup.py
:
>>> from numpy.distutils.misc_util import get_info >>> info = get_info('npymath') >>> _ = config.add_extension('foo', sources=['foo.c'], extra_info=info)
In other words, the usage of info
is exactly the same as when using
blas_info
and co.
When you are building with Meson, use:
# Note that this will get easier in the future, when Meson has
# support for numpy built in; most of this can then be replaced
# by `dependency('numpy')`.
incdir_numpy = run_command(py3,
[
'c',
'import os; os.chdir(".."); import numpy; print(numpy.get_include())'
],
check: true
).stdout().strip()
inc_np = include_directories(incdir_numpy)
cc = meson.get_compiler('c')
npymath_path = incdir_numpy / '..' / 'lib'
npymath_lib = cc.find_library('npymath', dirs: npymath_path)
py3.extension_module('module_name',
...
include_directories: inc_np,
dependencies: [npymath_lib],
Halfprecision functions#
The header file <numpy/halffloat.h>
provides functions to work with
IEEE 7542008 16bit floating point values. While this format is
not typically used for numerical computations, it is useful for
storing values which require floating point but do not need much precision.
It can also be used as an educational tool to understand the nature
of floating point roundoff error.
Like for other types, NumPy includes a typedef npy_half for the 16 bit float. Unlike for most of the other types, you cannot use this as a normal type in C, since it is a typedef for npy_uint16. For example, 1.0 looks like 0x3c00 to C, and if you do an equality comparison between the different signed zeros, you will get 0.0 != 0.0 (0x8000 != 0x0000), which is incorrect.
For these reasons, NumPy provides an API to work with npy_half values
accessible by including <numpy/halffloat.h>
and linking to npymath
.
For functions that are not provided directly, such as the arithmetic
operations, the preferred method is to convert to float
or double and back again, as in the following example.
npy_half sum(int n, npy_half *array) {
float ret = 0;
while(n) {
ret += npy_half_to_float(*array++);
}
return npy_float_to_half(ret);
}
External Links:

NPY_HALF_ZERO#
This macro is defined to positive zero.

NPY_HALF_PZERO#
This macro is defined to positive zero.

NPY_HALF_NZERO#
This macro is defined to negative zero.

NPY_HALF_ONE#
This macro is defined to 1.0.

NPY_HALF_NEGONE#
This macro is defined to 1.0.

NPY_HALF_PINF#
This macro is defined to +inf.

NPY_HALF_NINF#
This macro is defined to inf.

NPY_HALF_NAN#
This macro is defined to a NaN value, guaranteed to have its sign bit unset.

npy_half npy_float_to_half(float f)#
Converts a singleprecision float to a halfprecision float. The value is rounded to the nearest representable half, with ties going to the nearest even. If the value is too small or too big, the system’s floating point underflow or overflow bit will be set.

npy_half npy_double_to_half(double d)#
Converts a doubleprecision float to a halfprecision float. The value is rounded to the nearest representable half, with ties going to the nearest even. If the value is too small or too big, the system’s floating point underflow or overflow bit will be set.

int npy_half_eq_nonan(npy_half h1, npy_half h2)#
Compares two halfprecision floats that are known to not be NaN (h1 == h2). If a value is NaN, the result is undefined.

int npy_half_lt_nonan(npy_half h1, npy_half h2)#
Compares two halfprecision floats that are known to not be NaN (h1 < h2). If a value is NaN, the result is undefined.

int npy_half_le_nonan(npy_half h1, npy_half h2)#
Compares two halfprecision floats that are known to not be NaN (h1 <= h2). If a value is NaN, the result is undefined.

int npy_half_iszero(npy_half h)#
Tests whether the halfprecision float has a value equal to zero. This may be slightly faster than calling npy_half_eq(h, NPY_ZERO).

int npy_half_isfinite(npy_half h)#
Tests whether the halfprecision float is finite (not NaN or Inf).

npy_half npy_half_copysign(npy_half x, npy_half y)#
Returns the value of x with the sign bit copied from y. Works for any value, including Inf and NaN.

npy_half npy_half_spacing(npy_half h)#
This is the same for halfprecision float as npy_spacing and npy_spacingf described in the lowlevel floating point section.

npy_half npy_half_nextafter(npy_half x, npy_half y)#
This is the same for halfprecision float as npy_nextafter and npy_nextafterf described in the lowlevel floating point section.

npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)#
Lowlevel function which converts a 32bit singleprecision float, stored as a uint32, into a 16bit halfprecision float.

npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)#
Lowlevel function which converts a 64bit doubleprecision float, stored as a uint64, into a 16bit halfprecision float.

npy_uint32 npy_halfbits_to_floatbits(npy_uint16 h)#
Lowlevel function which converts a 16bit halfprecision float into a 32bit singleprecision float, stored as a uint32.

npy_uint64 npy_halfbits_to_doublebits(npy_uint16 h)#
Lowlevel function which converts a 16bit halfprecision float into a 64bit doubleprecision float, stored as a uint64.