tensor_ops.hpp File Reference

Element-wise arithmetic with NumPy's broadcasting, reductions, reshaping and creation of tensors. More...

#include <algorithm>
#include <array>
#include <bit>
#include <cmath>
#include <complex>
#include <concepts>
#include <cstddef>
#include <cstring>
#include <functional>
#include <limits>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
#include "simd.hpp"
#include "tensor.hpp"

Include dependency graph for tensor_ops.hpp:

Go to the source code of this file.

Classes
struct	numerical::detail::axis_walk
	The extents and strides that walk a tensor of shape with axis axis taken out: outer positions before it, inner after it, the axis itself extent long at stride inner. More...

Namespaces
namespace	numerical
	The NumericAL public API.
namespace	numerical::detail

Functions
std::string	numerical::detail::shape_text (const std::vector< std::size_t > &shape)
	A shape as NumPy prints it, (2, 3), for a message.
std::vector< std::size_t >	numerical::detail::broadcast_strides (const std::vector< std::size_t > &shape, std::size_t rank)
	The row-major strides of shape, with a stride of 0 on each axis of extent 1 (it broadcasts), aligned on the last axis of a result of rank rank (missing leading axes stride 0 too).
axis_walk	numerical::detail::walk_along (const std::vector< std::size_t > &shape, std::size_t axis)
	axis_walk of shape along axis; outer and inner 0 where the result is empty.
template<std::floating_point R>
R	numerical::detail::numpy_muladd (R x, R y, R z, bool fused)
	x · y + z: one fused multiply-add where fused, else the product rounded first (statements apart, so no compiler contracts them).
template<split_complex C>
C	numerical::detail::numpy_multiply (const C &a, const C &b, bool fused)
	a · b of complex numbers as NumPy's array loop computes it: ai·bi and ai·br rounded, then ar·br − ai·bi and ar·bi + ai·br, each by a fused multiply-add where fused (NumPy's vectors on arm64 and on x86-64 with FMA), else each product rounded first. No recovery of an infinity from a NaN, as std::complex's (C Annex G) makes: 1j · (inf + inf·j) is NaN, (1e308 + 1e308j) · (1e308 − 1e308j) inf + inf·j when fused.
template<split_complex C>
C	numerical::detail::numpy_divide (const C &a, const C &b, bool fused)
	a / b of complex numbers as NumPy's array loop computes it: Smith's division by the larger part of b (an exact zero divided part by part, to an infinity or a NaN), each x + y·z of it one fused multiply-add where fused (NumPy's build contracts them on arm64, not on x86-64).
bool	numerical::detail::numpy_fuses_products () noexcept
	Whether NumPy's complex product fuses its multiply-adds on this processor: on arm64, and on x86-64 where its AVX2 loop runs (FMA with it); elsewhere its loop rounds each product.
bool	numerical::detail::numpy_fuses_scalar_products () noexcept
	Whether NumPy's complex product fuses its multiply-adds in its scalar loop (where numpy_loop_is_scalar holds), as its PyPI wheels are built: that loop sits beside the vector one, and clang (macOS) and MSVC (Windows) contract it where the vector loop has FMA; GCC (Linux, ISO C) does not.
bool	numerical::detail::numpy_fuses_quotients () noexcept
	Whether NumPy's complex quotient fuses its multiply-adds, as its PyPI wheels are built: on macOS arm64, whose clang contracts them (FMA is arm64's baseline); not on Linux, whose GCC does not, nor on x86-64, whose baseline has no FMA.
template<typename T >
bool	numerical::detail::numpy_loop_is_scalar (const tensor< T > &a, const tensor< T > &b) noexcept
	Whether NumPy multiplies a by b in its scalar loop rather than its vector one: when the result is one element and an operand of rank above 0 is broadcast to a higher rank, NumPy's iterator hands its loop an output step of 0, which the vector loop refuses.
