//! The float types that rustymath's functions take, and the operations that its algorithms build on.
use crate::search::{Key, RadixKey};
use crate::{parallel::Slots, slice::{Arg, Ternary}};
/// Whether `mla` fuses its multiply and add into one operation with one rounding. It does where the target has a fused multiply-add instruction, as on AArch64 and on x86-64 with FMA. Elsewhere, as in wasm and on baseline x86-64, `mul_add` calls a software `fma`, which is slow and stops the loop around it vectorizing. There `mla` multiplies and adds with two roundings.
pub const FUSED: bool = cfg!(any(target_arch = "aarch64", target_feature = "fma"));
/// A float type that rustymath's functions take: `half::f16`, `f32` or `f64`.
pub trait Float: num_traits::Float + sealed::Sealed + crate::slice::KernelItem + crate::buffer::Source<Self> {
/// The type this float computes in, its working precision: `f32` for `f16`, and the type itself otherwise. A function of an `f16` computes in `f32` and rounds once. Sums keep their running totals in this type.
type Working: Native;
/// This float in its working type, which holds it exactly.
fn working(self) -> Self::Working;
/// `t` rounded to this type.
fn from_working(t: Self::Working) -> Self;
/// This float as an `f64`, which holds it exactly.
fn wide(self) -> f64;
/// `x` rounded to the nearest float of this type.
fn narrow(x: f64) -> Self;
/// The least float of this type above `self`, as `f64::next_up` gives it.
fn next_up(self) -> Self;
/// The greatest float of this type below `self`, as `f64::next_down` gives it.
fn next_down(self) -> Self;
/// Rounds to the nearest whole number, with ties to even.
fn round_ties_even(self) -> Self;
/// The unsigned type, of this float's width, that `order_key` gives.
type OrderKey: RadixKey + Key;
/// A key whose unsigned order is the order of `total_cmp`, after `¯0` becomes `0` and every NaN becomes one NaN. So `¯0` and `0` share a key, NaN's key follows every other, and two floats match exactly when their keys are equal.
fn order_key(self) -> Self::OrderKey;
/// A key whose unsigned order is the order of `total_cmp`, after every NaN becomes positive. So `¯0`'s key comes just before `0`'s, every NaN's follows `∞`'s, and `from_order_key` reads the float back exactly, apart from a NaN's sign. Sorts and selections compare these keys and write floats back from them.
fn total_key(self) -> Self::OrderKey;
/// The float whose `order_key` or `total_key` is `k`. A total key reads back exactly, apart from a NaN's sign. An order key reads `¯0` back as `0`, and every NaN as one NaN.
fn from_order_key(k: Self::OrderKey) -> Self;
/// Writes `x` where `mask`'s lanes aren't zero and `y` where they are into the next `n` slots of `out`, through [`Choose`](crate::Choose), and gives whether every result succeeded. Fused graphs pass a comparison's lanes of 1 and 0 as the mask.
fn choose(x: Arg<'_, Self>, y: Arg<'_, Self>, mask: Arg<'_, Self>, n: usize, out: &mut Slots<'_, Self>) -> bool;
}
/// A float type that rustymath's algorithms run in: `f32` or `f64`.
pub trait Native: Float<Working = Self> + sealed::Core + sealed::Folds {
/// `+` that may be reordered with other algebraic operations. A sum of them vectorizes.
fn algebraic_add(self, other: Self) -> Self;
/// `×` that may be reordered with other algebraic operations.
fn algebraic_mul(self, other: Self) -> Self;
}
/// Traits that only rustymath implements. `Core` holds the bit-level operations that the algorithms build on.
pub(crate) mod sealed {
use num_complex::Complex;
pub trait Sealed {}
macro_rules! sealed { ($($t:ty),+) => {$( impl Sealed for $t {} )+} }
sealed!(i8, i16, i32, i64, u8, u16, u32);
impl<T: super::Native> Sealed for Complex<T> {}
pub trait Core: Copy {
/// Whether this is `f64`. Algorithms pick their constants by it, and the compiler removes the branch.
const WIDE: bool;
/// 1.5 × 2^m, where m is the number of mantissa bits. Adding it to a float of magnitude below 2^(m-1) rounds that float to a whole number, which the low bits of the sum hold.
const SHIFT: Self;
/// The unsigned integer type with this float's bits.
type Bits: Copy;
/// `u - less` as a float, for a whole number u below 2^m, through the bits of 2^m + u. A loop that calls it vectorizes even where the target has no vector conversion for u's type. [`convert`](crate::convert) builds on it.
fn from_low(u: Self::Bits, less: Self) -> Self;
/// `self×b+c`, with one rounding when `FUSED` is true and two otherwise.
fn mla<const FUSED: bool>(self, b: Self, c: Self) -> Self;
/// 2^n, from `z` = `SHIFT + n`, for a whole n within the exponents of normal floats.
fn pow2(z: Self) -> Self;
/// `(m, e)` with `self` = m × 2^e, m from 0.75 to 1.5 and e whole, for a positive normal `self`.
fn split(self) -> (Self, Self);
/// Whether `self` is positive, normal and finite, from the top 32 bits of its bits, which hold its sign and exponent, through one subtraction and one unsigned comparison. A loop of them vectorizes on every target, four f64 values to a 128-bit register, where LLVM merges the comparisons `self >= MIN_POSITIVE` and `self < ∞` into a class test that it leaves scalar on x86-64.
fn positive_normal(self) -> bool;
/// `self` negated when the whole number n is odd, for `z` = `SHIFT + n`.
fn flip(self, z: Self) -> Self;
/// Whether the whole number n is odd, for `self` = `SHIFT + n`.
fn odd(self) -> bool;
/// `self` with the low half of its mantissa bits cleared. The product of two such floats is exact, which double-word products without FMA rely on.
fn upper(self) -> Self;
/// `(s, 1/s)` for the power of two s that brings `self` into [2, 4), after raising a smaller `self` to 2^(2 - bias). Both are normal for every finite `self`, and they come from its exponent bits with no branch.
fn scales(self) -> (Self, Self);
/// The greater of `self` and `y`, from a comparison and a select. Where they're equal, or either is NaN, it gives `y`, as x86-64's `maxsd` does. A loop of them vectorizes on every target.
fn maximum(self, y: Self) -> Self;
/// The lesser of `self` and `y`, as `maximum` gives the greater.
fn minimum(self, y: Self) -> Self;
}
/// The folds and scans of floats by maximum and minimum, and their running sums. Every route that folds by a maximum or minimum keeps one rule: no running result is ever NaN. Each fold or scan starts at `-∞` or `∞`. An item joins a running result only through `Core::maximum` or `step_max`, which take the running result as `y` and keep it beside a NaN item. Merges combine running results alone. The rule makes every route give the maximum of the items that aren't NaN, whatever its grouping, piece sizes, vector width or target. The one exception is an AArch64 scan with no start, which reads its first item as it is: a NaN first item is its first running result, and the next `fmaxnm` drops it. Where a step keeps a NaN running result, a scan turns each NaN item into its start as it reads it, so the groups that `grouped` starts at an item keep the rule. The defaults' fold loops vectorize. The default running sum is `grouped`'s. `reduce.rs` overrides them where a target has instructions that do better.
pub trait Folds: Core + num_traits::Float {
/// One step of a fold or scan by maximum: `Core::maximum` of the item `v` and the running result `t`, which gives `t` beside a NaN item. Chains, where each step waits for the one before, take it through `Reduction::chain`. AArch64 overrides it with `fmaxnm`, one instruction to the select's two. LLVM can't vectorize that override, so loops that vectorize call `Core::maximum` instead.
#[inline(always)]
fn step_max(t: Self, v: Self) -> Self { Self::maximum(v, t) }
/// One step of a fold or scan by minimum, as `step_max` is one by maximum.
#[inline(always)]
fn step_min(t: Self, v: Self) -> Self { Self::minimum(v, t) }
fn fold_max(x: &[Self]) -> Self { crate::reduce::folded(x, Self::neg_infinity(), Self::step_max, Self::maximum) }
fn fold_min(x: &[Self]) -> Self { crate::reduce::folded(x, Self::infinity(), Self::step_min, Self::minimum) }
fn scan_max(x: &[Self], start: Option<Self>, out: &mut Vec<Self>) {
crate::reduce::records(x, start, out, |v| Self::step_max(Self::neg_infinity(), v), Self::step_max, crate::reduce::below)
}
fn scan_min(x: &[Self], start: Option<Self>, out: &mut Vec<Self>) {
crate::reduce::records(x, start, out, |v| Self::step_min(Self::infinity(), v), Self::step_min, crate::reduce::above)
}
fn scan_sum(x: &[Self], start: Option<Self>, out: &mut Vec<Self>) { crate::reduce::grouped(x, start, out, |v| v, |a, b| a + b, |t| t, false) }
/// The rows `a` and `c` lane by lane, the lesser first, for items that aren't NaN. Equal items keep their places. The default takes one comparison and two selects, as `Core::minimum` and `Core::maximum` do, which x86-64 makes `minpd` and `maxpd`. `reduce.rs` overrides it where a row fills one vector, because LLVM leaves some of those selects scalar there.
#[inline(always)]
fn ordered_lanes<const W: usize>(a: [Self; W], c: [Self; W]) -> ([Self; W], [Self; W]) {
(std::array::from_fn(|j| c[j].minimum(a[j])), std::array::from_fn(|j| a[j].maximum(c[j])))
}
}
}
/// A float's bits as a key in its total order: the sign bit flipped, and a negative float's other bits flipped too. `Float::from_order_key` undoes it.
#[inline(always)]
fn ordered_bits<K: num_traits::PrimInt>(bits: K) -> K {
let top = K::zero().count_zeros() - 1;
// Loops that compute keys vectorize only without a branch here. The arithmetic shift spreads the sign bit over every bit.
bits ^ (bits.signed_shr(top) | K::one() << top as usize)
}
/// `Float`'s order-key items, for a float whose bits have type `$k`, through `ordered_bits`. An order key adds `¯0` to `0` first and gives every NaN the greatest key. A total key clears a NaN's sign bit, so every NaN follows `∞`.
macro_rules! order_key {
($k:ty) => {
type OrderKey = $k;
#[inline(always)]
fn order_key(self) -> $k {
let key = ordered_bits((self + <Self as num_traits::Zero>::zero()).to_bits());
if self.is_nan() { <$k>::MAX } else { key }
}
#[inline(always)]
fn total_key(self) -> $k {
let bits = self.to_bits();
ordered_bits(if self.is_nan() { bits & !(1 << (<$k>::BITS - 1)) } else { bits })
}
#[inline(always)]
fn from_order_key(k: $k) -> Self {
Self::from_bits(if k >> (<$k>::BITS - 1) == 1 { k & !(1 << (<$k>::BITS - 1)) } else { !k })
}
};
}
macro_rules! native {
($t:ident, $bits:ty, $mantissa:literal, $bias:literal) => {
impl sealed::Sealed for $t {}
impl Float for $t {
type Working = $t;
#[inline]
fn working(self) -> $t { self }
#[inline]
fn from_working(t: $t) -> $t { t }
#[inline]
fn wide(self) -> f64 { self as f64 }
#[inline]
fn narrow(x: f64) -> $t { x as $t }
#[inline]
fn choose(x: Arg<'_, $t>, y: Arg<'_, $t>, mask: Arg<'_, $t>, n: usize, out: &mut Slots<'_, $t>) -> bool { crate::Choose.many_into(x, y, mask, n, out) }
#[inline]
fn next_up(self) -> $t { $t::next_up(self) }
#[inline]
fn next_down(self) -> $t { $t::next_down(self) }
#[inline]
fn round_ties_even(self) -> Self { $t::round_ties_even(self) }
order_key!($bits);
}
impl Native for $t {
#[inline]
fn algebraic_add(self, other: $t) -> $t { $t::algebraic_add(self, other) }
#[inline]
fn algebraic_mul(self, other: $t) -> $t { $t::algebraic_mul(self, other) }
}
impl sealed::Core for $t {
const WIDE: bool = $mantissa == 52;
const SHIFT: $t = 1.5 * (1u64 << $mantissa) as $t;
type Bits = $bits;
#[inline]
fn from_low(u: $bits, less: $t) -> $t {
let two_m = (1 as $bits << $mantissa) as $t;
$t::from_bits(two_m.to_bits() | u) - (two_m + less)
}
#[inline]
fn mla<const FUSED: bool>(self, b: $t, c: $t) -> $t { if FUSED { self.mul_add(b, c) } else { self * b + c } }
// The low bits of `z` hold n. Adding the bias to them and shifting them into the exponent field gives 2^n. Integer adds and shifts vectorize on every target, where a conversion from a float to a 64-bit integer doesn't on x86-64.
#[inline]
fn pow2(z: $t) -> $t { crate::recipe::pow2(&crate::recipe::Floats, z) }
#[inline]
fn split(self) -> ($t, $t) {
crate::recipe::split(&crate::recipe::Floats, self)
}
// Subtracting the smallest normal float's top bits wraps zero and the subnormals around to the top, where negative floats' top bits already lie. Only positive normal floats then lie below infinity's top bits less the same amount.
#[inline(always)]
fn positive_normal(self) -> bool {
const SHIFT: u32 = <$bits>::BITS - 32;
((self.to_bits() >> SHIFT) as u32).wrapping_sub(1 << ($mantissa - SHIFT)) < (2 * $bias) << ($mantissa - SHIFT)
}
// The lowest bit of `z` is the parity of n. Shifted into the sign bit, it negates `self` with no branch.
#[inline]
fn flip(self, z: $t) -> $t { $t::from_bits(self.to_bits() ^ (z.to_bits() << (<$bits>::BITS - 1))) }
#[inline]
fn odd(self) -> bool { self.to_bits() & 1 == 1 }
#[inline]
fn upper(self) -> $t { $t::from_bits(self.to_bits() & !((1 << (($mantissa + 2) / 2)) - 1)) }
#[inline]
fn scales(self) -> ($t, $t) {
let bottom = $t::from_bits(2 << $mantissa);
let e = (if self < bottom { bottom } else { self }).to_bits() >> $mantissa;
($t::from_bits((2 * $bias + 1 - e) << $mantissa), $t::from_bits((e - 1) << $mantissa))
}
#[inline(always)]
fn maximum(self, y: $t) -> $t { if self > y { self } else { y } }
#[inline(always)]
fn minimum(self, y: $t) -> $t { if self < y { self } else { y } }
}
};
}
native!(f32, u32, 23, 127);
native!(f64, u64, 52, 1023);
impl sealed::Sealed for half::f16 {}
impl Float for half::f16 {
type Working = f32;
#[inline]
fn working(self) -> f32 { self.to_f32() }
#[inline]
fn from_working(t: f32) -> Self { half::f16::from_f32(t) }
#[inline]
fn wide(self) -> f64 { self.to_f64() }
#[inline]
fn narrow(x: f64) -> Self { half::f16::from_f64(x) }
#[inline]
fn choose(x: Arg<'_, Self>, y: Arg<'_, Self>, mask: Arg<'_, Self>, n: usize, out: &mut Slots<'_, Self>) -> bool {
crate::Choose.many_into(x, y, mask, n, out)
}
// The float above a positive float has the next bit pattern. The float above a negative float has the previous one. Both zeros step to the least positive subnormal.
#[inline]
fn next_up(self) -> Self {
let bits = self.to_bits();
let next = if bits & 0x7fff == 0 { 1 } else if bits >> 15 == 0 { bits + 1 } else { bits - 1 };
if self.is_nan() || self == half::f16::INFINITY { self } else { half::f16::from_bits(next) }
}
#[inline]
fn next_down(self) -> Self { -(-self).next_up() }
#[inline]
fn round_ties_even(self) -> Self { half::f16::from_f32(self.to_f32().round_ties_even()) }
order_key!(u16);
}
/// The polynomial with coefficients `c`, highest degree first, at `x`, by Horner's rule. The coefficients convert to `R` when the compiler inlines the call.
#[inline]
pub(crate) fn horner<R: Native, const FUSED: bool, const N: usize>(x: R, c: [f64; N]) -> R {
crate::recipe::horner(&crate::recipe::NativeOps::<FUSED>::DEFAULT, x, &c)
}
/// ln 2, split into a part whose product with any whole number up to 2^11 at f64 or 2^8 at f32 is exact, and the rest. That covers every exponent of a float and every n that `exp` reduces by.
#[inline]
pub(crate) fn ln2<R: Native>() -> (R, R) { crate::recipe::ln2(&crate::recipe::Floats) }
/// y × 2^n, for a whole n, in two factors. The first multiplication is exact, so a result below the normal floats rounds once, in the second.
#[inline]
pub(crate) fn pow2_times<R: Native, const FUSED: bool>(y: R, n: R) -> R { crate::recipe::ldexp(&crate::recipe::NativeOps::<FUSED>::DEFAULT, y, n) }
/// (m, e) with x = m × 2^e, m from 0.75 to 1.5 and e whole, for a positive x, subnormals included. A subnormal x scales up by 2^64 first, so that its bits split as a normal float's do.
#[inline]
pub(crate) fn split_positive<R: Native>(x: R) -> (R, R) { crate::recipe::frexp(&crate::recipe::Floats, x) }