NumPy core math library#

The numpy core math library (npymath) is a first step in this direction. Thislibrary contains most math-related C99 functionality, which can be used onplatforms where C99 is not well supported. The core math functions have thesame API as the C99 ones, except for thenpy_* prefix.

The available functions are defined in<numpy/npy_math.h> - please refer tothis header when in doubt.

Note

An effort is underway to makenpymath smaller (since C99 compatibilityof compilers has improved over time) and more easily vendorable or usable asa header-only dependency. That will avoid problems with shipping a staticlibrary built with a compiler which may not match the compiler used by adownstream package or end user. Seegh-20880 for details.

Floating point classification#

NPY_NAN#

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

NPY_INFINITY#

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

NPY_PZERO#

This macro is defined to positive zero. The corresponding single andextension 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). Thecorresponding single and extension precision macro are available with thesuffix F and L.

npy_isnan(x)#

This is an alias for C99 isnan: works for single, doubleand 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 aNaN nor an infinity.

npy_isinf(x)#

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

npy_signbit(x)#

This is an alias for C99 signbit: works for single, doubleand 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 signas y. Works for any value, including inf and nan. Single and extendedprecisions are available with suffix f and l.

Useful math constants#

The following math constants are available innpy_math.h. Singleand extended precision are also available by adding thef andl 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})\)

Low-level floating point manipulation#

Those can be useful for precise floating point comparison.

doublenpy_nextafter(doublex,doubley)#

This is an alias to C99 nextafter: return next representablefloating point value from x in the direction of y. Single and extendedprecisions are available with suffix f and l.

doublenpy_spacing(doublex)#

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

voidnpy_set_floatstatus_divbyzero()#

Set the divide by zero floating point exception

voidnpy_set_floatstatus_overflow()#

Set the overflow floating point exception

voidnpy_set_floatstatus_underflow()#

Set the underflow floating point exception

voidnpy_set_floatstatus_invalid()#

Set the invalid floating point exception

intnpy_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 thatnpy_get_floatstatus_barrier is preferable as it preventsaggressive compiler optimizations reordering the call relative tothe code setting the status, which could lead to incorrect results.

intnpy_get_floatstatus_barrier(char*)#

Get floating point status. A pointer to a local variable is passed in toprevent aggressive compiler optimizations from reordering this function callrelative to the code setting the status, which could lead to incorrectresults.

Returns a bitmask with following possible flags:

  • NPY_FPE_DIVIDEBYZERO

  • NPY_FPE_OVERFLOW

  • NPY_FPE_UNDERFLOW

  • NPY_FPE_INVALID

intnpy_clear_floatstatus()#

Clears the floating point status. Returns the previous status mask.

Note thatnpy_clear_floatstatus_barrier is preferable as itprevents aggressive compiler optimizations reordering the call relative tothe code setting the status, which could lead to incorrect results.

intnpy_clear_floatstatus_barrier(char*)#

Clears the floating point status. A pointer to a local variable is passed in toprevent aggressive compiler optimizations from reordering this function call.Returns the previous status mask.

Support for complex numbers#

C99-like complex functions have been added. Those can be used if you wish toimplement portable C extensions. Since NumPy 2.0 we use C99 complex types asthe underlying type:

typedefdouble_Complexnpy_cdouble;typedeffloat_Complexnpy_cfloat;typedeflongdouble_Complexnpy_clongdouble;

MSVC does not support the_Complex type itself, but has added support forthe C99complex.h header by providing its own implementation. Thus, underMSVC, the equivalent MSVC types will be used:

typedef_Dcomplexnpy_cdouble;typedef_Fcomplexnpy_cfloat;typedef_Lcomplexnpy_clongdouble;

Because MSVC still does not support C99 syntax for initializing a complexnumber, you need to restrict to C90-compatible syntax, e.g.:

/* a = 1 + 2i \*/npy_complexa=npy_cpack(1,2);npy_complexb;b=npy_log(a);

A few utilities have also been added innumpy/npy_math.h, in order to retrieve or set the real or the imaginarypart of a complex number:

npy_cdoublec;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 touse C99 complex types. Up until now the following was being used to representcomplex types:

typedefstruct{doublereal,imag;}npy_cdouble;typedefstruct{floatreal,imag;}npy_cfloat;typedefstruct{npy_longdoublereal,imag;}npy_clongdouble;

Using thestruct representation ensured that complex numbers could be usedon all platforms, even the ones without support for built-in complex types. Italso meant that a static library had to be shipped together with NumPy toprovide a C99 compatibility layer for downstream packages to use. In recentyears however, support for native complex types has been improved immensely,with MSVC adding built-in support for thecomplex.h header in 2019.

To ease cross-version compatibility, macros that use the new set APIs havebeen 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 innumpy/npy_2_complexcompat.h. Itchecks whether the macros exist, and falls back to the 1.x syntax in case theydon’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 copy-pastethe compatibility layer code into their own sources and use that, so thatthey can continue to support both NumPy 1.x and 2.x without issues. Note alsothat thecomplex.h header is included innumpy/npy_common.h, whichmakescomplex 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 ownPython extension, you need to add thenpymath compile and link options to yourextension. The exact steps to take will depend on the build system you are using.The generic steps to take are:

  1. Add the numpy include directory (= the value ofnp.get_include()) toyour include directories,

  2. Thenpymath static library resides in thelib directory right nextto numpy’s include directory (i.e.,pathlib.Path(np.get_include())/'..'/'lib'). Add that to your library search directories,

  3. Link withlibnpymath andlibm.