template<std::floating_point T>
T	numerical::detail::divided_by_count (T value, std::size_t n)
	A mean: value divided by the count n, in the wider of T and double and back to T, the count exact (to 2^53), as NumPy divides by its integer count; a float's count rounded to a float past 2^24 (2^24 + 1 ones averaged to 1, NumPy's 0.99999994).
template<std::floating_point T, typename Order >
T	numerical::detail::settled_zero (T value, const T *first, std::size_t count, std::size_t stride, const Order &larger)
	The extreme value of count elements stride apart from first, its zero's sign settled as IEEE 754-2019's maximum and minimum settle it, and NumPy's on arm64: no zero is larger than the other, so the first met stayed; max gives +0 where one is among them, min −0, whatever their order. A NaN, or any other value, as it is.
std::vector< std::size_t >	numerical::detail::without_axis (std::vector< std::size_t > shape, std::size_t axis)
	The shape with axis axis removed.
template<typename T , typename Reduce >
tensor< T >	numerical::detail::reduce_along (const tensor< T > &t, std::size_t axis, const Reduce &reduce)
	Reduces each line of t along axis by reduce, which takes the line's first element, its length and its stride.
template<typename T >
T	numerical::detail::square (const T &x)
	x · x.
template<typename U >
std::complex< U >	numerical::detail::square (const std::complex< U > &x)
	NumPy's complex square, (re·re − im·im) + (re·im + im·re)i, without the infinity handling of a general product (C Annex G), which costs a call per element.
template<typename T >
bool	numerical::detail::is_nan (const T &)
	A NaN? Never for a scalar without one.
bool	numerical::detail::is_nan (double x)
	A NaN?
template<typename T , typename Order >
T	numerical::detail::contiguous_extreme (const T *first, std::size_t count, const Order &larger)
	The largest (larger is std::greater) or smallest of count contiguous elements, or a NaN if one is among them, as NumPy's: the comparisons run without a branch over native vectors (or eight scalar lanes), each lane its own running best, and the NaN is tracked apart.
template<typename T , typename Order >
T	numerical::detail::extreme (const T *first, std::size_t count, std::size_t stride, const Order &larger)
	The largest (larger is std::greater) or smallest element of a line.
template<typename T , typename Order >
tensor< T >	numerical::detail::extremes_along (const tensor< T > &t, std::size_t axis, const Order &larger)
	The largest (larger is std::greater) or smallest elements along axis. Along an inner axis each line is contiguous (extreme); along an outer one the slabs are compared element by element, so the loop runs over contiguous memory rather than down a column.
std::vector< std::size_t >	numerical::broadcast_shape (const std::vector< std::size_t > &a, const std::vector< std::size_t > &b)
	The shape two shapes broadcast to, as NumPy's broadcast_shapes.
template<typename T , typename Op >
tensor< T >	numerical::broadcast_apply (const tensor< T > &a, const tensor< T > &b, const Op &op)
	op(a, b) element by element, the operands broadcast against each other.
template<detail::split_complex C>
tensor< C >	numerical::numpy_multiply (const tensor< C > &a, const tensor< C > &b)
	a · b of complex tensors, broadcast, each product as NumPy's array loop computes it (detail::numpy_multiply on this processor): the binding's *, where std::complex's operator (this library's *) recovers infinities NumPy's makes NaN.
template<detail::split_complex C>
tensor< C >	numerical::numpy_divide (const tensor< C > &a, const tensor< C > &b)
	a / b of complex tensors, broadcast, each quotient as NumPy's array loop computes it (detail::numpy_divide on this processor): the binding's /.
template<typename T >
tensor< T >	numerical::operator+ (const tensor< T > &a, const tensor< T > &b)
	Element-wise sum, broadcast.
template<typename T >
tensor< T >	numerical::operator- (const tensor< T > &a, const tensor< T > &b)
	Element-wise difference, broadcast.
