Skip to main content

thermite/element/float/
mod.rs

1use super::{SignedElement, SignedIntegerElement, UnsignedIntegerElement};
2use crate::LargeInt;
3// Named only by the `std` arm of the `cfg_if!` below.
4#[cfg(feature = "std")]
5use crate::register::FloatRegister;
6use crate::vector::SplatConst;
7use crate::vector::ops::MulAddExt;
8
9pub mod spec;
10
11pub(crate) mod algebraic;
12
13/// Marker type for a compile-time integer constant cast to a float element type.
14///
15/// Implements [`SplatConst<f32>`] and [`SplatConst<f64>`], enabling use with
16/// [`const_splat!`](crate::const_splat) and [`FloatElement::ConstInt`].
17pub struct IntConst<const N: crate::LargeInt>;
18
19/// Marker type for a compile-time rational constant (N/D) cast to a float element type.
20///
21/// Implements [`SplatConst<f32>`] and [`SplatConst<f64>`], enabling use with
22/// [`const_splat!`](crate::const_splat) and [`FloatElement::ConstRatio`].
23pub struct RatioConst<const N: crate::LargeInt, const D: crate::LargeInt>;
24
25impl<const N: crate::LargeInt> SplatConst<f32> for IntConst<N> {
26    const VALUE: f32 = N as f32;
27}
28
29impl<const N: crate::LargeInt> SplatConst<f64> for IntConst<N> {
30    const VALUE: f64 = N as f64;
31}
32
33impl<const N: crate::LargeInt, const D: crate::LargeInt> SplatConst<f32> for RatioConst<N, D> {
34    const VALUE: f32 = {
35        assert!(D != 0, "RatioConst: denominator must not be zero");
36        let (q, r) = (N / D, N % D);
37        (q as f32) + (r as f32) / (D as f32)
38    };
39}
40
41impl<const N: crate::LargeInt, const D: crate::LargeInt> SplatConst<f64> for RatioConst<N, D> {
42    const VALUE: f64 = {
43        assert!(D != 0, "RatioConst: denominator must not be zero");
44        let (q, r) = (N / D, N % D);
45        (q as f64) + (r as f64) / (D as f64)
46    };
47}
48
49/// A trait for float element types that can be used in SIMD operations.
50///
51/// This provides common scalar fallbacks, as well as float specifications for
52/// non-IEE 754 floating point formats
53pub trait FloatElement:
54    SignedElement
55    + crate::math::FloatConsts
56    + num_traits::NumOps
57    + core::ops::Neg<Output = Self>
58    + MulAddExt<Self, Self, Output = Self>
59{
60    /// Marker type for splatting a compile-time integer constant as this float type.
61    ///
62    /// Satisfies `SplatConst<Self>`, enabling const-folded splats via
63    /// [`const_splat!`](crate::const_splat).
64    type ConstInt<const N: crate::LargeInt>: SplatConst<Self>;
65
66    /// Marker type for splatting a compile-time rational constant (N/D) as this float type.
67    ///
68    /// Satisfies `SplatConst<Self>`, enabling const-folded splats via
69    /// [`const_splat!`](crate::const_splat).
70    type ConstRatio<const N: crate::LargeInt, const D: crate::LargeInt>: SplatConst<Self>;
71
72    /// Try to represent this LargeInt value as this float type,
73    /// returning None if it cannot be represented exactly.
74    fn try_from_int(value: LargeInt) -> Option<Self>;
75    fn try_from_ratio(n: LargeInt, d: LargeInt) -> Option<Self>;
76
77    cfg_if::cfg_if! {
78        if #[cfg(all(feature = "spirv", target_arch = "spirv"))] {
79            #[inline(always)]
80            fn from_int(value: LargeInt) -> Self {
81                Self::try_from_int(value).unwrap_or(Self::ZERO)
82            }
83
84            #[inline(always)]
85            fn from_ratio(n: LargeInt, d: LargeInt) -> Self {
86                Self::try_from_ratio(n, d).unwrap_or(Self::ZERO)
87            }
88        } else {
89            #[inline(always)]
90            fn from_int(value: LargeInt) -> Self {
91                #[cold]
92                fn _panic_int_overflow() -> ! {
93                    panic!("LargeInt value exceeds maximum exact representable value for this float type")
94                }
95
96                Self::try_from_int(value).unwrap_or_else(|| _panic_int_overflow())
97            }
98
99            #[inline(always)]
100            fn from_ratio(n: LargeInt, d: LargeInt) -> Self {
101                #[cold]
102                fn _panic_ratio_overflow() -> ! {
103                    panic!("LargeInt ratio exceeds maximum exact representable value for this float type")
104                }
105
106                Self::try_from_ratio(n, d).unwrap_or_else(|| _panic_ratio_overflow())
107            }
108        }
109    }
110
111    fn sqrt(value: Self) -> Self;
112    fn floor(value: Self) -> Self;
113    fn ceil(value: Self) -> Self;
114    fn round(value: Self) -> Self;
115    fn trunc(value: Self) -> Self;
116
117    #[inline(always)]
118    fn fract(value: Self) -> Self {
119        value - FloatElement::trunc(value) // fallback implementation
120    }
121
122    fn next_up(value: Self) -> Self;
123    fn next_down(value: Self) -> Self;
124
125    /// Does the format support Infinity?
126    /// If FALSE, overflow saturates to MAX_FINITE instead of INF.
127    /// (e.g., E4M3 = false, E5M2 = true)
128    const HAS_INFINITY: bool;
129
130    /// Does the format distinguish between +0 and -0?
131    /// (Usually true, but some integer-like quantizations might not)
132    const HAS_SIGNED_ZERO: bool;
133
134    /// Does the format support subnormal numbers?
135    /// If FALSE, any value smaller than MinNormal is flushed to zero (FTZ).
136    const HAS_SUBNORMALS: bool;
137}
138
139pub trait FloatElementWithBits: FloatElement {
140    type Bits: UnsignedIntegerElement<Unsigned = Self::Bits> + TryFrom<u32>;
141    type SignedBits: SignedIntegerElement<Signed = Self::SignedBits> + TryFrom<u32>;
142
143    const EXP_BITS: u32;
144    const MANTISSA_BITS: u32;
145    const EXP_BIAS: Self::SignedBits;
146
147    /// The specific bit pattern for NaN.
148    /// IEEE formats have a *range* of NaNs, but E4M3 has only *one* (0x7F).
149    const NAN_PATTERN: Option<Self::Bits>;
150
151    /// If !HAS_INFINITY, what is the max finite bit pattern?
152    /// Used for clamping overflow.
153    const MAX_FINITE_PATTERN: Self::Bits;
154
155    /// Largest positive subnormal value
156    const MAX_SUBNORMAL: Self::Bits;
157
158    /// Magic value for crushing denormals
159    const DENORMAL_TRICK: Self::Bits;
160
161    // /// Is there an implicit leading bit (1.xxx)?
162    // /// Almost always TRUE.
163    // /// Exception: x87 80-bit float (FALSE).
164    // const IMPLICIT_LEAD_BIT: bool = true;
165
166    // maximum unsigned integer that can be exactly represented in this float type without loss of precision
167    const MAX_LARGE_UINT: crate::LargeUInt;
168
169    const MAX_BIASED_EXP: Self::SignedBits;
170    const EXP_LSB_MASK: Self::Bits;
171    const SIGN_MANTISSA_MASK: Self::Bits;
172
173    const HALF_EXP_BITS: Self::Bits;
174    const FREXP_BIAS_OFFSET: Self::SignedBits;
175
176    /// Convert from f64 to this float type, potentially losing precision.
177    fn from_f64(value: f64) -> Self;
178
179    fn from_signed(value: Self::SignedBits) -> Self;
180}
181
182trait FloatElementInternal: FloatElement {
183    fn try_from_int(value: crate::LargeInt) -> Option<Self>;
184    fn try_from_ratio(n: crate::LargeInt, d: crate::LargeInt) -> Option<Self>;
185}
186
187macro_rules! impl_float_element {
188    (CONSTS $($const:ident: $const_ty:ty = $value:expr;)+) => {paste::paste! {
189        $(const $const: $const_ty = $value;)+
190
191        const FREXP_BIAS_OFFSET: Self::SignedBits = Self::EXP_BIAS - 1;
192        const HALF_EXP_BITS: Self::Bits = (Self::FREXP_BIAS_OFFSET << Self::MANTISSA_BITS) as _;
193        const MAX_LARGE_UINT: crate::LargeUInt = ((1 as crate::LargeUInt) << (Self::MANTISSA_BITS + 1)) as _;
194    }};
195
196    (COMMON) => {
197        #[inline(always)] fn try_from_int(value: crate::LargeInt) -> Option<Self> { FloatElementInternal::try_from_int(value) }
198        #[inline(always)] fn try_from_ratio(n: crate::LargeInt, d: crate::LargeInt) -> Option<Self> { FloatElementInternal::try_from_ratio(n, d) }
199
200        type ConstInt<const N: crate::LargeInt> = IntConst<N>;
201        type ConstRatio<const N: crate::LargeInt, const D: crate::LargeInt> = RatioConst<N, D>;
202
203        const HAS_INFINITY: bool = true;
204        const HAS_SIGNED_ZERO: bool = true;
205        const HAS_SUBNORMALS: bool = cfg!(not(feature = "ignore-denormals"));
206    };
207
208    (MUL_ADDE) => {
209        #[inline(always)] fn mul_adde(self, rhs: Self, acc: Self) -> Self { self * rhs + acc }
210        #[inline(always)] fn mul_sube(self, rhs: Self, acc: Self) -> Self { self * rhs - acc }
211        #[inline(always)] fn nmul_adde(self, rhs: Self, acc: Self) -> Self { acc - self * rhs }
212        #[inline(always)] fn nmul_sube(self, rhs: Self, acc: Self) -> Self { self * -rhs - acc }
213    };
214
215    ($t:ty $(: $f:ident)? => $bits:ty, $signed:ty { $($const:ident: $const_ty:ty = $value:expr;)* }) => {paste::paste! {
216        impl FloatElementInternal for $t {
217            #[inline(always)]
218            fn try_from_int(value: crate::LargeInt) -> Option<Self> {
219                if crate::likely(value.unsigned_abs() < Self::MAX_LARGE_UINT) {
220                    Some(value as $t) // safe to convert directly
221                } else {
222                    None
223                }
224            }
225
226            // This implementation is more accurate than simply doing n as f64 / d as f64,
227            // since that can lose precision when n and d are large but their ratio is small.
228            // Instead, we do integer division first, then add the fractional part. Although
229            // this is non-trivial, LLVM should optimize it down to a constant value when
230            // the inputs are known at compile time.
231            #[inline(always)]
232            fn try_from_ratio(n: crate::LargeInt, d: crate::LargeInt) -> Option<Self> {
233                if d == 0 {
234                    return None;
235                }
236
237                // fast path for values that both fit in the float exactly
238                if let (Some(n), Some(d)) = (<Self as FloatElementInternal>::try_from_int(n), <Self as FloatElementInternal>::try_from_int(d)) {
239                    return Some(n / d);
240                }
241
242                let (q, r) = (n / d, n % d);
243
244                let mut result = FloatElementInternal::try_from_int(q)?;
245
246                // d may not be exactly representable, but this will still scale it correctly
247                result += (r as $t) / (d as $t);
248
249                Some(result)
250            }
251        }
252
253        cfg_if::cfg_if! {
254            if #[cfg(feature = "std")] {
255                impl FloatElement for $t {
256                    #[inline(always)] fn sqrt(value: Self) -> Self { value.sqrt() }
257                    #[inline(always)] fn floor(value: Self) -> Self { value.floor() }
258                    #[inline(always)] fn ceil(value: Self) -> Self { value.ceil() }
259                    #[inline(always)] fn round(value: Self) -> Self { value.round() }
260                    #[inline(always)] fn trunc(value: Self) -> Self { value.trunc() }
261                    #[inline(always)] fn fract(value: Self) -> Self { value.fract() }
262                    #[inline(always)] fn next_up(value: Self) -> Self { value.next_up() }
263                    #[inline(always)] fn next_down(value: Self) -> Self { value.next_down() }
264
265                    impl_float_element!(COMMON);
266                }
267
268                impl MulAddExt for $t {
269                    type Output = Self;
270
271                    // trust the register implementation
272                    const HAS_TRUE_FMA: bool = <$t as FloatRegister>::HAS_TRUE_FMA;
273
274                    #[inline(always)] fn mul_add(self, rhs: Self, acc: Self) -> Self { <$t>::mul_add(self, rhs, acc) }
275                    #[inline(always)] fn mul_sub(self, rhs: Self, acc: Self) -> Self { <$t>::mul_add(self, rhs, -acc) }
276                    #[inline(always)] fn nmul_add(self, rhs: Self, acc: Self) -> Self { <$t>::mul_add(self, -rhs, acc) }
277                    #[inline(always)] fn nmul_sub(self, rhs: Self, acc: Self) -> Self { <$t>::mul_add(self, -rhs, -acc) }
278
279                    #[inline(always)] fn mul_adde(self, rhs: Self, acc: Self) -> Self { if !<Self as MulAddExt>::HAS_TRUE_FMA { self * rhs + acc } else { <$t>::mul_add(self, rhs, acc) } }
280                    #[inline(always)] fn mul_sube(self, rhs: Self, acc: Self) -> Self { if !<Self as MulAddExt>::HAS_TRUE_FMA { self * rhs - acc } else { <$t>::mul_add(self, rhs, -acc) } }
281                    #[inline(always)] fn nmul_adde(self, rhs: Self, acc: Self) -> Self { if !<Self as MulAddExt>::HAS_TRUE_FMA { acc - self * rhs } else { <$t>::mul_add(self, -rhs, acc) } }
282                    #[inline(always)] fn nmul_sube(self, rhs: Self, acc: Self) -> Self { if !<Self as MulAddExt>::HAS_TRUE_FMA { self * -rhs - acc } else { <$t>::mul_add(self, -rhs, -acc) } }
283                }
284            } else if #[cfg(all(feature = "spirv", target_arch = "spirv"))] {
285                impl FloatElement for $t {
286                    #[inline(always)] fn sqrt(value: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op1::<Self, Self, {crate::backend::spirv::arch::glsl::SQRT}, false>(value) } }
287                    #[inline(always)] fn floor(value: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op1::<Self, Self, {crate::backend::spirv::arch::glsl::FLOOR}, false>(value) } }
288                    #[inline(always)] fn ceil(value: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op1::<Self, Self, {crate::backend::spirv::arch::glsl::CEIL}, false>(value) } }
289                    #[inline(always)] fn round(value: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op1::<Self, Self, {crate::backend::spirv::arch::glsl::ROUND}, false>(value) } }
290                    #[inline(always)] fn trunc(value: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op1::<Self, Self, {crate::backend::spirv::arch::glsl::TRUNC}, false>(value) } }
291                    #[inline(always)] fn fract(value: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op1::<Self, Self, {crate::backend::spirv::arch::glsl::FRACT}, false>(value) } }
292                    #[inline(always)] fn next_up(value: Self) -> Self { value.next_up() }
293                    #[inline(always)] fn next_down(value: Self) -> Self { value.next_down() }
294
295                    impl_float_element!(COMMON);
296                }
297
298                impl MulAddExt for $t {
299                    type Output = Self;
300
301                    // GPU hardware always has FMA
302                    const HAS_TRUE_FMA: bool = true;
303
304                    #[inline(always)] fn mul_add(self, rhs: Self, acc: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op3::<Self, Self, Self, Self, {crate::backend::spirv::arch::glsl::FMA}, false>(self, rhs, acc) } }
305                    #[inline(always)] fn mul_sub(self, rhs: Self, acc: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op3::<Self, Self, Self, Self, {crate::backend::spirv::arch::glsl::FMA}, false>(self, rhs, -acc) } }
306                    #[inline(always)] fn nmul_add(self, rhs: Self, acc: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op3::<Self, Self, Self, Self, {crate::backend::spirv::arch::glsl::FMA}, false>(self, -rhs, acc) } }
307                    #[inline(always)] fn nmul_sub(self, rhs: Self, acc: Self) -> Self { unsafe { crate::backend::spirv::arch::glsl_op3::<Self, Self, Self, Self, {crate::backend::spirv::arch::glsl::FMA}, false>(self, -rhs, -acc) } }
308
309                    #[inline(always)] fn mul_adde(self, rhs: Self, acc: Self) -> Self { self.mul_add(rhs, acc) }
310                    #[inline(always)] fn mul_sube(self, rhs: Self, acc: Self) -> Self { self.mul_sub(rhs, acc) }
311                    #[inline(always)] fn nmul_adde(self, rhs: Self, acc: Self) -> Self { self.nmul_add(rhs, acc) }
312                    #[inline(always)] fn nmul_sube(self, rhs: Self, acc: Self) -> Self { self.nmul_sub(rhs, acc) }
313                }
314            } else if #[cfg(all(feature = "nightly", feature = "wasm", any(target_arch = "wasm32", target_arch = "wasm64")))] {
315                impl FloatElement for $t {
316                    // WASM has native scalar float ops for these
317                    #[inline(always)] fn sqrt(value: Self) -> Self { crate::backend::wasm::arch::[<$t _sqrt>](value) }
318                    #[inline(always)] fn floor(value: Self) -> Self { crate::backend::wasm::arch::[<$t _floor>](value) }
319                    #[inline(always)] fn ceil(value: Self) -> Self { crate::backend::wasm::arch::[<$t _ceil>](value) }
320                    #[inline(always)] fn trunc(value: Self) -> Self { crate::backend::wasm::arch::[<$t _trunc>](value) }
321                    #[inline(always)] fn fract(value: Self) -> Self { value - crate::backend::wasm::arch::[<$t _trunc>](value) }
322                    // WASM nearest() is banker's rounding (ties-to-even), not half-away-from-zero
323                    #[inline(always)] fn round(value: Self) -> Self { libm::[<round $($f)?>](value) }
324                    // No WASM scalar nextafter; fall back to libm
325                    #[inline(always)] fn next_up(value: Self) -> Self { libm::[<nextafter $($f)?>](value, Self::INFINITY) }
326                    #[inline(always)] fn next_down(value: Self) -> Self { libm::[<nextafter $($f)?>](value, Self::NEG_INFINITY) }
327
328                    impl_float_element!(COMMON);
329                }
330
331                impl MulAddExt for $t {
332                    type Output = Self;
333
334                    // No hardware scalar FMA on WASM; use libm for exact, separate ops for estimating
335                    const HAS_TRUE_FMA: bool = false;
336
337                    #[inline(always)] fn mul_add(self, rhs: Self, acc: Self) -> Self { libm::[<fma $($f)?>](self, rhs, acc) }
338                    #[inline(always)] fn mul_sub(self, rhs: Self, acc: Self) -> Self { libm::[<fma $($f)?>](self, rhs, -acc) }
339                    #[inline(always)] fn nmul_add(self, rhs: Self, acc: Self) -> Self { libm::[<fma $($f)?>](self, -rhs, acc) }
340                    #[inline(always)] fn nmul_sub(self, rhs: Self, acc: Self) -> Self { libm::[<fma $($f)?>](self, -rhs, -acc) }
341
342                    impl_float_element!(MUL_ADDE);
343                }
344            } else {
345                impl FloatElement for $t {
346                    #[inline(always)] fn sqrt(value: Self) -> Self { libm::[<sqrt $($f)?>](value) }
347                    #[inline(always)] fn floor(value: Self) -> Self { libm::[<floor $($f)?>](value) }
348                    #[inline(always)] fn ceil(value: Self) -> Self { libm::[<ceil $($f)?>](value) }
349                    #[inline(always)] fn round(value: Self) -> Self { libm::[<round $($f)?>](value) }
350                    #[inline(always)] fn trunc(value: Self) -> Self { libm::[<trunc $($f)?>](value) }
351                    #[inline(always)] fn next_up(value: Self) -> Self { libm::[<nextafter $($f)?>](value, Self::INFINITY) }
352                    #[inline(always)] fn next_down(value: Self) -> Self { libm::[<nextafter $($f)?>](value, Self::NEG_INFINITY) }
353
354                    impl_float_element!(COMMON);
355                }
356
357                impl MulAddExt for $t {
358                    type Output = Self;
359
360                    const HAS_TRUE_FMA: bool = false;
361
362                    #[inline(always)] fn mul_add(self, rhs: Self, acc: Self) -> Self { libm::[<fma $($f)?>](self, rhs, acc) }
363                    #[inline(always)] fn mul_sub(self, rhs: Self, acc: Self) -> Self { libm::[<fma $($f)?>](self, rhs, -acc) }
364                    #[inline(always)] fn nmul_add(self, rhs: Self, acc: Self) -> Self { libm::[<fma $($f)?>](self, -rhs, acc) }
365                    #[inline(always)] fn nmul_sub(self, rhs: Self, acc: Self) -> Self { libm::[<fma $($f)?>](self, -rhs, -acc) }
366
367                    impl_float_element!(MUL_ADDE);
368                }
369            }
370        }
371
372        impl FloatElementWithBits for $t {
373            type Bits = $bits;
374            type SignedBits = $signed;
375
376            impl_float_element!(CONSTS $($const: $const_ty = $value;)*);
377
378            #[inline(always)] fn from_f64(value: f64) -> Self { value as $t }
379            #[inline(always)] fn from_signed(value: Self::SignedBits) -> Self { value as $t }
380        }
381    }};
382}
383
384impl_float_element!(f32: f => u32, i32 {
385    EXP_BITS: u32 = 8;
386    MANTISSA_BITS: u32 = 23;
387    EXP_BIAS: i32 = 127;
388    MAX_BIASED_EXP: i32 = 255;
389
390    // 8 bits of exponent
391    EXP_LSB_MASK: u32 = 0xFF;
392
393    // Clear bits 23-30
394    SIGN_MANTISSA_MASK: u32 = 0x807F_FFFF;
395
396    // Canonical Quiet NaN: Sign=0, Exp=All 1s, Mantissa=100...0
397    // (Note: IEEE 754 allows many NaN patterns; this is just the standard "default")
398    NAN_PATTERN: Option<u32> = None; //0x7FC0_0000, but we don't want to use it
399
400    // Max Finite: Sign=0, Exp=254 (0xFE), Mantissa=All 1s
401    MAX_FINITE_PATTERN: u32 = 0x7F7F_FFFF;
402
403    MAX_SUBNORMAL: u32 = 0x007F_FFFF;
404    DENORMAL_TRICK: u32 = 0x0C800001;
405
406    // IMPLICIT_LEAD_BIT: bool = true;
407});
408
409impl_float_element!(f64 => u64, i64 {
410    EXP_BITS: u32 = 11;
411    MANTISSA_BITS: u32 = if cfg!(not(all(feature = "spirv", target_arch = "spirv", not(target_feature = "ext:Float64")))) { 52 } else { 0 };
412    EXP_BIAS: i64 = 1023;
413    MAX_BIASED_EXP: i64 = 2047;
414
415    // 11 bits of exponent
416    EXP_LSB_MASK: u64 = 0x7FF;
417
418    // Clear bits 52-62
419    SIGN_MANTISSA_MASK: u64 = 0x800F_FFFF_FFFF_FFFF;
420
421    // Canonical Quiet NaN: Sign=0, Exp=All 1s, Mantissa=100...0
422    NAN_PATTERN: Option<u64> = None; //0x7FF8_0000_0000_0000, but we don't want to use it
423
424    // Max Finite: Sign=0, Exp=2046 (0x7FE), Mantissa=All 1s
425    MAX_FINITE_PATTERN: u64 = 0x7FEF_FFFF_FFFF_FFFF;
426
427    MAX_SUBNORMAL: u64 = 0x000F_FFFF_FFFF_FFFF;
428    DENORMAL_TRICK: u64 = 0x0360000000000001;
429
430    // IMPLICIT_LEAD_BIT: bool = true;
431});
432
433#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
434pub enum RoundingMode {
435    NearestTiesToEven,
436    Truncate,
437
438    #[cfg(feature = "rand")]
439    Stochastic(u64),
440}