Note

Keep in mind that when you are cross compiling, you must use thenumpyfor the platform you are building for, not the native one for the buildmachine. Otherwise you pick up a static library built for the wrongarchitecture.

When you build withnumpy.distutils (deprecated), then use this in yoursetup.py:

>>>fromnumpy.distutils.misc_utilimportget_info>>>info=get_info('npymath')>>>_=config.add_extension('foo',sources=['foo.c'],extra_info=info)

In other words, the usage ofinfo is exactly the same as when usingblas_info and co.

When you are building withMeson, 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],

Half-precision functions#

The header file<numpy/halffloat.h> provides functions to work withIEEE 754-2008 16-bit floating point values. While this format isnot typically used for numerical computations, it is useful forstoring values which require floating point but do not need much precision.It can also be used as an educational tool to understand the natureof floating point round-off error.

Like for other types, NumPy includes a typedef npy_half for the 16 bitfloat. Unlike for most of the other types, you cannot use this as anormal 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 comparisonbetween 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 valuesaccessible by including<numpy/halffloat.h> and linking tonpymath.For functions that are not provided directly, such as the arithmeticoperations, the preferred method is to convert to floator double and back again, as in the following example.

npy_halfsum(intn,npy_half*array){floatret=0;while(n--){ret+=npy_half_to_float(*array++);}returnnpy_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.

floatnpy_half_to_float(npy_halfh)#

Converts a half-precision float to a single-precision float.

doublenpy_half_to_double(npy_halfh)#

Converts a half-precision float to a double-precision float.

npy_halfnpy_float_to_half(floatf)#

Converts a single-precision float to a half-precision float. Thevalue is rounded to the nearest representable half, with ties goingto the nearest even. If the value is too small or too big, thesystem’s floating point underflow or overflow bit will be set.

npy_halfnpy_double_to_half(doubled)#

Converts a double-precision float to a half-precision float. Thevalue is rounded to the nearest representable half, with ties goingto the nearest even. If the value is too small or too big, thesystem’s floating point underflow or overflow bit will be set.

intnpy_half_eq(npy_halfh1,npy_halfh2)#

Compares two half-precision floats (h1 == h2).

intnpy_half_ne(npy_halfh1,npy_halfh2)#

Compares two half-precision floats (h1 != h2).

intnpy_half_le(npy_halfh1,npy_halfh2)#

Compares two half-precision floats (h1 <= h2).

intnpy_half_lt(npy_halfh1,npy_halfh2)#

Compares two half-precision floats (h1 < h2).

intnpy_half_ge(npy_halfh1,npy_halfh2)#

Compares two half-precision floats (h1 >= h2).

intnpy_half_gt(npy_halfh1,npy_halfh2)#

Compares two half-precision floats (h1 > h2).

intnpy_half_eq_nonan(npy_halfh1,npy_halfh2)#

Compares two half-precision floats that are known to not be NaN (h1 == h2). Ifa value is NaN, the result is undefined.

intnpy_half_lt_nonan(npy_halfh1,npy_halfh2)#

Compares two half-precision floats that are known to not be NaN (h1 < h2). Ifa value is NaN, the result is undefined.

intnpy_half_le_nonan(npy_halfh1,npy_halfh2)#

Compares two half-precision floats that are known to not be NaN (h1 <= h2). Ifa value is NaN, the result is undefined.

intnpy_half_iszero(npy_halfh)#

Tests whether the half-precision float has a value equal to zero. This may be slightlyfaster than calling npy_half_eq(h, NPY_ZERO).

intnpy_half_isnan(npy_halfh)#

Tests whether the half-precision float is a NaN.

intnpy_half_isinf(npy_halfh)#

Tests whether the half-precision float is plus or minus Inf.

intnpy_half_isfinite(npy_halfh)#

Tests whether the half-precision float is finite (not NaN or Inf).

intnpy_half_signbit(npy_halfh)#

Returns 1 is h is negative, 0 otherwise.

npy_halfnpy_half_copysign(npy_halfx,npy_halfy)#

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

npy_halfnpy_half_spacing(npy_halfh)#

This is the same for half-precision float as npy_spacing and npy_spacingfdescribed in the low-level floating point section.

npy_halfnpy_half_nextafter(npy_halfx,npy_halfy)#

This is the same for half-precision float as npy_nextafter and npy_nextafterfdescribed in the low-level floating point section.

npy_uint16npy_floatbits_to_halfbits(npy_uint32f)#

Low-level function which converts a 32-bit single-precision float, storedas a uint32, into a 16-bit half-precision float.

npy_uint16npy_doublebits_to_halfbits(npy_uint64d)#

Low-level function which converts a 64-bit double-precision float, storedas a uint64, into a 16-bit half-precision float.

npy_uint32npy_halfbits_to_floatbits(npy_uint16h)#

Low-level function which converts a 16-bit half-precision floatinto a 32-bit single-precision float, stored as a uint32.

npy_uint64npy_halfbits_to_doublebits(npy_uint16h)#

Low-level function which converts a 16-bit half-precision floatinto a 64-bit double-precision float, stored as a uint64.