template<typename T >
tensor< T >	numerical::operator* (const tensor< T > &a, const tensor< T > &b)
	Element-wise product, broadcast.
template<typename T >
tensor< T >	numerical::operator/ (const tensor< T > &a, const tensor< T > &b)
	Element-wise quotient, broadcast.
template<typename T >
tensor< T >	numerical::operator- (const tensor< T > &a)
	Element-wise negation.
template<typename T , typename F >
auto	numerical::map (const tensor< T > &t, const F &f) -> tensor< std::decay_t< decltype(f(std::declval< const T & >()))> >
	f(x) for every element, in a tensor of the same shape and of f's result type.
template<typename T >
tensor< T >	numerical::pow (const tensor< T > &a, const tensor< T > &b)
	a to the power b element by element, broadcast, as std::pow; a single exponent of 2, 0.5 or -1 is computed as NumPy computes it, x·x, sqrt(x) or 1/x (exact where pow rounds: a complex square has no stray real part).
template<typename T >
T	numerical::sum (const tensor< T > &t)
	The sum of every element, by pairwise halving (0 for an empty tensor).
template<typename T >
tensor< T >	numerical::sum (const tensor< T > &t, std::size_t axis)
	The sums along axis, by pairwise halving.
template<typename T >
T	numerical::mean (const tensor< T > &t)
	The mean of every element. An empty tensor's mean is NaN, as NumPy's (which also warns).
template<typename T >
tensor< T >	numerical::mean (const tensor< T > &t, std::size_t axis)
	The means along axis.
template<std::totally_ordered T>
T	numerical::max (const tensor< T > &t)
	The largest element (a NaN, if there is one).
template<std::totally_ordered T>
T	numerical::min (const tensor< T > &t)
	The smallest element (a NaN, if there is one).
template<std::totally_ordered T>
tensor< T >	numerical::max (const tensor< T > &t, std::size_t axis)
	The largest elements along axis.
template<std::totally_ordered T>
tensor< T >	numerical::min (const tensor< T > &t, std::size_t axis)
	The smallest elements along axis.
template<typename T >
tensor< T >	numerical::reshape (const tensor< T > &t, std::vector< std::size_t > shape)
	The same elements in the same order under another shape.
template<typename T >
tensor< T >	numerical::full (std::vector< std::size_t > shape, const T &value)
	A tensor of shape with every element value.
template<typename T >
tensor< T >	numerical::zeros (std::vector< std::size_t > shape)
	A tensor of shape full of zeros.
template<typename T >
tensor< T >	numerical::ones (std::vector< std::size_t > shape)
	A tensor of shape full of ones.
template<typename T >
tensor< T >	numerical::eye (std::size_t n)
	The n × n identity.
tensor< double >	numerical::arange (double start, double stop, double step)
	start, start + step, … below stop (above it for a negative step), as NumPy's arange: its length is ceil((stop − start) / step) and element i is start + i·delta with delta = (start + step) − start, the step as the first two elements hold it, rounded once (a fused multiply-add, as NumPy's arm64 builds compute it; one without fusion rounds the product first and can differ in the last bits).
tensor< double >	numerical::linspace (double start, double stop, std::size_t count)
	count values evenly spaced from start to stop, both included, as NumPy's linspace: element i is start + i·step with step = (stop − start) / (count − 1), and the last is stop exactly.

Detailed Description

Element-wise arithmetic with NumPy's broadcasting, reductions, reshaping and creation of tensors.

Two shapes broadcast as NumPy's do: aligned on their last axes, each pair of extents equal or one of them 1 (a missing leading axis counts as 1), the result taking the larger. A sum adds by pairwise halving, as NumPy's reference algorithm does (NumPy's vector builds group the terms their own way, so the last bits can differ by platform): its rounding error grows with the logarithm of the count instead of the count. A maximum or a minimum is defined for an ordered scalar only, and an empty reduction without an identity throws.

Definition in file tensor_ops.hpp.