rustymath_core/float.rs
1//! The float types that rustymath's functions take, and the operations that its algorithms build on.
2use crate::search::{Key, RadixKey};
3use crate::{parallel::Slots, slice::{Arg, Ternary}};
4
5/// 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.
6pub const FUSED: bool = cfg!(any(target_arch = "aarch64", target_feature = "fma"));
7
8/// A float type that rustymath's functions take: `half::f16`, `f32` or `f64`.
9pub trait Float: num_traits::Float + sealed::Sealed + crate::slice::KernelItem + crate::buffer::Source<Self> {
10 /// 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.
11 type Working: Native;
12 /// This float in its working type, which holds it exactly.
13 fn working(self) -> Self::Working;
14 /// `t` rounded to this type.
15 fn from_working(t: Self::Working) -> Self;
16 /// This float as an `f64`, which holds it exactly.
17 fn wide(self) -> f64;
18 /// `x` rounded to the nearest float of this type.
19 fn narrow(x: f64) -> Self;
20 /// The least float of this type above `self`, as `f64::next_up` gives it.
21 fn next_up(self) -> Self;
22 /// The greatest float of this type below `self`, as `f64::next_down` gives it.
23 fn next_down(self) -> Self;
24 /// Rounds to the nearest whole number, with ties to even.
25 fn round_ties_even(self) -> Self;
26 /// The unsigned type, of this float's width, that `order_key` gives.
27 type OrderKey: RadixKey + Key;
28 /// 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.
29 fn order_key(self) -> Self::OrderKey;
30 /// 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.
31 fn total_key(self) -> Self::OrderKey;
32 /// 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.
33 fn from_order_key(k: Self::OrderKey) -> Self;
34 /// 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.
35 fn choose(x: Arg<'_, Self>, y: Arg<'_, Self>, mask: Arg<'_, Self>, n: usize, out: &mut Slots<'_, Self>) -> bool;
36}
37
38/// A float type that rustymath's algorithms run in: `f32` or `f64`.
39pub trait Native: Float<Working = Self> + sealed::Core + sealed::Folds {
40 /// `+` that may be reordered with other algebraic operations. A sum of them vectorizes.
41 fn algebraic_add(self, other: Self) -> Self;
42 /// `×` that may be reordered with other algebraic operations.
43 fn algebraic_mul(self, other: Self) -> Self;
44}
45
46/// Traits that only rustymath implements. `Core` holds the bit-level operations that the algorithms build on.
47pub(crate) mod sealed {
48 use num_complex::Complex;
49 pub trait Sealed {}
50 macro_rules! sealed { ($($t:ty),+) => {$( impl Sealed for $t {} )+} }
51 sealed!(i8, i16, i32, i64, u8, u16, u32);
52 impl<T: super::Native> Sealed for Complex<T> {}
53
54 pub trait Core: Copy {
55 /// Whether this is `f64`. Algorithms pick their constants by it, and the compiler removes the branch.
56 const WIDE: bool;
57 /// 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.
58 const SHIFT: Self;
59 /// The unsigned integer type with this float's bits.
60 type Bits: Copy;
61 /// `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.
62 fn from_low(u: Self::Bits, less: Self) -> Self;
63 /// `self×b+c`, with one rounding when `FUSED` is true and two otherwise.
64 fn mla<const FUSED: bool>(self, b: Self, c: Self) -> Self;
65 /// 2^n, from `z` = `SHIFT + n`, for a whole n within the exponents of normal floats.
66 fn pow2(z: Self) -> Self;
67 /// `(m, e)` with `self` = m × 2^e, m from 0.75 to 1.5 and e whole, for a positive normal `self`.
68 fn split(self) -> (Self, Self);
69 /// 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.
70 fn positive_normal(self) -> bool;
71 /// `self` negated when the whole number n is odd, for `z` = `SHIFT + n`.
72 fn flip(self, z: Self) -> Self;
73 /// Whether the whole number n is odd, for `self` = `SHIFT + n`.
74 fn odd(self) -> bool;
75 /// `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.
76 fn upper(self) -> Self;
77 /// `(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.
78 fn scales(self) -> (Self, Self);
79 /// 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.
80 fn maximum(self, y: Self) -> Self;
81 /// The lesser of `self` and `y`, as `maximum` gives the greater.
82 fn minimum(self, y: Self) -> Self;
83 }
84 /// 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.
85 pub trait Folds: Core + num_traits::Float {
86 /// 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.
87 #[inline(always)]
88 fn step_max(t: Self, v: Self) -> Self { Self::maximum(v, t) }
89 /// One step of a fold or scan by minimum, as `step_max` is one by maximum.
90 #[inline(always)]
91 fn step_min(t: Self, v: Self) -> Self { Self::minimum(v, t) }
92 fn fold_max(x: &[Self]) -> Self { crate::reduce::folded(x, Self::neg_infinity(), Self::step_max, Self::maximum) }
93 fn fold_min(x: &[Self]) -> Self { crate::reduce::folded(x, Self::infinity(), Self::step_min, Self::minimum) }
94 fn scan_max(x: &[Self], start: Option<Self>, out: &mut Vec<Self>) {
95 crate::reduce::records(x, start, out, |v| Self::step_max(Self::neg_infinity(), v), Self::step_max, crate::reduce::below)
96 }
97 fn scan_min(x: &[Self], start: Option<Self>, out: &mut Vec<Self>) {
98 crate::reduce::records(x, start, out, |v| Self::step_min(Self::infinity(), v), Self::step_min, crate::reduce::above)
99 }
100 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) }
101 /// 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.
102 #[inline(always)]
103 fn ordered_lanes<const W: usize>(a: [Self; W], c: [Self; W]) -> ([Self; W], [Self; W]) {
104 (std::array::from_fn(|j| c[j].minimum(a[j])), std::array::from_fn(|j| a[j].maximum(c[j])))
105 }
106 }
107}
108
109/// 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.
110#[inline(always)]
111fn ordered_bits<K: num_traits::PrimInt>(bits: K) -> K {
112 let top = K::zero().count_zeros() - 1;
113 // Loops that compute keys vectorize only without a branch here. The arithmetic shift spreads the sign bit over every bit.
114 bits ^ (bits.signed_shr(top) | K::one() << top as usize)
115}
116
117/// `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 `∞`.
118macro_rules! order_key {
119 ($k:ty) => {
120 type OrderKey = $k;
121 #[inline(always)]
122 fn order_key(self) -> $k {
123 let key = ordered_bits((self + <Self as num_traits::Zero>::zero()).to_bits());
124 if self.is_nan() { <$k>::MAX } else { key }
125 }
126 #[inline(always)]
127 fn total_key(self) -> $k {
128 let bits = self.to_bits();
129 ordered_bits(if self.is_nan() { bits & !(1 << (<$k>::BITS - 1)) } else { bits })
130 }
131 #[inline(always)]
132 fn from_order_key(k: $k) -> Self {
133 Self::from_bits(if k >> (<$k>::BITS - 1) == 1 { k & !(1 << (<$k>::BITS - 1)) } else { !k })
134 }
135 };
136}
137
138macro_rules! native {
139 ($t:ident, $bits:ty, $mantissa:literal, $bias:literal) => {
140 impl sealed::Sealed for $t {}
141 impl Float for $t {
142 type Working = $t;
143 #[inline]
144 fn working(self) -> $t { self }
145 #[inline]
146 fn from_working(t: $t) -> $t { t }
147 #[inline]
148 fn wide(self) -> f64 { self as f64 }
149 #[inline]
150 fn narrow(x: f64) -> $t { x as $t }
151 #[inline]
152 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) }
153 #[inline]
154 fn next_up(self) -> $t { $t::next_up(self) }
155 #[inline]
156 fn next_down(self) -> $t { $t::next_down(self) }
157 #[inline]
158 fn round_ties_even(self) -> Self { $t::round_ties_even(self) }
159 order_key!($bits);
160 }
161 impl Native for $t {
162 #[inline]
163 fn algebraic_add(self, other: $t) -> $t { $t::algebraic_add(self, other) }
164 #[inline]
165 fn algebraic_mul(self, other: $t) -> $t { $t::algebraic_mul(self, other) }
166 }
167 impl sealed::Core for $t {
168 const WIDE: bool = $mantissa == 52;
169 const SHIFT: $t = 1.5 * (1u64 << $mantissa) as $t;
170 type Bits = $bits;
171 #[inline]
172 fn from_low(u: $bits, less: $t) -> $t {
173 let two_m = (1 as $bits << $mantissa) as $t;
174 $t::from_bits(two_m.to_bits() | u) - (two_m + less)
175 }
176 #[inline]
177 fn mla<const FUSED: bool>(self, b: $t, c: $t) -> $t { if FUSED { self.mul_add(b, c) } else { self * b + c } }
178 // 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.
179 #[inline]
180 fn pow2(z: $t) -> $t { crate::recipe::pow2(&crate::recipe::Floats, z) }
181 #[inline]
182 fn split(self) -> ($t, $t) {
183 crate::recipe::split(&crate::recipe::Floats, self)
184 }
185 // 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.
186 #[inline(always)]
187 fn positive_normal(self) -> bool {
188 const SHIFT: u32 = <$bits>::BITS - 32;
189 ((self.to_bits() >> SHIFT) as u32).wrapping_sub(1 << ($mantissa - SHIFT)) < (2 * $bias) << ($mantissa - SHIFT)
190 }
191 // The lowest bit of `z` is the parity of n. Shifted into the sign bit, it negates `self` with no branch.
192 #[inline]
193 fn flip(self, z: $t) -> $t { $t::from_bits(self.to_bits() ^ (z.to_bits() << (<$bits>::BITS - 1))) }
194 #[inline]
195 fn odd(self) -> bool { self.to_bits() & 1 == 1 }
196 #[inline]
197 fn upper(self) -> $t { $t::from_bits(self.to_bits() & !((1 << (($mantissa + 2) / 2)) - 1)) }
198 #[inline]
199 fn scales(self) -> ($t, $t) {
200 let bottom = $t::from_bits(2 << $mantissa);
201 let e = (if self < bottom { bottom } else { self }).to_bits() >> $mantissa;
202 ($t::from_bits((2 * $bias + 1 - e) << $mantissa), $t::from_bits((e - 1) << $mantissa))
203 }
204 #[inline(always)]
205 fn maximum(self, y: $t) -> $t { if self > y { self } else { y } }
206 #[inline(always)]
207 fn minimum(self, y: $t) -> $t { if self < y { self } else { y } }
208 }
209 };
210}
211native!(f32, u32, 23, 127);
212native!(f64, u64, 52, 1023);
213
214impl sealed::Sealed for half::f16 {}
215impl Float for half::f16 {
216 type Working = f32;
217 #[inline]
218 fn working(self) -> f32 { self.to_f32() }
219 #[inline]
220 fn from_working(t: f32) -> Self { half::f16::from_f32(t) }
221 #[inline]
222 fn wide(self) -> f64 { self.to_f64() }
223 #[inline]
224 fn narrow(x: f64) -> Self { half::f16::from_f64(x) }
225 #[inline]
226 fn choose(x: Arg<'_, Self>, y: Arg<'_, Self>, mask: Arg<'_, Self>, n: usize, out: &mut Slots<'_, Self>) -> bool {
227 crate::Choose.many_into(x, y, mask, n, out)
228 }
229 // 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.
230 #[inline]
231 fn next_up(self) -> Self {
232 let bits = self.to_bits();
233 let next = if bits & 0x7fff == 0 { 1 } else if bits >> 15 == 0 { bits + 1 } else { bits - 1 };
234 if self.is_nan() || self == half::f16::INFINITY { self } else { half::f16::from_bits(next) }
235 }
236 #[inline]
237 fn next_down(self) -> Self { -(-self).next_up() }
238 #[inline]
239 fn round_ties_even(self) -> Self { half::f16::from_f32(self.to_f32().round_ties_even()) }
240 order_key!(u16);
241}
242
243/// 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.
244#[inline]
245pub(crate) fn horner<R: Native, const FUSED: bool, const N: usize>(x: R, c: [f64; N]) -> R {
246 crate::recipe::horner(&crate::recipe::NativeOps::<FUSED>::DEFAULT, x, &c)
247}
248
249/// 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.
250#[inline]
251pub(crate) fn ln2<R: Native>() -> (R, R) { crate::recipe::ln2(&crate::recipe::Floats) }
252
253/// 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.
254#[inline]
255pub(crate) fn pow2_times<R: Native, const FUSED: bool>(y: R, n: R) -> R { crate::recipe::ldexp(&crate::recipe::NativeOps::<FUSED>::DEFAULT, y, n) }
256
257/// (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.
258#[inline]
259pub(crate) fn split_positive<R: Native>(x: R) -> (R, R) { crate::recipe::frexp(&crate::recipe::Floats, x) }