feat: add display UI with splash screen, dashboard, and RSSI sparkline

Implement full SSD1306 OLED display support: - Splash screen on boot (1.5s) with project name, version, board - Dashboard showing radio state, config, counters, RSSI/SNR - RSSI history sparkline (64-sample ring buffer, 2px-wide bars) - DisplayOn/DisplayOff host commands routed via dedicated channel - Pure render functions against generic DrawTarget (no hardware coupling) - Display on/off control bypasses radio task entirely Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
2026-03-29 20:42:54 -06:00
1 parent 7064110293
commit 3102fc747a
@@ -0,0 +1,279 @@
+//! `f32` extension
+
+use crate::float::F32;
+
+/// `f32` extension providing various arithmetic approximations and polyfills
+/// for `std` functionality.
+pub trait F32Ext: Sized {
+    /// Compute absolute value with a constant-time, data-independent
+    /// implementation.
+    fn abs(self) -> f32;
+
+    /// Approximates `acos(x)` in radians in the range `[0, pi]`
+    fn acos(self) -> f32;
+
+    /// Approximates `asin(x)` in radians in the range `[-pi/2, pi/2]`.
+    fn asin(self) -> f32;
+
+    /// Approximates `atan(x)` in radians with a maximum error of `0.002`.
+    fn atan(self) -> f32;
+
+    /// Approximates `atan(x)` normalized to the `[−1,1]` range with a maximum
+    /// error of `0.1620` degrees.
+    fn atan_norm(self) -> f32;
+
+    /// Approximates the four quadrant arctangent `atan2(x)` in radians, with
+    /// a maximum error of `0.002`.
+    fn atan2(self, other: f32) -> f32;
+
+    /// Approximates the four quadrant arctangent.
+    /// Normalized to the `[0,4)` range with a maximum error of `0.1620` degrees.
+    fn atan2_norm(self, other: f32) -> f32;
+
+    /// Approximates floating point ceiling.
+    fn ceil(self) -> f32;
+
+    /// Copies the sign from one number to another and returns it.
+    fn copysign(self, sign: f32) -> f32;
+
+    /// Approximates cosine in radians with a maximum error of `0.002`.
+    fn cos(self) -> f32;
+
+    /// Calculates Euclidean division, the matching method for `rem_euclid`.
+    fn div_euclid(self, other: f32) -> f32;
+
+    /// Approximates `e^x`.
+    fn exp(self) -> f32;
+
+    /// Approximates floating point floor.
+    fn floor(self) -> f32;
+
+    /// Retrieve the fractional part of floating point with sign.
+    fn fract(self) -> f32;
+
+    /// Approximates the length of the hypotenuse of a right-angle triangle given
+    /// legs of length `x` and `y`.
+    fn hypot(self, other: f32) -> f32;
+
+    /// Approximates `1/x` with an average deviation of ~8%.
+    fn inv(self) -> f32;
+
+    /// Approximates inverse square root with an average deviation of ~5%.
+    fn invsqrt(self) -> f32;
+
+    /// Approximates `ln(x)`.
+    fn ln(self) -> f32;
+
+    /// Approximates `log` with an arbitrary base.
+    fn log(self, base: f32) -> f32;
+
+    /// Approximates `log2`.
+    fn log2(self) -> f32;
+
+    /// Approximates `log10`.
+    fn log10(self) -> f32;
+
+    /// Computes `(self * a) + b`.
+    fn mul_add(self, a: f32, b: f32) -> f32;
+
+    /// Approximates `self^n`.
+    fn powf(self, n: f32) -> f32;
+
+    /// Approximates `self^n` where n is an `i32`
+    fn powi(self, n: i32) -> f32;
+
+    /// Returns the reciprocal (inverse) of a number, `1/x`.
+    fn recip(self) -> f32;
+
+    /// Calculates the least nonnegative remainder of `self (mod other)`.
+    fn rem_euclid(self, other: f32) -> f32;
+
+    /// Round the number part of floating point with sign.
+    fn round(self) -> f32;
+
+    /// Returns a number that represents the sign of `self`.
+    fn signum(self) -> f32;
+
+    /// Approximates sine in radians with a maximum error of `0.002`.
+    fn sin(self) -> f32;
+
+    /// Simultaneously computes the sine and cosine of the number, `x`.
+    /// Returns `(sin(x), cos(x))`.
+    fn sin_cos(self) -> (f32, f32);
+
+    /// Approximates square root with an average deviation of ~5%.
+    fn sqrt(self) -> f32;
+
+    /// Approximates `tan(x)` in radians with a maximum error of `0.6`.
+    fn tan(self) -> f32;
+
+    /// Retrieve whole number part of floating point with sign.
+    fn trunc(self) -> f32;
+}
+
+impl F32Ext for f32 {
+    #[inline]
+    fn abs(self) -> f32 {
+        F32(self).abs().0
+    }
+
+    #[inline]
+    fn acos(self) -> f32 {
+        F32(self).acos().0
+    }
+
+    #[inline]
+    fn asin(self) -> f32 {
+        F32(self).asin().0
+    }
+
+    #[inline]
+    fn atan(self) -> f32 {
+        F32(self).atan().0
+    }
+
+    #[inline]
+    fn atan_norm(self) -> f32 {
+        F32(self).atan_norm().0
+    }
+
+    #[inline]
+    fn atan2(self, other: f32) -> f32 {
+        F32(self).atan2(F32(other)).0
+    }
+
+    #[inline]
+    fn atan2_norm(self, other: f32) -> f32 {
+        F32(self).atan2_norm(F32(other)).0
+    }
+
+    #[inline]
+    fn ceil(self) -> f32 {
+        F32(self).ceil().0
+    }
+
+    #[inline]
+    fn copysign(self, sign: f32) -> f32 {
+        F32(self).copysign(F32(sign)).0
+    }
+
+    #[inline]
+    fn cos(self) -> f32 {
+        F32(self).cos().0
+    }
+
+    #[inline]
+    fn div_euclid(self, other: f32) -> f32 {
+        F32(self).div_euclid(F32(other)).0
+    }
+
+    #[inline]
+    fn exp(self) -> f32 {
+        F32(self).exp().0
+    }
+
+    #[inline]
+    fn floor(self) -> f32 {
+        F32(self).floor().0
+    }
+
+    #[inline]
+    fn fract(self) -> f32 {
+        F32(self).fract().0
+    }
+
+    #[inline]
+    fn hypot(self, other: f32) -> f32 {
+        F32(self).hypot(other.into()).0
+    }
+
+    #[inline]
+    fn inv(self) -> f32 {
+        F32(self).inv().0
+    }
+
+    #[inline]
+    fn invsqrt(self) -> f32 {
+        F32(self).invsqrt().0
+    }
+
+    #[inline]
+    fn ln(self) -> f32 {
+        F32(self).ln().0
+    }
+
+    #[inline]
+    fn log(self, base: f32) -> f32 {
+        F32(self).log(F32(base)).0
+    }
+
+    #[inline]
+    fn log2(self) -> f32 {
+        F32(self).log2().0
+    }
+
+    #[inline]
+    fn log10(self) -> f32 {
+        F32(self).log10().0
+    }
+
+    #[inline]
+    fn mul_add(self, a: f32, b: f32) -> f32 {
+        F32(self).mul_add(F32(a), F32(b)).0
+    }
+
+    #[inline]
+    fn powf(self, n: f32) -> f32 {
+        F32(self).powf(F32(n)).0
+    }
+
+    #[inline]
+    fn powi(self, n: i32) -> f32 {
+        F32(self).powi(n).0
+    }
+
+    #[inline]
+    fn recip(self) -> f32 {
+        F32(self).recip().0
+    }
+
+    #[inline]
+    fn rem_euclid(self, other: f32) -> f32 {
+        F32(self).rem_euclid(F32(other)).0
+    }
+
+    #[inline]
+    fn round(self) -> f32 {
+        F32(self).round().0
+    }
+
+    #[inline]
+    fn signum(self) -> f32 {
+        F32(self).signum().0
+    }
+
+    #[inline]
+    fn sin(self) -> f32 {
+        F32(self).sin().0
+    }
+
+    #[inline]
+    fn sin_cos(self) -> (f32, f32) {
+        (F32(self).sin().0, F32(self).cos().0)
+    }
+
+    #[inline]
+    fn sqrt(self) -> f32 {
+        F32(self).sqrt().0
+    }
+
+    #[inline]
+    fn tan(self) -> f32 {
+        F32(self).tan().0
+    }
+
+    #[inline]
+    fn trunc(self) -> f32 {
+        F32(self).trunc().0
+    }
+}
@@ -0,0 +1,649 @@
+//! Floating point operations
+
+pub(crate) mod abs;
+pub(crate) mod acos;
+pub(crate) mod asin;
+pub(crate) mod atan;
+pub(crate) mod atan2;
+pub(crate) mod ceil;
+pub(crate) mod copysign;
+pub(crate) mod cos;
+pub(crate) mod div_euclid;
+pub(crate) mod exp;
+pub(crate) mod floor;
+pub(crate) mod fract;
+pub(crate) mod hypot;
+pub(crate) mod inv;
+pub(crate) mod invsqrt;
+pub(crate) mod ln;
+pub(crate) mod log;
+pub(crate) mod log10;
+pub(crate) mod log2;
+pub(crate) mod mul_add;
+pub(crate) mod powf;
+pub(crate) mod powi;
+pub(crate) mod recip;
+pub(crate) mod rem_euclid;
+pub(crate) mod round;
+pub(crate) mod signum;
+pub(crate) mod sin;
+pub(crate) mod sin_cos;
+pub(crate) mod sqrt;
+pub(crate) mod tan;
+pub(crate) mod trunc;
+
+use core::{
+    cmp::Ordering,
+    fmt::{self, Display, LowerExp, UpperExp},
+    iter::{Product, Sum},
+    num::ParseFloatError,
+    ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign},
+    str::FromStr,
+};
+
+#[cfg(feature = "num-traits")]
+use num_traits::{Inv, Num, One, Zero};
+
+/// Sign mask.
+pub(crate) const SIGN_MASK: u32 = 0b1000_0000_0000_0000_0000_0000_0000_0000;
+
+/// Exponent mask.
+pub(crate) const EXPONENT_MASK: u32 = 0b0111_1111_1000_0000_0000_0000_0000_0000;
+
+/// Mantissa mask.
+pub(crate) const MANTISSA_MASK: u32 = 0b0000_0000_0111_1111_1111_1111_1111_1111;
+
+/// Exponent mask.
+pub(crate) const EXPONENT_BIAS: u32 = 127;
+
+/// Mantissa bits.
+///
+/// Note: `MANTISSA_DIGITS` is available in `core::f32`, but the actual bits taken up are 24 - 1.
+pub(crate) const MANTISSA_BITS: u32 = 23;
+
+/// 32-bit floating point wrapper which implements fast approximation-based
+/// operations.
+#[derive(Copy, Clone, Debug, Default, PartialEq, PartialOrd)]
+pub struct F32(pub f32);
+
+impl F32 {
+    /// The value `0.0`.
+    pub const ZERO: Self = Self(0.0);
+
+    /// The value `1.0`.
+    pub const ONE: Self = Self(1.0);
+
+    /// The radix or base of the internal representation of `f32`.
+    pub const RADIX: u32 = f32::RADIX;
+
+    /// Number of significant digits in base 2.
+    pub const MANTISSA_DIGITS: u32 = f32::MANTISSA_DIGITS;
+
+    /// Approximate number of significant digits in base 10.
+    pub const DIGITS: u32 = f32::DIGITS;
+
+    /// [Machine epsilon] value for `f32`.
+    ///
+    /// This is the difference between `1.0` and the next larger representable number.
+    ///
+    /// [Machine epsilon]: https://en.wikipedia.org/wiki/Machine_epsilon
+    pub const EPSILON: Self = Self(f32::EPSILON);
+
+    /// Smallest finite `f32` value.
+    pub const MIN: Self = Self(f32::MIN);
+
+    /// Smallest positive normal `f32` value.
+    pub const MIN_POSITIVE: Self = Self(f32::MIN_POSITIVE);
+
+    /// Largest finite `f32` value.
+    pub const MAX: Self = Self(f32::MAX);
+
+    /// One greater than the minimum possible normal power of 2 exponent.
+    pub const MIN_EXP: i32 = f32::MIN_EXP;
+
+    /// Maximum possible power of 2 exponent.
+    pub const MAX_EXP: i32 = f32::MAX_EXP;
+
+    /// Minimum possible normal power of 10 exponent.
+    pub const MIN_10_EXP: i32 = f32::MIN_10_EXP;
+
+    /// Maximum possible power of 10 exponent.
+    pub const MAX_10_EXP: i32 = f32::MAX_10_EXP;
+
+    /// Not a Number (NaN).
+    pub const NAN: Self = Self(f32::NAN);
+
+    /// Infinity (∞).
+    pub const INFINITY: Self = Self(f32::INFINITY);
+
+    /// Negative infinity (−∞).
+    pub const NEG_INFINITY: Self = Self(f32::NEG_INFINITY);
+
+    /// Returns `true` if this value is `NaN`.
+    #[inline]
+    pub fn is_nan(self) -> bool {
+        self.0.is_nan()
+    }
+
+    /// Returns `true` if this value is positive infinity or negative infinity, and
+    /// `false` otherwise.
+    #[inline]
+    pub fn is_infinite(self) -> bool {
+        self.0.is_infinite()
+    }
+
+    /// Returns `true` if this number is neither infinite nor `NaN`.
+    #[inline]
+    pub fn is_finite(self) -> bool {
+        self.0.is_finite()
+    }
+
+    /// Returns `true` if `self` has a positive sign, including `+0.0`, `NaN`s with
+    /// positive sign bit and positive infinity.
+    #[inline]
+    pub fn is_sign_positive(self) -> bool {
+        self.0.is_sign_positive()
+    }
+
+    /// Returns `true` if `self` has a negative sign, including `-0.0`, `NaN`s with
+    /// negative sign bit and negative infinity.
+    #[inline]
+    pub fn is_sign_negative(self) -> bool {
+        self.0.is_sign_negative()
+    }
+
+    /// Raw transmutation to `u32`.
+    ///
+    /// This is currently identical to `transmute::<f32, u32>(self)` on all platforms.
+    ///
+    /// See [`F32::from_bits`] for some discussion of the portability of this operation
+    /// (there are almost no issues).
+    #[inline]
+    pub fn to_bits(self) -> u32 {
+        self.0.to_bits()
+    }
+
+    /// Raw transmutation from `u32`.
+    ///
+    /// This is currently identical to `transmute::<u32, f32>(v)` on all platforms.
+    /// It turns out this is incredibly portable, for two reasons:
+    ///
+    /// - Floats and Ints have the same endianness on all supported platforms.
+    /// - IEEE-754 very precisely specifies the bit layout of floats.
+    ///
+    /// See [`f32::from_bits`] for more information.
+    #[inline]
+    pub fn from_bits(v: u32) -> Self {
+        Self(f32::from_bits(v))
+    }
+
+    /// Extract exponent bits.
+    pub(crate) fn extract_exponent_bits(self) -> u32 {
+        (self.to_bits() & EXPONENT_MASK)
+            .overflowing_shr(MANTISSA_BITS)
+            .0
+    }
+
+    /// Extract the exponent of a float's value.
+    pub(crate) fn extract_exponent_value(self) -> i32 {
+        (self.extract_exponent_bits() as i32) - EXPONENT_BIAS as i32
+    }
+
+    /// Remove sign.
+    pub(crate) fn without_sign(self) -> Self {
+        Self::from_bits(self.to_bits() & !SIGN_MASK)
+    }
+
+    /// Set the exponent to the given value.
+    pub(crate) fn set_exponent(self, exponent: i32) -> Self {
+        debug_assert!(exponent <= 127 && exponent >= -128);
+        let without_exponent: u32 = self.to_bits() & !EXPONENT_MASK;
+        let only_exponent: u32 = ((exponent + EXPONENT_BIAS as i32) as u32)
+            .overflowing_shl(MANTISSA_BITS)
+            .0;
+
+        Self::from_bits(without_exponent | only_exponent)
+    }
+
+    /// Is this floating point value equivalent to an integer?
+    pub(crate) fn is_integer(&self) -> bool {
+        let exponent = self.extract_exponent_value();
+        let self_bits = self.to_bits();
+
+        // if exponent is negative we shouldn't remove anything, this stops an opposite shift.
+        let exponent_clamped = i32::max(exponent, 0) as u32;
+
+        // find the part of the fraction that would be left over
+        let fractional_part = (self_bits).overflowing_shl(exponent_clamped).0 & MANTISSA_MASK;
+
+        // if fractional part contains anything, we know it *isn't* an integer.
+        // if zero there will be nothing in the fractional part
+        // if it is whole, there will be nothing in the fractional part
+        fractional_part == 0
+    }
+
+    /// Is this floating point value even?
+    fn is_even(&self) -> bool {
+        // any floating point value that doesn't fit in an i32 range is even,
+        // and will loose 1's digit precision at exp values of 23+
+        if self.extract_exponent_value() >= 31 {
+            true
+        } else {
+            (self.0 as i32) % 2 == 0
+        }
+    }
+}
+
+impl Add for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn add(self, rhs: F32) -> F32 {
+        F32(self.0 + rhs.0)
+    }
+}
+
+impl Add<f32> for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn add(self, rhs: f32) -> F32 {
+        F32(self.0 + rhs)
+    }
+}
+
+impl Add<F32> for f32 {
+    type Output = F32;
+
+    #[inline]
+    fn add(self, rhs: F32) -> F32 {
+        F32(self + rhs.0)
+    }
+}
+
+impl AddAssign for F32 {
+    #[inline]
+    fn add_assign(&mut self, rhs: F32) {
+        self.0 += rhs.0;
+    }
+}
+
+impl AddAssign<f32> for F32 {
+    #[inline]
+    fn add_assign(&mut self, rhs: f32) {
+        self.0 += rhs;
+    }
+}
+
+impl AddAssign<F32> for f32 {
+    #[inline]
+    fn add_assign(&mut self, rhs: F32) {
+        *self += rhs.0;
+    }
+}
+
+impl Display for F32 {
+    #[inline]
+    fn fmt(&self, fmt: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(fmt, "{}", self.0)
+    }
+}
+
+impl Div for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn div(self, rhs: F32) -> F32 {
+        F32(self.0 / rhs.0)
+    }
+}
+
+impl Div<f32> for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn div(self, rhs: f32) -> F32 {
+        F32(self.0 / rhs)
+    }
+}
+
+impl Div<F32> for f32 {
+    type Output = F32;
+
+    #[inline]
+    fn div(self, rhs: F32) -> F32 {
+        F32(self / rhs.0)
+    }
+}
+
+impl DivAssign for F32 {
+    #[inline]
+    fn div_assign(&mut self, rhs: F32) {
+        self.0 /= rhs.0;
+    }
+}
+
+impl DivAssign<f32> for F32 {
+    #[inline]
+    fn div_assign(&mut self, rhs: f32) {
+        self.0 /= rhs;
+    }
+}
+
+impl DivAssign<F32> for f32 {
+    #[inline]
+    fn div_assign(&mut self, rhs: F32) {
+        *self /= rhs.0;
+    }
+}
+
+impl From<f32> for F32 {
+    #[inline]
+    fn from(n: f32) -> F32 {
+        F32(n)
+    }
+}
+
+impl From<F32> for f32 {
+    #[inline]
+    fn from(n: F32) -> f32 {
+        n.0
+    }
+}
+
+impl From<i8> for F32 {
+    #[inline]
+    fn from(n: i8) -> F32 {
+        F32(n.into())
+    }
+}
+
+impl From<i16> for F32 {
+    #[inline]
+    fn from(n: i16) -> F32 {
+        F32(n.into())
+    }
+}
+
+impl From<u8> for F32 {
+    #[inline]
+    fn from(n: u8) -> F32 {
+        F32(n.into())
+    }
+}
+
+impl From<u16> for F32 {
+    #[inline]
+    fn from(n: u16) -> F32 {
+        F32(n.into())
+    }
+}
+
+impl FromStr for F32 {
+    type Err = ParseFloatError;
+
+    #[inline]
+    fn from_str(src: &str) -> Result<F32, ParseFloatError> {
+        f32::from_str(src).map(F32)
+    }
+}
+
+impl LowerExp for F32 {
+    #[inline]
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "{:e}", self.0)
+    }
+}
+
+impl Mul for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn mul(self, rhs: F32) -> F32 {
+        F32(self.0 * rhs.0)
+    }
+}
+
+impl Mul<f32> for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn mul(self, rhs: f32) -> F32 {
+        F32(self.0 * rhs)
+    }
+}
+
+impl Mul<F32> for f32 {
+    type Output = F32;
+
+    #[inline]
+    fn mul(self, rhs: F32) -> F32 {
+        F32(self * rhs.0)
+    }
+}
+
+impl MulAssign for F32 {
+    #[inline]
+    fn mul_assign(&mut self, rhs: F32) {
+        self.0 *= rhs.0;
+    }
+}
+
+impl MulAssign<f32> for F32 {
+    #[inline]
+    fn mul_assign(&mut self, rhs: f32) {
+        self.0 *= rhs;
+    }
+}
+
+impl MulAssign<F32> for f32 {
+    #[inline]
+    fn mul_assign(&mut self, rhs: F32) {
+        *self *= rhs.0;
+    }
+}
+
+impl Neg for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn neg(self) -> F32 {
+        F32(-self.0)
+    }
+}
+
+impl PartialEq<f32> for F32 {
+    fn eq(&self, other: &f32) -> bool {
+        self.0.eq(other)
+    }
+}
+
+impl PartialEq<F32> for f32 {
+    fn eq(&self, other: &F32) -> bool {
+        self.eq(&other.0)
+    }
+}
+
+impl PartialOrd<f32> for F32 {
+    fn partial_cmp(&self, other: &f32) -> Option<Ordering> {
+        self.0.partial_cmp(other)
+    }
+}
+
+impl PartialOrd<F32> for f32 {
+    fn partial_cmp(&self, other: &F32) -> Option<Ordering> {
+        self.partial_cmp(&other.0)
+    }
+}
+
+impl Product for F32 {
+    #[inline]
+    fn product<I>(iter: I) -> Self
+    where
+        I: Iterator<Item = F32>,
+    {
+        F32(f32::product(iter.map(f32::from)))
+    }
+}
+
+impl Rem for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn rem(self, rhs: F32) -> F32 {
+        F32(self.0 % rhs.0)
+    }
+}
+
+impl Rem<f32> for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn rem(self, rhs: f32) -> F32 {
+        F32(self.0 % rhs)
+    }
+}
+
+impl Rem<F32> for f32 {
+    type Output = F32;
+
+    #[inline]
+    fn rem(self, rhs: F32) -> F32 {
+        F32(self % rhs.0)
+    }
+}
+
+impl RemAssign for F32 {
+    #[inline]
+    fn rem_assign(&mut self, rhs: F32) {
+        self.0 %= rhs.0;
+    }
+}
+
+impl RemAssign<f32> for F32 {
+    #[inline]
+    fn rem_assign(&mut self, rhs: f32) {
+        self.0 %= rhs;
+    }
+}
+
+impl Sub for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn sub(self, rhs: F32) -> F32 {
+        F32(self.0 - rhs.0)
+    }
+}
+
+impl Sub<f32> for F32 {
+    type Output = F32;
+
+    #[inline]
+    fn sub(self, rhs: f32) -> F32 {
+        F32(self.0 - rhs)
+    }
+}
+
+impl Sub<F32> for f32 {
+    type Output = F32;
+
+    #[inline]
+    fn sub(self, rhs: F32) -> F32 {
+        F32(self - rhs.0)
+    }
+}
+
+impl SubAssign for F32 {
+    #[inline]
+    fn sub_assign(&mut self, rhs: F32) {
+        self.0 -= rhs.0;
+    }
+}
+
+impl SubAssign<f32> for F32 {
+    #[inline]
+    fn sub_assign(&mut self, rhs: f32) {
+        self.0 -= rhs;
+    }
+}
+
+impl SubAssign<F32> for f32 {
+    #[inline]
+    fn sub_assign(&mut self, rhs: F32) {
+        *self -= rhs.0;
+    }
+}
+
+impl Sum for F32 {
+    #[inline]
+    fn sum<I>(iter: I) -> Self
+    where
+        I: Iterator<Item = F32>,
+    {
+        F32(f32::sum(iter.map(f32::from)))
+    }
+}
+
+impl UpperExp for F32 {
+    #[inline]
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "{:E}", self.0)
+    }
+}
+
+#[cfg(feature = "num-traits")]
+#[cfg_attr(docsrs, doc(cfg(feature = "num-traits")))]
+impl Zero for F32 {
+    fn zero() -> Self {
+        Self::ZERO
+    }
+
+    fn is_zero(&self) -> bool {
+        Self::ZERO == *self
+    }
+}
+
+#[cfg(feature = "num-traits")]
+#[cfg_attr(docsrs, doc(cfg(feature = "num-traits")))]
+impl One for F32 {
+    fn one() -> Self {
+        Self::ONE
+    }
+
+    fn is_one(&self) -> bool {
+        Self::ONE == *self
+    }
+}
+
+#[cfg(feature = "num-traits")]
+#[cfg_attr(docsrs, doc(cfg(feature = "num-traits")))]
+impl Num for F32 {
+    type FromStrRadixErr = num_traits::ParseFloatError;
+
+    fn from_str_radix(str: &str, radix: u32) -> Result<Self, Self::FromStrRadixErr> {
+        f32::from_str_radix(str, radix).map(Self)
+    }
+}
+
+#[cfg(feature = "num-traits")]
+#[cfg_attr(docsrs, doc(cfg(feature = "num-traits")))]
+impl Inv for F32 {
+    type Output = Self;
+
+    fn inv(self) -> Self {
+        self.inv()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    #[allow(unused_imports)] // remove when we have more tests
+    use super::F32;
+
+    #[cfg(feature = "num-traits")]
+    #[test]
+    fn inv_trait() {
+        assert_eq!(num_traits::Inv::inv(F32(2.0)), F32(0.5));
+    }
+}
@@ -0,0 +1,31 @@
+//! Compute the absolute value of a single-precision float.
+//!
+//! Method described at: <https://bits.stephan-brumme.com/absFloat.html>
+
+use super::{F32, SIGN_MASK};
+
+impl F32 {
+    /// Computes the absolute value of `self`.
+    ///
+    /// Returns [`Self::NAN`] if the number is [`Self::NAN`].
+    pub fn abs(self) -> Self {
+        Self::from_bits(self.to_bits() & !SIGN_MASK)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32::ONE.abs(), 1.0);
+        assert_eq!(F32::ZERO.abs(), 0.0);
+        assert_eq!(F32(-1.0).abs(), 1.0);
+    }
+
+    #[test]
+    fn nan() {
+        assert!(F32::NAN.abs().is_nan());
+    }
+}
@@ -0,0 +1,70 @@
+//! arccos approximation for a single-precision float.
+//!
+//! Method described at:
+//! <https://math.stackexchange.com/questions/2908908/express-arccos-in-terms-of-arctan>
+
+use super::F32;
+use core::f32::consts::PI;
+
+impl F32 {
+    /// Computes `acos(x)` approximation in radians in the range `[0, pi]`.
+    pub(crate) fn acos(self) -> Self {
+        if self > 0.0 {
+            ((Self::ONE - self * self).sqrt() / self).atan()
+        } else if self == 0.0 {
+            Self(PI / 2.0)
+        } else {
+            ((Self::ONE - self * self).sqrt() / self).atan() + PI
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+    use core::f32::consts;
+
+    const MAX_ERROR: f32 = 0.03;
+
+    #[test]
+    fn sanity_check() {
+        // Arccosine test vectors - `(input, output)`
+        let test_vectors: &[(f32, f32)] = &[
+            (2.000, f32::NAN),
+            (1.000, 0.0),
+            (0.866, consts::FRAC_PI_6),
+            (0.707, consts::FRAC_PI_4),
+            (0.500, consts::FRAC_PI_3),
+            (f32::EPSILON, consts::FRAC_PI_2),
+            (0.000, consts::FRAC_PI_2),
+            (-f32::EPSILON, consts::FRAC_PI_2),
+            (-0.500, 2.0 * consts::FRAC_PI_3),
+            (-0.707, 3.0 * consts::FRAC_PI_4),
+            (-0.866, 5.0 * consts::FRAC_PI_6),
+            (-1.000, consts::PI),
+            (-2.000, f32::NAN),
+        ];
+
+        for &(x, expected) in test_vectors {
+            let actual = F32(x).acos();
+            if expected.is_nan() {
+                assert!(
+                    actual.is_nan(),
+                    "acos({}) returned {}, should be NAN",
+                    x,
+                    actual
+                );
+            } else {
+                let delta = (actual - expected).abs();
+
+                assert!(
+                    delta <= MAX_ERROR,
+                    "delta {} too large: {} vs {}",
+                    delta,
+                    actual,
+                    expected
+                );
+            }
+        }
+    }
+}
@@ -0,0 +1,25 @@
+//! arcsin approximation for a single-precision float.
+//!
+//! Method described at:
+//! <https://dsp.stackexchange.com/questions/25770/looking-for-an-arcsin-algorithm>
+
+use super::F32;
+
+impl F32 {
+    /// Computes `asin(x)` approximation in radians in the range `[-pi/2, pi/2]`.
+    pub fn asin(self) -> Self {
+        (self * (Self::ONE - self * self).invsqrt()).atan()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+    use core::f32::consts::FRAC_PI_2;
+
+    #[test]
+    fn sanity_check() {
+        let difference = F32(FRAC_PI_2).sin().asin() - FRAC_PI_2;
+        assert!(difference.abs() <= F32::EPSILON);
+    }
+}
@@ -0,0 +1,80 @@
+//! arctangent approximation for a single-precision float.
+//!
+//! Method described at:
+//! <https://ieeexplore.ieee.org/document/6375931>
+
+use super::F32;
+use core::f32::consts::FRAC_PI_2;
+
+impl F32 {
+    /// Approximates `atan(x)` approximation in radians with a maximum error of
+    /// `0.002`.
+    ///
+    /// Returns [`Self::NAN`] if the number is [`Self::NAN`].
+    pub fn atan(self) -> Self {
+        FRAC_PI_2 * self.atan_norm()
+    }
+
+    /// Approximates `atan(x)` normalized to the `[−1,1]` range with a maximum
+    /// error of `0.1620` degrees.
+    pub fn atan_norm(self) -> Self {
+        const SIGN_MASK: u32 = 0x8000_0000;
+        const B: f32 = 0.596_227;
+
+        // Extract the sign bit
+        let ux_s = SIGN_MASK & self.to_bits();
+
+        // Calculate the arctangent in the first quadrant
+        let bx_a = (B * self).abs();
+        let n = bx_a + self * self;
+        let atan_1q = n / (1.0 + bx_a + n);
+
+        // Restore the sign bit and convert to float
+        Self::from_bits(ux_s | atan_1q.to_bits())
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+    use core::f32::consts;
+
+    /// 0.1620 degrees in radians
+    const MAX_ERROR: f32 = 0.003;
+
+    #[test]
+    fn sanity_check() {
+        // Arctangent test vectors - `(input, output)`
+        let test_vectors: &[(f32, f32)] = &[
+            (3.0_f32.sqrt() / 3.0, consts::FRAC_PI_6),
+            (1.0, consts::FRAC_PI_4),
+            (3.0_f32.sqrt(), consts::FRAC_PI_3),
+            (-(3.0_f32.sqrt()) / 3.0, -consts::FRAC_PI_6),
+            (-1.0, -consts::FRAC_PI_4),
+            (-(3.0_f32.sqrt()), -consts::FRAC_PI_3),
+        ];
+
+        for &(x, expected) in test_vectors {
+            let actual = F32(x).atan().0;
+            let delta = actual - expected;
+
+            assert!(
+                delta <= MAX_ERROR,
+                "delta {} too large: {} vs {}",
+                delta,
+                actual,
+                expected
+            );
+        }
+    }
+
+    #[test]
+    fn zero() {
+        assert_eq!(F32::ZERO.atan(), F32::ZERO);
+    }
+
+    #[test]
+    fn nan() {
+        assert!(F32::NAN.atan().is_nan());
+    }
+}
@@ -0,0 +1,79 @@
+//! Four quadrant arctangent approximation for a single-precision float.
+//!
+//! Method described at: <https://ieeexplore.ieee.org/document/6375931>
+
+use super::F32;
+use core::f32::consts::PI;
+
+impl F32 {
+    /// Approximates the four quadrant arctangent of `self` (`y`) and
+    /// `rhs` (`x`) in radians with a maximum error of `0.002`.
+    ///
+    /// - `x = 0`, `y = 0`: `0`
+    /// - `x >= 0`: `arctan(y/x)` -> `[-pi/2, pi/2]`
+    /// - `y >= 0`: `arctan(y/x) + pi` -> `(pi/2, pi]`
+    /// - `y < 0`: `arctan(y/x) - pi` -> `(-pi, -pi/2)`
+    pub fn atan2(self, rhs: Self) -> Self {
+        let n = self.atan2_norm(rhs);
+        PI / 2.0 * if n > 2.0 { n - 4.0 } else { n }
+    }
+
+    /// Approximates `atan2(y,x)` normalized to the `[0, 4)` range with a maximum
+    /// error of `0.1620` degrees.
+    pub(crate) fn atan2_norm(self, rhs: Self) -> Self {
+        const SIGN_MASK: u32 = 0x8000_0000;
+        const B: f32 = 0.596_227;
+
+        let y = self;
+        let x = rhs;
+
+        // Extract sign bits from floating point values
+        let ux_s = SIGN_MASK & x.to_bits();
+        let uy_s = SIGN_MASK & y.to_bits();
+
+        // Determine quadrant offset
+        let q = ((!ux_s & uy_s) >> 29 | ux_s >> 30) as f32;
+
+        // Calculate arctangent in the first quadrant
+        let bxy_a = (B * x * y).abs();
+        let n = bxy_a + y * y;
+        let atan_1q = n / (x * x + bxy_a + n);
+
+        // Translate it to the proper quadrant
+        let uatan_2q = (ux_s ^ uy_s) | atan_1q.to_bits();
+        Self(q) + Self::from_bits(uatan_2q)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+    use core::f32::consts::PI;
+
+    /// 0.1620 degrees in radians
+    const MAX_ERROR: f32 = 0.003;
+
+    #[test]
+    fn sanity_check() {
+        let test_vectors: &[(f32, f32, f32)] = &[
+            (0.0, 1.0, 0.0),
+            (0.0, -1.0, PI),
+            (3.0, 2.0, (3.0f32 / 2.0).atan()),
+            (2.0, -1.0, (2.0f32 / -1.0).atan() + PI),
+            (-2.0, -1.0, (-2.0f32 / -1.0).atan() - PI),
+        ];
+
+        for &(y, x, expected) in test_vectors {
+            let actual = F32(y).atan2(F32(x)).0;
+            let delta = actual - expected;
+
+            assert!(
+                delta <= MAX_ERROR,
+                "delta {} too large: {} vs {}",
+                delta,
+                actual,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,25 @@
+//! Floating point ceiling approximation for a single-precision float.
+
+use super::F32;
+
+impl F32 {
+    /// Returns the smallest integer greater than or equal to a number.
+    pub fn ceil(self) -> Self {
+        -(-self).floor()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32(-1.1).ceil().0, -1.0);
+        assert_eq!(F32(-0.1).ceil().0, 0.0);
+        assert_eq!(F32(0.0).ceil().0, 0.0);
+        assert_eq!(F32(1.0).ceil().0, 1.0);
+        assert_eq!(F32(1.1).ceil().0, 2.0);
+        assert_eq!(F32(2.9).ceil().0, 3.0);
+    }
+}
@@ -0,0 +1,33 @@
+//! Copy the sign over from another number.
+
+use super::{F32, SIGN_MASK};
+
+impl F32 {
+    /// Returns a number composed of the magnitude of `self` and the sign of
+    /// `sign`.
+    pub fn copysign(self, sign: Self) -> Self {
+        let source_bits = sign.to_bits();
+        let source_sign = source_bits & SIGN_MASK;
+        let signless_destination_bits = self.to_bits() & !SIGN_MASK;
+        Self::from_bits(signless_destination_bits | source_sign)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32(1.0).copysign(F32(-1.0)).0, -1.0);
+        assert_eq!(F32(-1.0).copysign(F32(1.0)).0, 1.0);
+        assert_eq!(F32(1.0).copysign(F32(1.0)).0, 1.0);
+        assert_eq!(F32(-1.0).copysign(F32(-1.0)).0, -1.0);
+
+        let large_float = F32(100_000_000.13425345345);
+        assert_eq!(large_float.copysign(-large_float), -large_float);
+        assert_eq!((-large_float).copysign(large_float), large_float);
+        assert_eq!(large_float.copysign(large_float), large_float);
+        assert_eq!((-large_float).copysign(-large_float), -large_float);
+    }
+}
@@ -0,0 +1,308 @@
+//! Cosine function implemented using a parabolic approximation.
+//!
+//! Note that while this method is inspired by a [Taylor series] approximation,
+//! it provides both better performance and accuracy.
+//!
+//! Most of the other trigonometric functions in this crate are implemented in
+//! terms of this approximation.
+//!
+//! Method described: <https://stackoverflow.com/posts/28050328/revisions>
+//!
+//! Originally based on this approximation by Nick Capen:
+//!
+//! <https://web.archive.org/web/20180105200343if_/http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648>
+//!
+//! Post quoted below for posterity as the original site is down:
+//!
+//! > In some situations you need a good approximation of sine and cosine that
+//! > runs at very high performance. One example is to implement dynamic
+//! > subdivision of round surfaces, comparable to those in Quake 3. Or to
+//! > implement wave motion, in case there are no vertex shaders 2.0 available.
+//! >
+//! > The standard C `sinf()` and `cosf()` functions are terribly slow and offer
+//! > much more precision than what we actually need. What we really want is an
+//! > approximation that offers the best compromise between precision and
+//! > performance. The most well known approximation method is to use a
+//! > [Taylor series] about 0 (also known as a Maclaurin series), which for
+//! > sine becomes:
+//! >
+//! > ```text
+//! > x - 1/6 x^3 + 1/120 x^5 - 1/5040 x^7 + ...
+//! > ```
+//! >
+//! > When we plot this we get: taylor.gif. The green line is the real sine,
+//! > the red line is the first four terms of the Taylor series. This seems
+//! > like an acceptable approximation, but lets have a closer look:
+//! > taylor_zoom.gif. It behaves very nicely up till `pi/2`, but after that it
+//! > quickly deviates. At pi, it evaluates to -0.075 instead of 0. Using this
+//! > for wave simulation will result in jerky motion which is unacceptable.
+//! >
+//! > We could add another term, which in fact reduces the error significantly,
+//! > but this makes the formula pretty lengthy. We already need 7
+//! > multiplications and 3 additions for the 4-term version. The Taylor series
+//! > just can't give us the precision and performance we're looking for.
+//! >
+//! > We did however note that we need `sin(pi) = 0`. And there's another thing
+//! > we can see from taylor_zoom.gif: this looks very much like a parabola!
+//! > So let's try to find the formula of a parabola that matches it as closely
+//! > as possible. The generic formula for a parabola is `A + B x + C x²`. So
+//! > this gives us three degrees of freedom. Obvious choices are that we want
+//! > `sin(0) = 0`, `sin(pi/2) = 1` and `sin(pi) = 0`. This gives us the following
+//! > three equations:
+//! >
+//! > ```text
+//! > A + B 0 + C 0² = 0
+//! > A + B pi/2 + C (pi/2)² = 1
+//! > A + B pi + C pi² = 0
+//! > ```
+//! >
+//! > Which has as solution `A = 0, B = 4/pi, C = -4/pi²`. So our parabola
+//! > approximation becomes `4/pi x - 4/pi² x²`. Plotting this we get:
+//! > parabola.gif. This looks worse than the 4-term Taylor series, right?
+//! > Wrong! The maximum absolute error is 0.056. Furthermore, this
+//! > approximation will give us smooth wave motion, and can be calculated
+//! > in only 3 multiplications and 1 addition!
+//! >
+//! > Unfortunately it's not very practical yet. This is what we get in the
+//! > `[-pi, pi]` range: negative.gif. It's quite obvious we want at least a
+//! > full period. But it's also clear that it's just another parabola,
+//! > mirrored around the origin. The formula for it is `4/pi x + 4/pi² x²`.
+//! > So the straightforward (pseudo-C) solution is:
+//! >
+//! > ```text
+//! > if(x > 0)
+//! > {
+//! >     y = 4/pi x - 4/pi² x²;
+//! > }
+//! > else
+//! > {
+//! >     y = 4/pi x + 4/pi² x²;
+//! > }
+//! > ```
+//! >
+//! > Adding a branch is not a good idea though. It makes the code
+//! > significantly slower. But look at how similar the two parts really are.
+//! > A subtraction becomes an addition depending on the sign of x. In a naive
+//! > first attempt to eliminate the branch we can 'extract' the sign of x
+//! > using `x / abs(x)`. The division is very expensive but look at the resulting
+//! > formula: `4/pi x - x / abs(x) 4/pi² x²`. By inverting the division we can
+//! > simplify this to a very nice and clean `4/pi x - 4/pi² x abs(x)`. So for
+//! > just one extra operation we get both halves of our sine approximation!
+//! > Here's the graph of this formula confirming the result: abs.gif.
+//! >
+//! > Now let's look at cosine. Basic trigonometry tells us that
+//! > `cos(x) = sin(pi/2 + x)`. Is that all, add pi/2 to x? No, we actually get
+//! > the unwanted part of the parabola again: shift_sine.gif. What we need to
+//! > do is 'wrap around' when x > pi/2. This can be done by subtracting 2 pi.
+//! > So the code becomes:
+//! >
+//! > ```text
+//! > x += pi/2;
+//! >
+//! > if(x > pi)   // Original x > pi/2
+//! > {
+//! >     x -= 2 * pi;   // Wrap: cos(x) = cos(x - 2 pi)
+//! > }
+//! >
+//! > y = sin(x);
+//! > ```
+//! >
+//! > Yet another branch. To eliminate it, we can use binary logic, like this:
+//! >
+//! > ```text
+//! > x -= (x > pi) & (2 * pi);
+//! > ```
+//! >
+//! > Note that this isn't valid C code at all. But it should clarify how this
+//! > can work. When x > pi is false the & operation zeroes the right hand part
+//! > so the subtraction does nothing, which is perfectly equivalent. I'll
+//! > leave it as an excercise to the reader to create working C code for this
+//! > (or just read on). Clearly the cosine requires a few more operations than
+//! > the sine but there doesn't seem to be any other way and it's still
+//! > extremely fast.
+//! >
+//! > Now, a maximum error of 0.056 is nice but clearly the 4-term Taylor series
+//! > still has a smaller error on average. Recall what our sine looked like:
+//! > abs.gif. So isn't there anything we can do to further increase precision
+//! > at minimal cost? Of course the current version is already applicable for
+//! > many situations where something that looks like a sine is just as good as
+//! > the real sine. But for other situations that's just not good enough.
+//! >
+//! > Looking at the graphs you'll notice that our approximation always
+//! > overestimates the real sine, except at 0, pi/2 and pi. So what we need is
+//! > to 'scale it down' without touching these important points. The solution
+//! > is to use the squared parabola, which looks like this: squared.gif. Note
+//! > how it preserves those important points but it's always lower than the
+//! > real sine. So we can use a weighted average of the two to get a better
+//! > approximation:
+//! >
+//! > ```text
+//! > Q (4/pi x - 4/pi² x²) + P (4/pi x - 4/pi² x²)²
+//! > ```
+//! >
+//! > With `Q + P = 1`. You can use exact minimization of either the absolute
+//! > or relative error but I'll save you the math. The optimal weights are
+//! > `Q = 0.775`, `P = 0.225` for the absolute error and `Q = 0.782`,
+//! > `P = 0.218` for the relative error. I will use the former. The resulting
+//! > graph is: average.gif. Where did the red line go? It's almost entirely
+//! > covered by the green line, which instantly shows how good this
+//! > approximation really is. The maximum error is about 0.001, a 50x
+//! > improvement! The formula looks lengthy but the part between parenthesis
+//! > is the same value from the parabola, which needs to be computed only
+//! > once. In fact only 2 extra multiplications and 2 additions are required
+//! > to achieve this precision boost.
+//! >
+//! > It shouldn't come as a big surprise that to make it work for negative x
+//! > as well we need a second abs() operation. The final C code for the sine
+//! > becomes:
+//! >
+//! > ```text
+//! > float sin(float x)
+//! > {
+//! >     const float B = 4/pi;
+//! >     const float C = -4/(pi*pi);
+//! >
+//! >     float y = B * x + C * x * abs(x);
+//! >
+//! >     #ifdef EXTRA_PRECISION
+//! >     //  const float Q = 0.775;
+//! >         const float P = 0.225;
+//! >
+//! >         y = P * (y * abs(y) - y) + y;   // Q * y + P * y * abs(y)
+//! >     #endif
+//! > }
+//! > ```
+//! >
+//! > So we need merely 5 multiplications and 3 additions; still faster than
+//! > the 4-term Taylor if we neglect the abs(), and much more precise! The
+//! > cosine version just needs the extra shift and wrap operations on x.
+//! >
+//! > Last but not least, I wouldn't be Nick if I didn't include the SIMD
+//! > optimized assembly version. It allows to perform the wrap operation very
+//! > efficiently so I'll give you the cosine:
+//! >
+//! > ```asm
+//! > // cos(x) = sin(x + pi/2)
+//! > addps xmm0, PI_2
+//! > movaps xmm1, xmm0
+//! > cmpnltps xmm1, PI
+//! > andps xmm1, PIx2
+//! > subps xmm0, xmm1
+//! >
+//! > // Parabola
+//! > movaps xmm1, xmm0
+//! > andps xmm1, abs
+//! > mulps xmm1, xmm0
+//! > mulps xmm0, B
+//! > mulps xmm1, C
+//! > addps xmm0, xmm1
+//! >
+//! > // Extra precision
+//! > movaps xmm1, xmm0
+//! > andps xmm1, abs
+//! > mulps xmm1, xmm0
+//! > subps xmm1, xmm0
+//! > mulps xmm1, P
+//! > addps xmm0, xmm1
+//! > ```
+//! >
+//! > This code computes four cosines in parallel, resulting in a peak
+//! > performance of about 9 clock cycles per cosine for most CPU architectures.
+//! > Sines take only 6 clock cycles ideally. The lower precision version can
+//! > even run at 3 clock cycles per sine... And lets not forget that all input
+//! > between -pi and pi is valid and the formula is exact at -pi, -pi/2, 0,
+//! > pi/2 and pi.
+//! >
+//! > So, the conclusion is don't ever again use a Taylor series to approximate
+//! > a sine or cosine! To add a useful discussion to this article, I'd love to
+//! > hear if anyone knows good approximations for other transcendental functions
+//! > like the exponential, logarithm and power functions.
+
+use super::F32;
+use core::f32::consts::FRAC_1_PI;
+
+impl F32 {
+    /// Approximates `cos(x)` in radians with a maximum error of `0.002`.
+    pub fn cos(self) -> Self {
+        let mut x = self;
+        x *= FRAC_1_PI / 2.0;
+        x -= 0.25 + (x + 0.25).floor().0;
+        x *= 16.0 * (x.abs() - 0.5);
+        x += 0.225 * x * (x.abs() - 1.0);
+        x
+    }
+}
+
+#[cfg(test)]
+pub(crate) mod tests {
+    use super::F32;
+
+    /// Maximum error in radians
+    pub(crate) const MAX_ERROR: f32 = 0.002;
+
+    /// Cosine test vectors - `(input, output)`
+    const TEST_VECTORS: &[(f32, f32)] = &[
+        (0.000, 1.000),
+        (0.140, 0.990),
+        (0.279, 0.961),
+        (0.419, 0.914),
+        (0.559, 0.848),
+        (0.698, 0.766),
+        (0.838, 0.669),
+        (0.977, 0.559),
+        (1.117, 0.438),
+        (1.257, 0.309),
+        (1.396, 0.174),
+        (1.536, 0.035),
+        (1.676, -0.105),
+        (1.815, -0.242),
+        (1.955, -0.375),
+        (2.094, -0.500),
+        (2.234, -0.616),
+        (2.374, -0.719),
+        (2.513, -0.809),
+        (2.653, -0.883),
+        (2.793, -0.940),
+        (2.932, -0.978),
+        (3.072, -0.998),
+        (3.211, -0.998),
+        (3.351, -0.978),
+        (3.491, -0.940),
+        (3.630, -0.883),
+        (3.770, -0.809),
+        (3.910, -0.719),
+        (4.049, -0.616),
+        (4.189, -0.500),
+        (4.328, -0.375),
+        (4.468, -0.242),
+        (4.608, -0.105),
+        (4.747, 0.035),
+        (4.887, 0.174),
+        (5.027, 0.309),
+        (5.166, 0.438),
+        (5.306, 0.559),
+        (5.445, 0.669),
+        (5.585, 0.766),
+        (5.725, 0.848),
+        (5.864, 0.914),
+        (6.004, 0.961),
+        (6.144, 0.990),
+        (6.283, 1.000),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        for &(x, expected) in TEST_VECTORS {
+            let cos_x = F32(x).cos();
+            let delta = (cos_x - expected).abs();
+
+            assert!(
+                delta <= MAX_ERROR,
+                "delta {} too large: {} vs {}",
+                delta,
+                cos_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,34 @@
+//! Calculates Euclidian division for a single-precision float.
+
+use super::F32;
+
+impl F32 {
+    /// Calculates Euclidean division, the matching method for `rem_euclid`.
+    pub fn div_euclid(self, rhs: Self) -> Self {
+        let q = (self / rhs).trunc();
+
+        if self % rhs >= Self::ZERO {
+            q
+        } else if rhs > Self::ZERO {
+            q - Self::ONE
+        } else {
+            q + Self::ONE
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        let a = F32(7.0);
+        let b = F32(4.0);
+
+        assert_eq!(a.div_euclid(b), F32(1.0));
+        assert_eq!((-a).div_euclid(b), F32(-2.0));
+        assert_eq!(a.div_euclid(-b), F32(-1.0));
+        assert_eq!((-a).div_euclid(-b), F32(2.0));
+    }
+}
@@ -0,0 +1,121 @@
+//! Exp approximation for a single-precision float.
+//!
+//! Method described at: <https://stackoverflow.com/a/6985769/2036035>
+
+use super::{EXPONENT_BIAS, F32};
+use core::f32::consts;
+
+impl F32 {
+    /// Returns `e^(self)`, (the exponential function).
+    #[inline]
+    pub fn exp(self) -> Self {
+        self.exp_ln2_approx(4)
+    }
+
+    /// Exp approximation for `f32`.
+    pub(crate) fn exp_ln2_approx(self, partial_iter: u32) -> Self {
+        if self == Self::ZERO {
+            return Self::ONE;
+        }
+
+        if (self - Self::ONE).abs() < f32::EPSILON {
+            return consts::E.into();
+        }
+
+        if (self - (-Self::ONE)).abs() < f32::EPSILON {
+            return Self::ONE / consts::E;
+        }
+
+        // log base 2(E) == 1/ln(2)
+        // x_fract + x_whole = x/ln2_recip
+        // ln2*(x_fract + x_whole) = x
+        let x_ln2recip = self * consts::LOG2_E;
+        let x_fract = x_ln2recip.fract();
+        let x_trunc = x_ln2recip.trunc();
+
+        //guaranteed to be 0 < x < 1.0
+        let x_fract = x_fract * consts::LN_2;
+        let fract_exp = x_fract.exp_smallx(partial_iter);
+
+        //need the 2^n portion, we can just extract that from the whole number exp portion
+        let fract_exponent: i32 = fract_exp
+            .extract_exponent_value()
+            .saturating_add(x_trunc.0 as i32);
+
+        if fract_exponent < -(EXPONENT_BIAS as i32) {
+            return Self::ZERO;
+        }
+
+        if fract_exponent > ((EXPONENT_BIAS + 1) as i32) {
+            return Self::INFINITY;
+        }
+
+        fract_exp.set_exponent(fract_exponent)
+    }
+
+    /// if x is between 0.0 and 1.0, we can approximate it with the a series
+    ///
+    /// Series from here:
+    /// <https://stackoverflow.com/a/6984495>
+    ///
+    /// e^x ~= 1 + x(1 + x/2(1 + (x?
+    #[inline]
+    pub(crate) fn exp_smallx(self, iter: u32) -> Self {
+        let mut total = 1.0;
+
+        for i in (1..=iter).rev() {
+            total = 1.0 + ((self.0 / (i as f32)) * total);
+        }
+
+        Self(total)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    pub(crate) const MAX_ERROR: f32 = 0.001;
+
+    /// exp test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
+        (1e-07, 1.0000001),
+        (1e-06, 1.000001),
+        (1e-05, 1.00001),
+        (1e-04, 1.0001),
+        (0.001, 1.0010005),
+        (0.01, 1.0100502),
+        (0.1, 1.105171),
+        (1.0, 2.7182817),
+        (10.0, 22026.465),
+        (-1e-08, 1.0),
+        (-1e-07, 0.9999999),
+        (-1e-06, 0.999999),
+        (-1e-05, 0.99999),
+        (-1e-04, 0.9999),
+        (-0.001, 0.9990005),
+        (-0.01, 0.99004984),
+        (-0.1, 0.9048374),
+        (-1.0, 0.36787945),
+        (-10.0, 4.539_993e-5),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32(-1000000.0).exp(), F32::ZERO);
+
+        for &(x, expected) in TEST_VECTORS {
+            let exp_x = F32(x).exp();
+            let relative_error = (exp_x - expected).abs() / expected;
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large for input {} : {} vs {}",
+                relative_error,
+                x,
+                exp_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,31 @@
+//! Floating point floor approximation for a single-precision float.
+
+use super::F32;
+
+impl F32 {
+    /// Returns the largest integer less than or equal to a number.
+    pub fn floor(self) -> Self {
+        let mut res = (self.0 as i32) as f32;
+
+        if self.0 < res {
+            res -= 1.0;
+        }
+
+        Self(res)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32(-1.1).floor().0, -2.0);
+        assert_eq!(F32(-0.1).floor().0, -1.0);
+        assert_eq!(F32(0.0).floor().0, 0.0);
+        assert_eq!(F32(1.0).floor().0, 1.0);
+        assert_eq!(F32(1.1).floor().0, 1.0);
+        assert_eq!(F32(2.9).floor().0, 2.0);
+    }
+}
@@ -0,0 +1,60 @@
+//! Floating point fractional number for a single-precision float.
+
+use super::F32;
+use super::{EXPONENT_BIAS, MANTISSA_BITS, MANTISSA_MASK};
+
+impl F32 {
+    /// Returns the fractional part of a number with sign.
+    pub fn fract(self) -> Self {
+        let x_bits = self.to_bits();
+        let exponent = self.extract_exponent_value();
+
+        // we know it is *only* fraction
+        if exponent < 0 {
+            return self;
+        }
+
+        // find the part of the fraction that would be left over
+        let fractional_part = x_bits.overflowing_shl(exponent as u32).0 & MANTISSA_MASK;
+
+        // if there isn't a fraction we can just return 0
+        if fractional_part == 0 {
+            // TODO: most people don't actually care about -0.0,
+            // so would it be better to just not copysign?
+            return Self(0.0).copysign(self);
+        }
+
+        // Note: alternatively this could use -1.0, but it's assumed subtraction would be more costly
+        // example: 'let new_exponent_bits = 127_u32.overflowing_shl(23_u32).0)) - 1.0'
+        let exponent_shift: u32 = (fractional_part.leading_zeros() - (32 - MANTISSA_BITS)) + 1;
+
+        let fractional_normalized: u32 =
+            fractional_part.overflowing_shl(exponent_shift).0 & MANTISSA_MASK;
+
+        let new_exponent_bits = (EXPONENT_BIAS - (exponent_shift))
+            .overflowing_shl(MANTISSA_BITS)
+            .0;
+
+        Self::from_bits(fractional_normalized | new_exponent_bits).copysign(self)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        // fraction check actually won't be the same, though technically exactly accurate
+        // so we test by adding back the number removed.
+        assert_eq!(F32(2.9).fract().0 + 2.0, 2.9);
+        assert_eq!(F32(-1.1).fract().0 - 1.0, -1.1);
+        assert_eq!(F32(-0.1).fract().0, -0.1);
+        assert_eq!(F32(0.0).fract().0, 0.0);
+        assert_eq!(F32(1.0).fract().0 + 1.0, 1.0);
+        assert_eq!(F32(1.1).fract().0 + 1.0, 1.1);
+
+        assert_eq!(F32(-100_000_000.13425345345).fract().0, -0.0);
+        assert_eq!(F32(100_000_000.13425345345).fract().0, 0.0);
+    }
+}
@@ -0,0 +1,23 @@
+//! Calculate length of the hypotenuse of a right triangle.
+
+use super::F32;
+
+impl F32 {
+    /// Calculate the length of the hypotenuse of a right-angle triangle.
+    pub fn hypot(self, rhs: Self) -> Self {
+        (self * self + rhs * rhs).sqrt()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        let x = F32(3.0);
+        let y = F32(4.0);
+        let difference = x.hypot(y) - F32(25.0).sqrt();
+        assert!(difference.abs() <= F32::EPSILON);
+    }
+}
@@ -0,0 +1,39 @@
+//! Fast approximation of `1/x`.
+//!
+//! Method described at: <https://bits.stephan-brumme.com/inverse.html>
+
+use super::F32;
+
+impl F32 {
+    /// Fast approximation of `1/x`.
+    pub fn inv(self) -> Self {
+        Self(f32::from_bits(0x7f00_0000 - self.0.to_bits()))
+    }
+}
+
+#[cfg(test)]
+pub(crate) mod tests {
+    use super::F32;
+
+    /// Deviation from the actual value (8%)
+    pub(crate) const MAX_ERROR: f32 = 0.08;
+
+    #[test]
+    fn sanity_check() {
+        for x in 0..100 {
+            let x = F32(x as f32);
+            let inv_x = x.inv().0;
+            let expected = 1.0 / x;
+            let allowed_delta = x * MAX_ERROR;
+            let actual_delta = inv_x - expected;
+
+            assert!(
+                actual_delta <= allowed_delta,
+                "delta {} too large: {} vs {}",
+                actual_delta,
+                inv_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,41 @@
+//! Inverse square root approximation function for a single-precision float.
+//!
+//! Method described at: <https://bits.stephan-brumme.com/invSquareRoot.html>
+
+use super::F32;
+
+impl F32 {
+    /// Approximate inverse square root with an average deviation of ~5%.
+    pub fn invsqrt(self) -> Self {
+        Self::from_bits(0x5f37_5a86 - (self.to_bits() >> 1))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+    use crate::float::sqrt::tests::TEST_VECTORS;
+
+    /// Deviation from the actual value (5%)
+    const MAX_ERROR: f32 = 0.05;
+
+    #[test]
+    fn sanity_check() {
+        for (x, expected) in TEST_VECTORS {
+            // The tests vectors are for sqrt(x), so invert the expected value
+            let expected = 1.0 / expected;
+
+            let invsqrt_x = F32(*x).invsqrt().0;
+            let allowed_delta = x * MAX_ERROR;
+            let actual_delta = invsqrt_x - expected;
+
+            assert!(
+                actual_delta <= allowed_delta,
+                "delta {} too large: {} vs {}",
+                actual_delta,
+                invsqrt_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,121 @@
+//! Natural log (ln) approximation for a single-precision float.
+//!
+//! Method described at: <https://stackoverflow.com/a/44232045/>
+//!
+//! Modified to not be restricted to int range and only values of x above 1.0.
+//! Also got rid of most of the slow conversions. Should work for all positive values of x.
+
+use super::{EXPONENT_MASK, F32};
+use core::f32::consts::LN_2;
+
+impl F32 {
+    /// Approximates the natural logarithm of the number.
+    // Note: excessive precision ignored because it hides the origin of the numbers used for the
+    // ln(1.0->2.0) polynomial
+    #[allow(clippy::excessive_precision)]
+    pub fn ln(self) -> Self {
+        // x may essentially be 1.0 but, as clippy notes, these kinds of
+        // floating point comparisons can fail when the bit pattern is not the sames
+        if (self - Self::ONE).abs() < f32::EPSILON {
+            return Self::ZERO;
+        }
+
+        let x_less_than_1 = self < 1.0;
+
+        // Note: we could use the fast inverse approximation here found in super::inv::inv_approx, but
+        // the precision of such an approximation is assumed not good enough.
+        let x_working = if x_less_than_1 { self.inv() } else { self };
+
+        // according to the SO post ln(x) = ln((2^n)*y)= ln(2^n) + ln(y) = ln(2) * n + ln(y)
+        // get exponent value
+        let base2_exponent = x_working.extract_exponent_value() as u32;
+        let divisor = f32::from_bits(x_working.to_bits() & EXPONENT_MASK);
+
+        // supposedly normalizing between 1.0 and 2.0
+        let x_working = x_working / divisor;
+
+        // approximate polynomial generated from maple in the post using Remez Algorithm:
+        // https://en.wikipedia.org/wiki/Remez_algorithm
+        let ln_1to2_polynomial = -1.741_793_9
+            + (2.821_202_6
+                + (-1.469_956_8 + (0.447_179_55 - 0.056_570_851 * x_working) * x_working)
+                    * x_working)
+                * x_working;
+
+        // ln(2) * n + ln(y)
+        let result = (base2_exponent as f32) * LN_2 + ln_1to2_polynomial;
+
+        if x_less_than_1 {
+            -result
+        } else {
+            result
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    pub(crate) const MAX_ERROR: f32 = 0.001;
+
+    /// ln(x) test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
+        (1e-20, -46.0517),
+        (1e-19, -43.749115),
+        (1e-18, -41.446533),
+        (1e-17, -39.143948),
+        (1e-16, -36.841362),
+        (1e-15, -34.538776),
+        (1e-14, -32.23619),
+        (1e-13, -29.933607),
+        (1e-12, -27.631021),
+        (1e-11, -25.328436),
+        (1e-10, -23.02585),
+        (1e-09, -20.723267),
+        (1e-08, -18.420681),
+        (1e-07, -16.118095),
+        (1e-06, -13.815511),
+        (1e-05, -11.512925),
+        (1e-04, -9.2103405),
+        (0.001, -6.9077554),
+        (0.01, -4.6051702),
+        (0.1, -2.3025851),
+        (10.0, 2.3025851),
+        (100.0, 4.6051702),
+        (1000.0, 6.9077554),
+        (10000.0, 9.2103405),
+        (100000.0, 11.512925),
+        (1000000.0, 13.815511),
+        (10000000.0, 16.118095),
+        (100000000.0, 18.420681),
+        (1000000000.0, 20.723267),
+        (10000000000.0, 23.02585),
+        (100000000000.0, 25.328436),
+        (1000000000000.0, 27.631021),
+        (10000000000000.0, 29.933607),
+        (100000000000000.0, 32.23619),
+        (1000000000000000.0, 34.538776),
+        (1e+16, 36.841362),
+        (1e+17, 39.143948),
+        (1e+18, 41.446533),
+        (1e+19, 43.749115),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32::ONE.ln(), F32::ZERO);
+        for &(x, expected) in TEST_VECTORS {
+            let ln_x = F32(x).ln();
+            let relative_error = (ln_x - expected).abs() / expected;
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large: {} vs {}",
+                relative_error,
+                ln_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,192 @@
+//! log_b(a) approximation for a single-precision float.
+
+use super::F32;
+
+impl F32 {
+    /// Approximates the logarithm of the number with respect to an arbitrary base.
+    pub fn log(self, base: Self) -> Self {
+        (Self::ONE / base.ln()) * self.ln()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    pub(crate) const MAX_ERROR: f32 = 0.001;
+
+    /// log3(x) test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS_BASE3: &[(f32, f32)] = &[
+        (1e-20, -41.918_064),
+        (1e-19, -39.822_16),
+        (1e-18, -37.726_26),
+        (1e-17, -35.630_356),
+        (1e-16, -33.534_454),
+        (1e-15, -31.438_549),
+        (1e-14, -29.342_646),
+        (1e-13, -27.246_744),
+        (1e-12, -25.150_839),
+        (1e-11, -23.054_935),
+        (1e-10, -20.959_032),
+        (1e-09, -18.863_13),
+        (1e-08, -16.767_227),
+        (1e-07, -14.671_323),
+        (1e-06, -12.575_419),
+        (1e-05, -10.479_516),
+        (1e-04, -8.383_614),
+        (0.001, -6.287_709_7),
+        (0.01, -4.191_807),
+        (0.1, -2.095_903_4),
+        (10.0, 2.095_903_4),
+        (100.0, 4.191_807),
+        (1000.0, 6.287_709_7),
+        (10000.0, 8.383_614),
+        (100000.0, 10.479_516),
+        (1000000.0, 12.575_419),
+        (10000000.0, 14.671_323),
+        (100000000.0, 16.767_227),
+        (1000000000.0, 18.863_13),
+        (10000000000.0, 20.959_032),
+        (100000000000.0, 23.054_935),
+        (1000000000000.0, 25.150_839),
+        (10000000000000.0, 27.246_744),
+        (100000000000000.0, 29.342_646),
+        (1000000000000000.0, 31.438_549),
+        (1e+16, 33.534_454),
+        (1e+17, 35.630_356),
+        (1e+18, 37.726_26),
+        (1e+19, 39.822_16),
+    ];
+
+    /// log5.5(x) test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS_BASE5_5: &[(f32, f32)] = &[
+        (1e-20, -27.013_786),
+        (1e-19, -25.663_097),
+        (1e-18, -24.312_408),
+        (1e-17, -22.961_72),
+        (1e-16, -21.611_03),
+        (1e-15, -20.260_34),
+        (1e-14, -18.909_65),
+        (1e-13, -17.558_962),
+        (1e-12, -16.208_273),
+        (1e-11, -14.857_583),
+        (1e-10, -13.506_893),
+        (1e-09, -12.156_204),
+        (1e-08, -10.805_515),
+        (1e-07, -9.454_825),
+        (1e-06, -8.104_136),
+        (1e-05, -6.753_446_6),
+        (1e-04, -5.402_757_6),
+        (0.001, -4.052_068),
+        (0.01, -2.701_378_8),
+        (0.1, -1.350_689_4),
+        (10.0, 1.350_689_4),
+        (100.0, 2.701_378_8),
+        (1000.0, 4.052_068),
+        (10000.0, 5.402_757_6),
+        (100000.0, 6.753_446_6),
+        (1000000.0, 8.104_136),
+        (10000000.0, 9.454_825),
+        (100000000.0, 10.805_515),
+        (1000000000.0, 12.156_204),
+        (10000000000.0, 13.506_893),
+        (100000000000.0, 14.857_583),
+        (1000000000000.0, 16.208_273),
+        (10000000000000.0, 17.558_962),
+        (100000000000000.0, 18.909_65),
+        (1000000000000000.0, 20.260_34),
+        (1e+16, 21.611_03),
+        (1e+17, 22.961_72),
+        (1e+18, 24.312_408),
+        (1e+19, 25.663_097),
+    ];
+
+    /// log12.7(x) test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS_BASE12_7: &[(f32, f32)] = &[
+        (1e-20, -18.119_164),
+        (1e-19, -17.213_205),
+        (1e-18, -16.307_247),
+        (1e-17, -15.401_289),
+        (1e-16, -14.495_331),
+        (1e-15, -13.589_373),
+        (1e-14, -12.683_414),
+        (1e-13, -11.777_456),
+        (1e-12, -10.871_498),
+        (1e-11, -9.965_54),
+        (1e-10, -9.059_582),
+        (1e-09, -8.153_624),
+        (1e-08, -7.247_665_4),
+        (1e-07, -6.341_707),
+        (1e-06, -5.435_749),
+        (1e-05, -4.529_791),
+        (1e-04, -3.623_832_7),
+        (0.001, -2.717_874_5),
+        (0.01, -1.811_916_4),
+        (0.1, -0.905_958_2),
+        (10.0, 0.905_958_2),
+        (100.0, 1.811_916_4),
+        (1000.0, 2.717_874_5),
+        (10000.0, 3.623_832_7),
+        (100000.0, 4.529_791),
+        (1000000.0, 5.435_749),
+        (10000000.0, 6.341_707),
+        (100000000.0, 7.247_665_4),
+        (1000000000.0, 8.153_624),
+        (10000000000.0, 9.059_582),
+        (100000000000.0, 9.965_54),
+        (1000000000000.0, 10.871_498),
+        (10000000000000.0, 11.777_456),
+        (100000000000000.0, 12.683_414),
+        (1000000000000000.0, 13.589_373),
+        (1e+16, 14.495_331),
+        (1e+17, 15.401_289),
+        (1e+18, 16.307_247),
+        (1e+19, 17.213_205),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32::ONE.log(F32(3.0)), F32::ZERO);
+        assert_eq!(F32::ONE.log(F32(5.5)), F32::ZERO);
+        assert_eq!(F32::ONE.log(F32(12.7)), F32::ZERO);
+
+        for &(x, expected) in TEST_VECTORS_BASE3 {
+            let log_x = F32(x).log(F32(3.0));
+            let relative_error = (log_x - expected).abs() / expected;
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large: {} vs {}",
+                relative_error,
+                log_x,
+                expected
+            );
+        }
+
+        for &(x, expected) in TEST_VECTORS_BASE5_5 {
+            let log_x = F32(x).log(F32(5.5));
+            let relative_error = (log_x - expected).abs() / expected;
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large: {} vs {}",
+                relative_error,
+                log_x,
+                expected
+            );
+        }
+
+        for &(x, expected) in TEST_VECTORS_BASE12_7 {
+            let log_x = F32(x).log(F32(12.7));
+            let relative_error = (log_x - expected).abs() / expected;
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large: {} vs {}",
+                relative_error,
+                log_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,78 @@
+//! log base 2 approximation for a single-precision float.
+
+use super::F32;
+use core::f32::consts::LOG10_E;
+
+impl F32 {
+    /// Approximates the base 10 logarithm of the number.
+    pub fn log10(self) -> Self {
+        self.ln() * LOG10_E
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+    pub(crate) const MAX_ERROR: f32 = 0.001;
+
+    /// log10(x) test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
+        (1e-20, -20.0),
+        (1e-19, -19.0),
+        (1e-18, -18.0),
+        (1e-17, -17.0),
+        (1e-16, -16.0),
+        (1e-15, -15.0),
+        (1e-14, -14.0),
+        (1e-13, -13.0),
+        (1e-12, -12.0),
+        (1e-11, -11.0),
+        (1e-10, -10.0),
+        (1e-09, -9.0),
+        (1e-08, -8.0),
+        (1e-07, -7.0),
+        (1e-06, -6.0),
+        (1e-05, -5.0),
+        (1e-04, -4.0),
+        (0.001, -3.0),
+        (0.01, -2.0),
+        (0.1, -1.0),
+        (10.0, 1.0),
+        (100.0, 2.0),
+        (1000.0, 3.0),
+        (10000.0, 4.0),
+        (100000.0, 5.0),
+        (1000000.0, 6.0),
+        (10000000.0, 7.0),
+        (100000000.0, 8.0),
+        (1000000000.0, 9.0),
+        (10000000000.0, 10.0),
+        (100000000000.0, 11.0),
+        (1000000000000.0, 12.0),
+        (10000000000000.0, 13.0),
+        (100000000000000.0, 14.0),
+        (1000000000000000.0, 15.0),
+        (1e+16, 16.0),
+        (1e+17, 17.0),
+        (1e+18, 18.0),
+        (1e+19, 19.0),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32::ONE.log10(), F32::ZERO);
+
+        for &(x, expected) in TEST_VECTORS {
+            let ln_x = F32(x).log10();
+            let relative_error = (ln_x - expected).abs() / expected;
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large: {} vs {}",
+                relative_error,
+                ln_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,79 @@
+//! log base 2 approximation for a single-precision float.
+
+use super::F32;
+use core::f32::consts::LOG2_E;
+
+impl F32 {
+    /// Approximates the base 2 logarithm of the number.
+    pub fn log2(self) -> Self {
+        self.ln() * LOG2_E
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    pub(crate) const MAX_ERROR: f32 = 0.001;
+
+    /// log2(x) test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
+        (1e-20, -66.43856),
+        (1e-19, -63.116634),
+        (1e-18, -59.794704),
+        (1e-17, -56.47278),
+        (1e-16, -53.15085),
+        (1e-15, -49.828922),
+        (1e-14, -46.506992),
+        (1e-13, -43.185066),
+        (1e-12, -39.863136),
+        (1e-11, -36.54121),
+        (1e-10, -33.21928),
+        (1e-09, -29.897352),
+        (1e-08, -26.575424),
+        (1e-07, -23.253496),
+        (1e-06, -19.931568),
+        (1e-05, -16.60964),
+        (1e-04, -13.287712),
+        (0.001, -9.965784),
+        (0.01, -6.643856),
+        (0.1, -3.321928),
+        (10.0, 3.321928),
+        (100.0, 6.643856),
+        (1000.0, 9.965784),
+        (10000.0, 13.287712),
+        (100000.0, 16.60964),
+        (1000000.0, 19.931568),
+        (10000000.0, 23.253496),
+        (100000000.0, 26.575424),
+        (1000000000.0, 29.897352),
+        (10000000000.0, 33.21928),
+        (100000000000.0, 36.54121),
+        (1000000000000.0, 39.863136),
+        (10000000000000.0, 43.185066),
+        (100000000000000.0, 46.506992),
+        (1000000000000000.0, 49.828922),
+        (1e+16, 53.15085),
+        (1e+17, 56.47278),
+        (1e+18, 59.794704),
+        (1e+19, 63.116634),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32::ONE.log2(), F32::ZERO);
+
+        for &(x, expected) in TEST_VECTORS {
+            let ln_x = F32(x).log2().0;
+            let relative_error = (ln_x - expected).abs() / expected;
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large: {} vs {}",
+                relative_error,
+                ln_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,22 @@
+//! Fused multiply-add. Computes `(self * a) + b`
+
+use super::F32;
+
+impl F32 {
+    /// Computes `(self * a) + b`.
+    pub fn mul_add(self, a: Self, b: Self) -> Self {
+        self * a + b
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32(0.0).mul_add(F32(0.0), F32(1.0)), F32(1.0));
+        assert_eq!(F32(1.0).mul_add(F32(2.0), F32(3.5)), F32(5.5));
+        assert_eq!(F32(1.0).mul_add(F32(-1.0), F32(0.0)), F32(-1.0));
+    }
+}
@@ -0,0 +1,209 @@
+//! `x^n` with fractional `n` approximation for a single-precision float.
+
+use super::F32;
+
+impl F32 {
+    /// Approximates a number raised to a floating point power.
+    pub fn powf(self, n: Self) -> Self {
+        // using x^n = exp(ln(x^n)) = exp(n*ln(x))
+        if self >= Self::ZERO {
+            (n * self.ln()).exp()
+        } else if !n.is_integer() {
+            Self::NAN
+        } else if n.is_even() {
+            // if n is even, then we know that the result will have no sign, so we can remove it
+            (n * self.without_sign().ln()).exp()
+        } else {
+            // if n isn't even, we need to multiply by -1.0 at the end.
+            -(n * self.without_sign().ln()).exp()
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    /// error builds up from both exp and ln approximation, so we double the error allowed.
+    pub(crate) const MAX_ERROR: f32 = 0.002;
+
+    ///  powf(3,x) test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS_POW3: &[(f32, f32)] = &[
+        (-1e-20, 1.0),
+        (-1e-19, 1.0),
+        (-1e-18, 1.0),
+        (-1e-17, 1.0),
+        (-1e-16, 0.9999999999999999),
+        (-1e-15, 0.9999999999999989),
+        (-1e-14, 0.999999999999989),
+        (-1e-13, 0.9999999999998901),
+        (-1e-12, 0.9999999999989014),
+        (-1e-11, 0.9999999999890139),
+        (-1e-10, 0.9999999998901388),
+        (-1e-09, 0.9999999989013877),
+        (-1e-08, 0.9999999890138772),
+        (-1e-07, 0.999_999_9),
+        (-1e-06, 0.999_998_9),
+        (-1e-05, 0.999_989_03),
+        (-1e-04, 0.999_890_15),
+        (-0.001, 0.998_901_96),
+        (-0.01, 0.989_074),
+        (-0.1, 0.895_958_5),
+        (-1.0, 0.333_333_34),
+        (-10.0, 1.693_508_8e-5),
+        (-100.0, 0e0),
+        (-1000.0, 0.0),
+        (1e-20, 1.0),
+        (1e-19, 1.0),
+        (1e-18, 1.0),
+        (1e-17, 1.0),
+        (1e-16, 1.0),
+        (1e-15, 1.000000000000001),
+        (1e-14, 1.0000000000000109),
+        (1e-13, 1.00000000000011),
+        (1e-12, 1.0000000000010987),
+        (1e-11, 1.000000000010986),
+        (1e-10, 1.0000000001098612),
+        (1e-09, 1.0000000010986123),
+        (1e-08, 1.000000010986123),
+        (1e-07, 1.000_000_1),
+        (1e-06, 1.000_001_1),
+        (1e-05, 1.000_011),
+        (1e-04, 1.000_109_9),
+        (0.001, 1.001_099_2),
+        (0.01, 1.011_046_6),
+        (0.1, 1.116_123_2),
+        (1.0, 3.0),
+        (10.0, 59049.0),
+    ];
+
+    ///  powf(150,x) test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS_POW150: &[(f32, f32)] = &[
+        (-1e-20, 1.0),
+        (-1e-19, 1.0),
+        (-1e-18, 1.0),
+        (-1e-17, 1.0),
+        (-1e-16, 0.9999999999999994),
+        (-1e-15, 0.999999999999995),
+        (-1e-14, 0.9999999999999499),
+        (-1e-13, 0.999999999999499),
+        (-1e-12, 0.9999999999949893),
+        (-1e-11, 0.9999999999498936),
+        (-1e-10, 0.9999999994989365),
+        (-1e-09, 0.9999999949893649),
+        (-1e-08, 0.999_999_94),
+        (-1e-07, 0.999_999_5),
+        (-1e-06, 0.999_995),
+        (-1e-05, 0.999_949_9),
+        (-1e-04, 0.999_499_1),
+        (-0.001, 0.995_001_9),
+        (-0.01, 0.951_128_24),
+        (-0.1, 0.605_885_9),
+        (-1.0, 0.006_666_667),
+        (-10.0, 1.734_153e-22),
+        (-100.0, 0e0),
+        (-1000.0, 0.0),
+        (-10000.0, 0.0),
+        (-100000.0, 0.0),
+        (-1000000.0, 0.0),
+        (-10000000.0, 0.0),
+        (-100000000.0, 0.0),
+        (-1000000000.0, 0.0),
+        (-10000000000.0, 0.0),
+        (-100000000000.0, 0.0),
+        (-1000000000000.0, 0.0),
+        (-10000000000000.0, 0.0),
+        (-100000000000000.0, 0.0),
+        (-1000000000000000.0, 0.0),
+        (-1e+16, 0.0),
+        (-1e+17, 0.0),
+        (-1e+18, 0.0),
+        (-1e+19, 0.0),
+        (1e-20, 1.0),
+        (1e-19, 1.0),
+        (1e-18, 1.0),
+        (1e-17, 1.0),
+        (1e-16, 1.0000000000000004),
+        (1e-15, 1.000000000000005),
+        (1e-14, 1.0000000000000502),
+        (1e-13, 1.0000000000005012),
+        (1e-12, 1.0000000000050107),
+        (1e-11, 1.0000000000501064),
+        (1e-10, 1.0000000005010636),
+        (1e-09, 1.0000000050106352),
+        (1e-08, 1.000000050106354),
+        (1e-07, 1.000_000_5),
+        (1e-06, 1.000_005),
+        (1e-05, 1.000_050_1),
+        (1e-04, 1.000_501_2),
+        (0.001, 1.005_023_2),
+        (0.01, 1.051_382_9),
+        (0.1, 1.650_475_6),
+        (1.0, 150.0),
+        (10.0, 5.766_504e21),
+    ];
+
+    /// misc powf(x,n) test vectors - `(base_input, power_input, output)`
+    pub(crate) const TEST_VECTORS_MISC: &[(f32, f32, f32)] = &[
+        (-0.5881598, 2.0, 0.345_931_95),
+        (-0.5881598, 3.2, f32::NAN),
+        (-0.5881598, 3.0, -0.203_463_27),
+        (-1000000.0, 4.0, 1e+24),
+    ];
+
+    fn calc_relative_error(experimental: F32, expected: f32) -> F32 {
+        if experimental.is_nan() && expected.is_nan() {
+            F32::ZERO
+        } else if expected != 0.0 {
+            (experimental - expected) / expected
+        } else {
+            (experimental - expected) / (expected + 1.0e-20)
+        }
+    }
+
+    #[test]
+    fn sanity_check() {
+        for &(x, expected) in TEST_VECTORS_POW3 {
+            let exp_x = F32(3.0).powf(F32(x));
+            let relative_error = calc_relative_error(exp_x, expected);
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large for input {} : {} vs {}",
+                relative_error,
+                x,
+                exp_x,
+                expected
+            );
+        }
+
+        for &(x, expected) in TEST_VECTORS_POW150 {
+            let exp_x = F32(150.0).powf(F32(x));
+            let relative_error = calc_relative_error(exp_x, expected);
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large for input {} : {} vs {}",
+                relative_error,
+                x,
+                exp_x,
+                expected
+            );
+        }
+
+        for &(base_input, power_input, expected) in TEST_VECTORS_MISC {
+            let exp_x = F32(base_input).powf(F32(power_input));
+            let relative_error = calc_relative_error(exp_x, expected);
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large for input {}.powf({}) : {} vs {}",
+                relative_error,
+                base_input,
+                power_input,
+                exp_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,59 @@
+//! Takes the reciprocal (inverse) of a number, `1/x`.
+
+use super::F32;
+
+impl F32 {
+    /// Returns the reciprocal (inverse) of a number, `1/x`.
+    pub fn recip(self) -> Self {
+        let mut x = self;
+
+        let sx = if x < 0.0 { F32(-1.0) } else { F32(1.0) };
+        x *= sx;
+
+        let mut v = F32(f32::from_bits(0x7EF1_27EAu32.wrapping_sub(x.to_bits())));
+        let w = x * v;
+
+        // v.0 *= 2.0 - w;
+        // v.0 *= 4.0 + w * (-6.0 + w * (4.0 - w));
+        v.0 *=
+            8.0 + w * (-28.0 + w * (56.0 + w * (-70.0 + w * (56.0 + w * (-28.0 + w * (8.0 - w))))));
+
+        v.0 * sx
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    pub(crate) const MAX_ERROR: f32 = 1e-5;
+
+    pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
+        (0.00001, 100000.0),
+        (1.0, 1.0),
+        (2.0, 0.5),
+        (0.25, 4.0),
+        (-0.5, -2.0),
+        (core::f32::consts::PI, 1.0 / core::f32::consts::PI),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32(0.0).recip(), F32(core::f32::INFINITY));
+        assert_eq!(F32(-0.0).recip(), F32(core::f32::NEG_INFINITY));
+
+        for &(x, expected) in TEST_VECTORS {
+            let recip_x = F32(x).recip();
+            let relative_error = (recip_x - expected).abs() / expected;
+
+            assert!(
+                relative_error <= MAX_ERROR,
+                "relative_error {} too large for input {} : {} vs {}",
+                relative_error,
+                x,
+                recip_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,32 @@
+//! Calculate Euclidian remainder for a single-precision float.
+
+use super::F32;
+
+impl F32 {
+    /// Calculates the least non-negative remainder of `self (mod rhs)`.
+    pub fn rem_euclid(self, rhs: Self) -> Self {
+        let r = self % rhs;
+
+        if r >= Self::ZERO {
+            r
+        } else {
+            r + rhs.abs()
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        let a = F32(7.0);
+        let b = F32(4.0);
+
+        assert_eq!(a.rem_euclid(b), F32(3.0));
+        assert_eq!((-a).rem_euclid(b), F32(1.0));
+        assert_eq!(a.rem_euclid(-b), F32(3.0));
+        assert_eq!((-a).rem_euclid(-b), F32(1.0));
+    }
+}
@@ -0,0 +1,33 @@
+//! Round a single-precision float.
+
+use super::F32;
+
+impl F32 {
+    /// Returns the nearest integer to a number.
+    pub fn round(self) -> Self {
+        Self(((self.0 + Self(0.5).copysign(self).0) as i32) as f32)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32(0.0).round(), F32(0.0));
+        assert_eq!(F32(-0.0).round(), F32(-0.0));
+
+        assert_eq!(F32(0.49999).round(), F32(0.0));
+        assert_eq!(F32(-0.49999).round(), F32(-0.0));
+
+        assert_eq!(F32(0.5).round(), F32(1.0));
+        assert_eq!(F32(-0.5).round(), F32(-1.0));
+
+        assert_eq!(F32(9999.499).round(), F32(9999.0));
+        assert_eq!(F32(-9999.499).round(), F32(-9999.0));
+
+        assert_eq!(F32(9999.5).round(), F32(10000.0));
+        assert_eq!(F32(-9999.5).round(), F32(-10000.0));
+    }
+}
@@ -0,0 +1,33 @@
+//! Returns a number that represents the sign of `self`.
+
+use super::F32;
+
+impl F32 {
+    /// Returns a number that represents the sign of `self`.
+    ///
+    /// * `1.0` if the number is positive, `+0.0` or `INFINITY`
+    /// * `-1.0` if the number is negative, `-0.0` or `NEG_INFINITY`
+    /// * `NAN` if the number is `NAN`
+    pub fn signum(self) -> Self {
+        if self.is_nan() {
+            Self::NAN
+        } else {
+            F32(1.0).copysign(self)
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32::INFINITY.signum(), F32(1.0));
+        assert_eq!(F32(0.0).signum(), F32(1.0));
+        assert_eq!(F32(1.0).signum(), F32(1.0));
+        assert_eq!(F32::NEG_INFINITY.signum(), F32(-1.0));
+        assert_eq!(F32(-0.0).signum(), F32(-1.0));
+        assert_eq!(F32(-1.0).signum(), F32(-1.0));
+    }
+}
@@ -0,0 +1,83 @@
+//! Sine approximation, implemented in terms of `cos(x)`.
+
+use super::F32;
+use core::f32::consts::PI;
+
+impl F32 {
+    /// Approximates `sin(x)` in radians with a maximum error of `0.002`.
+    pub fn sin(self) -> Self {
+        (self - PI / 2.0).cos()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+    use crate::float::cos::tests::MAX_ERROR;
+
+    /// Sine test vectors - `(input, output)`
+    const TEST_VECTORS: &[(f32, f32)] = &[
+        (0.000, 0.000),
+        (0.140, 0.139),
+        (0.279, 0.276),
+        (0.419, 0.407),
+        (0.559, 0.530),
+        (0.698, 0.643),
+        (0.838, 0.743),
+        (0.977, 0.829),
+        (1.117, 0.899),
+        (1.257, 0.951),
+        (1.396, 0.985),
+        (1.536, 0.999),
+        (1.676, 0.995),
+        (1.815, 0.970),
+        (1.955, 0.927),
+        (2.094, 0.866),
+        (2.234, 0.788),
+        (2.374, 0.695),
+        (2.513, 0.588),
+        (2.653, 0.469),
+        (2.793, 0.342),
+        (2.932, 0.208),
+        (3.072, 0.070),
+        (3.211, -0.070),
+        (3.351, -0.208),
+        (3.491, -0.342),
+        (3.630, -0.469),
+        (3.770, -0.588),
+        (3.910, -0.695),
+        (4.049, -0.788),
+        (4.189, -0.866),
+        (4.328, -0.927),
+        (4.468, -0.970),
+        (4.608, -0.995),
+        (4.747, -0.999),
+        (4.887, -0.985),
+        (5.027, -0.951),
+        (5.166, -0.899),
+        (5.306, -0.829),
+        (5.445, -0.743),
+        (5.585, -0.643),
+        (5.725, -0.530),
+        (5.864, -0.407),
+        (6.004, -0.276),
+        (6.144, -0.139),
+        (6.283, 0.000),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        for &(x, expected) in TEST_VECTORS {
+            let sin_x = F32(x).sin();
+            let delta = (sin_x - expected).abs();
+
+            assert!(
+                delta <= MAX_ERROR,
+                "delta {} too large: {} vs {}",
+                delta,
+                sin_x,
+                expected
+            );
+        }
+    }
+}
@@ -0,0 +1,33 @@
+//! Simultaneously computes the sine and cosine of the number, `x`.
+
+use super::F32;
+
+impl F32 {
+    /// Simultaneously computes the sine and cosine of the number, `x`.
+    /// Returns `(sin(x), cos(x))`.
+    pub fn sin_cos(self) -> (Self, Self) {
+        (self.sin(), self.cos())
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    const TEST_VECTORS: &[f32] = &[
+        0.000, 0.140, 0.279, 0.419, 0.559, 0.698, 0.838, 0.977, 1.117, 1.257, 1.396, 1.536, 1.676,
+        1.815, 1.955, 2.094, 2.234, 2.374, 2.513, 2.653, 2.793, 2.932, 3.072, 3.211, 3.351, 3.491,
+        3.630, 3.770, 3.910, 4.049, 4.189, 4.328, 4.468, 4.608, 4.747, 4.887, 5.027, 5.166, 5.306,
+        5.445, 5.585, 5.725, 5.864, 6.004, 6.144, 6.283,
+    ];
+
+    #[test]
+    fn sanity_check() {
+        for &x in TEST_VECTORS {
+            let sin_x = F32(x).sin();
+            let cos_x = F32(x).cos();
+
+            assert_eq!(F32(x).sin_cos(), (sin_x, cos_x));
+        }
+    }
+}
@@ -0,0 +1,71 @@
+//! Square root approximation function for a single-precision float.
+//!
+//! Method described at: <https://bits.stephan-brumme.com/squareRoot.html>
+
+use super::F32;
+
+impl F32 {
+    /// Approximates the square root of a number with an average deviation of ~5%.
+    ///
+    /// Returns [`Self::NAN`] if `self` is a negative number.
+    pub fn sqrt(self) -> Self {
+        if self >= Self::ZERO {
+            Self::from_bits((self.to_bits() + 0x3f80_0000) >> 1)
+        } else {
+            Self::NAN
+        }
+    }
+}
+
+#[cfg(test)]
+pub(crate) mod tests {
+    use super::F32;
+
+    /// Deviation from the actual value (5%)
+    pub(crate) const MAX_ERROR: f32 = 0.05;
+
+    /// Square root test vectors - `(input, output)`
+    pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
+        (1.0, 1.0),
+        (2.0, 1.414),
+        (3.0, 1.732),
+        (4.0, 2.0),
+        (5.0, 2.236),
+        (10.0, 3.162),
+        (100.0, 10.0),
+        (250.0, 15.811),
+        (500.0, 22.36),
+        (1000.0, 31.622),
+        (2500.0, 50.0),
+        (5000.0, 70.710),
+        (1000000.0, 1000.0),
+        (2500000.0, 1581.138),
+        (5000000.0, 2236.067),
+        (10000000.0, 3162.277),
+        (25000000.0, 5000.0),
+        (50000000.0, 7071.067),
+        (100000000.0, 10000.0),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        for &(x, expected) in TEST_VECTORS {
+            let sqrt_x = F32(x).sqrt();
+            let allowed_delta = x * MAX_ERROR;
+            let actual_delta = sqrt_x - expected;
+
+            assert!(
+                actual_delta <= allowed_delta,
+                "delta {} too large: {} vs {}",
+                actual_delta,
+                sqrt_x,
+                expected
+            );
+        }
+    }
+
+    #[test]
+    fn negative_is_nan() {
+        assert!(F32(-1.0).sqrt().is_nan());
+    }
+}
@@ -0,0 +1,95 @@
+//! Tangent approximation for a single-precision float.
+
+use super::F32;
+
+impl F32 {
+    /// Approximates `tan(x)` in radians with a maximum error of `0.6`.
+    pub fn tan(self) -> Self {
+        self.sin() / self.cos()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    /// Maximum error in radians
+    // TODO(tarcieri): this is kinda bad, find a better approximation?
+    const MAX_ERROR: f32 = 0.6;
+
+    /// Tangent test vectors - `(input, output)`
+    const TEST_VECTORS: &[(f32, f32)] = &[
+        (0.000, 0.000),
+        (0.140, 0.141),
+        (0.279, 0.287),
+        (0.419, 0.445),
+        (0.559, 0.625),
+        (0.698, 0.839),
+        (0.838, 1.111),
+        (0.977, 1.483),
+        (1.117, 2.050),
+        (1.257, 3.078),
+        (1.396, 5.671),
+        (1.536, 28.636),
+        (1.676, -9.514),
+        (1.815, -4.011),
+        (1.955, -2.475),
+        (2.094, -1.732),
+        (2.234, -1.280),
+        (2.374, -0.966),
+        (2.513, -0.727),
+        (2.653, -0.532),
+        (2.793, -0.364),
+        (2.932, -0.213),
+        (3.072, -0.070),
+        (3.211, 0.070),
+        (3.351, 0.213),
+        (3.491, 0.364),
+        (3.630, 0.532),
+        (3.770, 0.727),
+        (3.910, 0.966),
+        (4.049, 1.280),
+        (4.189, 1.732),
+        (4.328, 2.475),
+        (4.468, 4.011),
+        (4.608, 9.514),
+        (4.747, -28.636),
+        (4.887, -5.671),
+        (5.027, -3.078),
+        (5.166, -2.050),
+        (5.306, -1.483),
+        (5.445, -1.111),
+        (5.585, -0.839),
+        (5.725, -0.625),
+        (5.864, -0.445),
+        (6.004, -0.287),
+        (6.144, -0.141),
+        (6.283, 0.000),
+    ];
+
+    #[test]
+    fn sanity_check() {
+        for &(x, expected) in TEST_VECTORS {
+            let tan_x = F32(x).tan();
+            let delta = (tan_x - expected).abs();
+
+            assert!(
+                delta <= MAX_ERROR,
+                "delta {} too large: {} vs {}",
+                delta,
+                tan_x,
+                expected
+            );
+        }
+    }
+
+    #[test]
+    fn zero() {
+        assert_eq!(F32::ZERO.tan(), F32::ZERO);
+    }
+
+    #[test]
+    fn nan() {
+        assert!(F32::NAN.tan().is_nan());
+    }
+}
@@ -0,0 +1,48 @@
+//! Floating point whole number for a single-precision float.
+
+use super::{F32, MANTISSA_MASK};
+
+impl F32 {
+    /// Returns the integer part of a number.
+    pub fn trunc(self) -> Self {
+        let x_bits = self.to_bits();
+        let exponent = self.extract_exponent_value();
+
+        // exponent is negative, there is no whole number, just return zero
+        if exponent < 0 {
+            return F32::ZERO.copysign(self);
+        }
+
+        let exponent_clamped = i32::max(exponent, 0) as u32;
+
+        // find the part of the fraction that would be left over
+        let fractional_part = x_bits.overflowing_shl(exponent_clamped).0 & MANTISSA_MASK;
+
+        // if there isn't a fraction we can just return the whole thing.
+        if fractional_part == 0_u32 {
+            return self;
+        }
+
+        let fractional_mask = fractional_part.overflowing_shr(exponent_clamped).0;
+
+        Self::from_bits(x_bits & !fractional_mask)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::F32;
+
+    #[test]
+    fn sanity_check() {
+        assert_eq!(F32(-1.1).trunc(), F32(-1.0));
+        assert_eq!(F32(-0.1).trunc(), F32(-0.0));
+        assert_eq!(F32(0.0).trunc(), F32(0.0));
+        assert_eq!(F32(1.0).trunc(), F32(1.0));
+        assert_eq!(F32(1.1).trunc(), F32(1.0));
+        assert_eq!(F32(2.9).trunc(), F32(2.0));
+
+        assert_eq!(F32(-100_000_000.13425345345).trunc(), F32(-100_000_000.0));
+        assert_eq!(F32(100_000_000.13425345345).trunc(), F32(100_000_000.0));
+    }
+}
@@ -0,0 +1,104 @@
+//! Embedded-friendly (i.e. `#![no_std]`) math library featuring fast, safe
+//! floating point approximations for common arithmetic operations, as well as
+//! 2D and 3D vector types, statistical analysis functions, and quaternions.
+//!
+//! ## Floating point approximations: `F32` and `F32Ext`
+//!
+//! `micromath` supports approximating many arithmetic operations on `f32`
+//! using bitwise operations, providing great performance and small code size
+//! at the cost of precision. For use cases like graphics and signal
+//! processing, these approximations are often sufficient and the performance
+//! gains worth the lost precision.
+//!
+//! These approximations are defined on the [`F32`] newtype wrapper.
+//!
+//! ### `F32Ext` extension trait
+//!
+//! Floating point approximations can used via the the [`F32Ext`] trait which
+//! is impl'd for `f32`, providing a drop-in `std`-compatible API.
+//!
+//! ```
+//! use micromath::F32Ext;
+//!
+//! let n = 2.0.sqrt();
+//! assert_eq!(n, 1.5); // close enough
+//! ```
+//!
+//! Since the `F32Ext` trait provides methods which are already defined in
+//! `std`, in cases where your crate links `std` the `F32Ext` versions of
+//! the same methods will not be used, in which case you will get an unused
+//! import warning for `F32Ext`.
+//!
+//! If you encounter this, add an `#[allow(unused_imports)]` above the import.
+//!
+//! ```
+//! #[allow(unused_imports)]
+//! use micromath::F32Ext;
+//! ```
+//!
+//! ## Vector types
+//!
+//! See the [`vector`] module for more information on vector types.
+//!
+//! The following vector types are available, all of which have `pub x` and
+//! `pub y` (and on 3D vectors, `pub z`) members:
+//!
+//! | Rust  | 2D      | 3D      |
+//! |-------|---------|---------|
+//! | `i8`  | `I8x2`  | `I8x3`  |
+//! | `i16` | `I16x2` | `I16x3` |
+//! | `i32` | `I32x2` | `I32x3` |
+//! | `u8`  | `U8x2`  | `U8x3`  |
+//! | `u16` | `U16x2` | `U16x3` |
+//! | `u32` | `U32x2` | `U32x3` |
+//! | `f32` | `F32x2` | `F32x3` |
+//!
+//! ## Statistical analysis
+//!
+//! See the [`statistics`] module for more information on statistical analysis
+//! traits and functionality.
+//!
+//! The following traits are available and impl'd for slices and iterators of
+//! `f32` (and can be impl'd for other types):
+//!
+//! - [`Mean`][`statistics::Mean`] - compute arithmetic mean with the `mean()` method.
+//! - [`StdDev`][`statistics::StdDev`] - compute standard deviation with the `stddev()` method.
+//! - [`Trim`][`statistics::Trim`] - cull outliers from a sample slice with the `trim()` method.
+//! - [`Variance`][`statistics::Variance`] - compute variance with the `variance()` method.
+
+#![no_std]
+#![cfg_attr(docsrs, feature(doc_cfg))]
+#![doc(
+    html_logo_url = "https://raw.githubusercontent.com/tarcieri/micromath/main/img/micromath-sq.png",
+    html_root_url = "https://docs.rs/micromath/2.0.0"
+)]
+#![forbid(unsafe_code)]
+#![warn(
+    missing_docs,
+    rust_2018_idioms,
+    trivial_casts,
+    trivial_numeric_casts,
+    unused_qualifications
+)]
+
+#[cfg(feature = "statistics")]
+#[cfg_attr(docsrs, doc(cfg(feature = "statistics")))]
+pub mod statistics;
+
+#[cfg(feature = "vector")]
+#[cfg_attr(docsrs, doc(cfg(feature = "vector")))]
+pub mod vector;
+
+mod f32ext;
+mod float;
+#[cfg(feature = "quaternion")]
+mod quaternion;
+
+pub use crate::{f32ext::F32Ext, float::F32};
+
+#[cfg(feature = "quaternion")]
+pub use crate::quaternion::Quaternion;
+
+#[cfg(feature = "num-traits")]
+#[cfg_attr(docsrs, doc(cfg(feature = "num-traits")))]
+pub use num_traits;
@@ -0,0 +1,448 @@
+//! Adapted from the `madgwick` crate: <https://github.com/japaric/madgwick>
+//! Copyright (c) 2018 Jorge Aparicio
+//!
+//! Original sources dual licensed under your choice of the Apache 2.0
+//! and/or MIT licenses, which matches this crate's licensing terms.
+//!
+//! See toplevel LICENSE-MIT for more information on the MIT license.
+//! Apache 2.0 license follows:
+//!
+//! Licensed under the Apache License, Version 2.0 (the "License");
+//! you may not use this file except in compliance with the License.
+//! You may obtain a copy of the License at:
+//!
+//! <https://www.apache.org/licenses/LICENSE-2.0>
+//!
+//! Unless required by applicable law or agreed to in writing, software
+//! distributed under the License is distributed on an "AS IS" BASIS,
+//! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+//! See the License for the specific language governing permissions and
+//! limitations under the License.
+
+use crate::F32;
+use core::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
+
+#[cfg(feature = "vector")]
+use crate::vector::{Component, F32x3, Vector3d};
+
+/// Quaternions are a number system that extends the complex numbers which can
+/// be used for efficiently computing spatial rotations.
+///
+/// They're computed as the quotient of two directed lines in a
+/// three-dimensional space, or equivalently as the quotient of two vectors.
+///
+/// For given real numbers `a`, `b`, `c`, and `d`, they take the form:
+///
+/// `a + bi + cj + dk`
+///
+/// where `i`, `j`, and `k` are the fundamental quaternion units:
+///
+/// `i² = j² = k² = i*j*k = -1`
+///
+/// Quaternion multiplication is non-commutative:
+///
+/// | x | 1  | i  | j  | k  |
+/// |---|----|----|----|----|
+/// | 1 | 1  | i  | j  | k  |
+/// | i | i  | -1 | k  | -j |
+/// | j | j  | -k | -1 | i  |
+/// | k | k  | j  | -i | -1 |
+#[cfg_attr(docsrs, doc(cfg(feature = "quaternion")))]
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Quaternion(f32, f32, f32, f32);
+
+impl Quaternion {
+    /// Identity quaternion.
+    pub const IDENTITY: Self = Self(1.0, 0.0, 0.0, 0.0);
+
+    /// Create a new quaternion.
+    pub const fn new(a: f32, b: f32, c: f32, d: f32) -> Self {
+        Self(a, b, c, d)
+    }
+
+    /// Get the quaternion that represents the smallest rotation between two vectors.
+    #[cfg(feature = "vector")]
+    pub fn from_two_vectors<C>(u: Vector3d<C>, v: Vector3d<C>) -> Self
+    where
+        C: Component + Into<f32>,
+    {
+        // Implementation from http://lolengine.net/blog/2014/02/24/quaternion-from-two-vectors-final
+        use crate::vector::Vector;
+
+        let n_uv = F32(u.dot(u).into() * v.dot(v).into()).sqrt();
+        let mut realpart = n_uv + u.dot(v).into();
+
+        let w = if realpart < 1e-6 * n_uv {
+            realpart = F32(0.);
+
+            if F32(u.x.into()).abs() > F32(u.z.into()).abs() {
+                Vector3d {
+                    x: -u.y.into(),
+                    y: u.x.into(),
+                    z: 0.,
+                }
+            } else {
+                Vector3d {
+                    x: 0.,
+                    y: -u.z.into(),
+                    z: u.y.into(),
+                }
+            }
+        } else {
+            F32x3 {
+                x: u.x.into(),
+                y: u.y.into(),
+                z: u.z.into(),
+            } * F32x3 {
+                x: v.x.into(),
+                y: v.y.into(),
+                z: v.z.into(),
+            }
+        };
+
+        let q = Quaternion(realpart.0, w.x, w.y, w.z);
+        let n = F32(q.norm()).invsqrt();
+        q * n.0
+    }
+
+    /// Returns the conjugate of this quaternion.
+    pub fn conj(self) -> Self {
+        Quaternion(self.0, -self.1, -self.2, -self.3)
+    }
+
+    /// Returns the dot product of this quaternion.
+    pub fn dot(self, rhs: Self) -> f32 {
+        self.0 * rhs.0 + self.1 * rhs.1 + self.2 * rhs.2 + self.3 * rhs.3
+    }
+
+    /// Compute the inverse of this quaternion.
+    ///
+    /// Panics if [`Quaternion::norm`] is zero.
+    pub fn inv(self) -> Self {
+        let norm = self.norm();
+        assert_ne!(norm, 0.0, "quaternion norm is zero");
+        self.conj() * F32(norm).inv().0
+    }
+
+    /// Compute the magnitude (a.k.a length) of this quaternion.
+    pub fn magnitude(self) -> f32 {
+        F32(self.norm()).sqrt().0
+    }
+
+    /// Returns the norm of this quaternion, i.e. `a²+b²+c²+d²`.
+    ///
+    /// <https://www.mathworks.com/help/aeroblks/quaternionnorm.html>
+    pub fn norm(self) -> f32 {
+        self.0 * self.0 + self.1 * self.1 + self.2 * self.2 + self.3 * self.3
+    }
+
+    /// Compute a quaternion for the given axis vector and angle.
+    #[cfg_attr(docsrs, doc(cfg(feature = "vector")))]
+    pub fn axis_angle<C>(v: Vector3d<C>, theta: C) -> Self
+    where
+        C: Component + Into<f32>,
+    {
+        let half_theta = F32(theta.into() * 0.5);
+
+        // TODO(tarcieri): refactor `Quaternion` to be (f32 + F32x3)
+        let v = F32x3 {
+            x: v.x.into(),
+            y: v.y.into(),
+            z: v.z.into(),
+        } * half_theta.sin().0;
+
+        Self(half_theta.cos().0, v.x, v.y, v.z)
+    }
+
+    /// Rotate a 3D vector using this quaternion, assumes the quaternion is of unit length.
+    #[cfg_attr(docsrs, doc(cfg(feature = "vector")))]
+    pub fn rotate<C>(self, v: Vector3d<C>) -> F32x3
+    where
+        C: Component + Into<f32>,
+    {
+        let Quaternion(qw, qx, qy, qz) = self;
+
+        let qwsq = qw * qw;
+        let qxsq = qx * qx;
+        let qysq = qy * qy;
+        let qzsq = qz * qz;
+
+        let x = (qwsq + qxsq - qysq - qzsq) * v.x.into()
+            + (2. * (qx * qy - qw * qz)) * v.y.into()
+            + (2. * (qx * qz + qw * qy)) * v.z.into();
+
+        let y = (2. * (qx * qy + qw * qz)) * v.x.into()
+            + (qwsq - qxsq + qysq - qzsq) * v.y.into()
+            + (2. * (qy * qz - qw * qx)) * v.z.into();
+
+        let z = (2. * (qx * qz - qw * qy)) * v.x.into()
+            + (2. * (qy * qz + qw * qx)) * v.y.into()
+            + (qwsq - qxsq - qysq + qzsq) * v.z.into();
+
+        F32x3 { x, y, z }
+    }
+
+    /// Scale by a scalar.
+    pub fn scale<S>(self, scalar: S) -> Self
+    where
+        S: Into<f32>,
+    {
+        let k = scalar.into();
+        Self(self.0 * k, self.1 * k, self.2 * k, self.3 * k)
+    }
+
+    /// Normalize the quaternion.
+    pub fn normalize(self) -> Self {
+        let norm = self.norm();
+        assert_ne!(norm, 0.0, "quaternion norm is zero");
+        let n = F32(norm).invsqrt();
+
+        self.scale(n)
+    }
+
+    /// Get the (roll, pitch, yaw) Euler angles, assumes the quaternion is normalized.
+    pub fn to_euler(&self) -> (f32, f32, f32) {
+        let r = F32(2. * (self.0 * self.1 + self.2 * self.3))
+            .atan2(F32(1. - 2. * (self.1 * self.1 + self.2 * self.2)));
+        let p = F32(2. * (self.0 * self.2 - self.1 * self.3)).asin();
+        let y = F32(2. * (self.0 * self.3 + self.1 * self.2))
+            .atan2(F32(1. - 2. * (self.2 * self.2 + self.3 * self.3)));
+
+        (r.0, p.0, y.0)
+    }
+
+    /// Convert this quaternion into an array.
+    pub fn to_array(&self) -> [f32; 4] {
+        [self.0, self.1, self.2, self.3]
+    }
+
+    /// Access the `w` / `a` real component
+    pub fn w(&self) -> f32 {
+        self.0
+    }
+
+    /// Access the `x` / `b` / `i` imaginary component
+    pub fn x(&self) -> f32 {
+        self.1
+    }
+
+    /// Access the `y` / `c` / `j` imaginary component
+    pub fn y(&self) -> f32 {
+        self.2
+    }
+
+    /// Access the `z` / `d` / `k` imaginary component
+    pub fn z(&self) -> f32 {
+        self.3
+    }
+}
+
+impl Add for Quaternion {
+    type Output = Self;
+
+    fn add(self, rhs: Self) -> Self {
+        Self(
+            self.0 + rhs.0,
+            self.1 + rhs.1,
+            self.2 + rhs.2,
+            self.3 + rhs.3,
+        )
+    }
+}
+
+impl AddAssign for Quaternion {
+    fn add_assign(&mut self, rhs: Self) {
+        *self = *self + rhs;
+    }
+}
+
+impl Default for Quaternion {
+    fn default() -> Self {
+        Self::IDENTITY
+    }
+}
+
+impl From<(f32, f32, f32, f32)> for Quaternion {
+    fn from(q: (f32, f32, f32, f32)) -> Quaternion {
+        Self::new(q.0, q.1, q.2, q.3)
+    }
+}
+
+impl From<Quaternion> for (f32, f32, f32, f32) {
+    fn from(q: Quaternion) -> (f32, f32, f32, f32) {
+        (q.0, q.1, q.2, q.3)
+    }
+}
+
+impl Mul for Quaternion {
+    type Output = Self;
+
+    fn mul(self, other: Self) -> Self {
+        Self(
+            self.0 * other.0 - self.1 * other.1 - self.2 * other.2 - self.3 * other.3,
+            self.0 * other.1 + self.1 * other.0 + self.2 * other.3 - self.3 * other.2,
+            self.0 * other.2 - self.1 * other.3 + self.2 * other.0 + self.3 * other.1,
+            self.0 * other.3 + self.1 * other.2 - self.2 * other.1 + self.3 * other.0,
+        )
+    }
+}
+
+impl Mul<f32> for Quaternion {
+    type Output = Self;
+
+    fn mul(self, k: f32) -> Self {
+        self.scale(k)
+    }
+}
+
+impl Mul<Quaternion> for f32 {
+    type Output = Quaternion;
+
+    fn mul(self, q: Quaternion) -> Quaternion {
+        q.scale(self)
+    }
+}
+
+impl MulAssign<f32> for Quaternion {
+    fn mul_assign(&mut self, k: f32) {
+        *self = *self * k;
+    }
+}
+
+impl Sub for Quaternion {
+    type Output = Self;
+
+    fn sub(self, rhs: Self) -> Self {
+        Self(
+            self.0 - rhs.0,
+            self.1 - rhs.1,
+            self.2 - rhs.2,
+            self.3 - rhs.3,
+        )
+    }
+}
+
+impl SubAssign for Quaternion {
+    fn sub_assign(&mut self, rhs: Self) {
+        *self = *self - rhs;
+    }
+}
+
+#[cfg(feature = "vector")]
+#[cfg_attr(docsrs, doc(cfg(feature = "vector")))]
+impl<C> From<Vector3d<C>> for Quaternion
+where
+    C: Component + Into<f32>,
+{
+    fn from(v: Vector3d<C>) -> Quaternion {
+        Self(0.0, v.x.into(), v.y.into(), v.z.into())
+    }
+}
+
+#[cfg(feature = "vector")]
+#[cfg_attr(docsrs, doc(cfg(feature = "vector")))]
+impl<C> Mul<Vector3d<C>> for Quaternion
+where
+    C: Component + Into<f32>,
+{
+    type Output = F32x3;
+
+    fn mul(self, v: Vector3d<C>) -> F32x3 {
+        self.rotate(v)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::Quaternion;
+    use crate::vector::Vector3d;
+    use crate::F32Ext;
+
+    const MAX_ERROR: f32 = 0.05;
+
+    #[test]
+    fn conj_test() {
+        let q = Quaternion(1.0, 2.0, 3.0, 4.0);
+        assert_eq!(q.conj(), Quaternion(1.0, -2.0, -3.0, -4.0));
+    }
+
+    #[test]
+    fn norm_test() {
+        let q = Quaternion(1.0, 2.0, 3.0, 4.0);
+        assert_eq!(q.norm(), 30.0);
+
+        let n = q.norm().invsqrt();
+        let r = q * n;
+
+        // The magnitude of the norm should be 1.0
+        let allowed_delta = 1.0 * MAX_ERROR;
+        let actual_delta = (r.norm() - 1.0).abs();
+
+        assert!(
+            actual_delta <= allowed_delta,
+            "delta {} too large: {} vs {}",
+            actual_delta,
+            r.norm(),
+            1.0
+        );
+    }
+
+    #[test]
+    fn add_assign() {
+        let mut q = Quaternion(1.0, 2.0, 3.0, 4.0);
+        q += Quaternion(1.0, 2.0, 3.0, 4.0);
+        assert_eq!(q, Quaternion(2.0, 4.0, 6.0, 8.0));
+    }
+
+    #[test]
+    fn mul_assign() {
+        let mut q = Quaternion(1.0, 2.0, 3.0, 4.0);
+        q *= 2.0;
+        assert_eq!(q, Quaternion(2.0, 4.0, 6.0, 8.0));
+    }
+
+    #[test]
+    fn sub_assign() {
+        let mut q = Quaternion(2.0, 4.0, 6.0, 8.0);
+        q -= Quaternion(1.0, 2.0, 3.0, 4.0);
+        assert_eq!(q, Quaternion(1.0, 2.0, 3.0, 4.0));
+    }
+
+    #[test]
+    fn mul_quaternion() {
+        let q = Quaternion(1.0, 2.0, 3.0, 4.0);
+        let r = Quaternion(4.0, 3.0, 2.0, 1.0);
+        assert_eq!(q * r, Quaternion(-12.0, 6.0, 24.0, 12.0));
+    }
+
+    #[test]
+    fn mul_f32() {
+        let q = Quaternion(1.0, 2.0, 3.0, 4.0);
+        let r = 2.0 * q;
+        assert_eq!(r, Quaternion(2.0, 4.0, 6.0, 8.0));
+
+        let s = r * 0.5;
+        assert_eq!(s, Quaternion(1.0, 2.0, 3.0, 4.0));
+    }
+
+    #[test]
+    fn smallest_rot() {
+        use crate::vector::Vector;
+
+        let v1 = Vector3d {
+            x: 0.,
+            y: 0.,
+            z: 1.,
+        };
+        let v2 = Vector3d {
+            x: 0.,
+            y: 1.,
+            z: 0.,
+        };
+
+        let q = Quaternion::from_two_vectors(v1, v2);
+        let v1_r = q.rotate(v1);
+
+        assert!((v1_r - v2).magnitude() < 1e-1);
+    }
+}
@@ -0,0 +1,23 @@
+//! Statistical analysis support.
+//!
+//! The `statistics` Cargo feature must be enabled to use this functionality.
+//!
+//! The following traits are available and impl'd for slices and iterators of
+//! `f32` (and can be impl'd for other types):
+//!
+//! - [Mean] - compute arithmetic mean with the `mean()` method.
+//! - [StdDev] - compute standard deviation with the `stddev()` method
+//! - [Trim] - cull outliers from a sample slice with the `trim()` method.
+//! - [Variance] - compute variance with the `variance() method.
+//!
+//! [Mean]: https://docs.rs/micromath/latest/micromath/statistics/trait.Mean.html
+//! [StdDev]: https://docs.rs/micromath/latest/micromath/statistics/trait.StdDev.html
+//! [Trim]: https://docs.rs/micromath/latest/micromath/statistics/trim/trait.Trim.html
+//! [Variance]: https://docs.rs/micromath/latest/micromath/statistics/trait.Variance.html
+
+mod mean;
+mod stddev;
+pub mod trim;
+mod variance;
+
+pub use self::{mean::Mean, stddev::StdDev, trim::Trim, variance::Variance};
@@ -0,0 +1,37 @@
+/// Arithmetic mean
+pub trait Mean {
+    /// Result type
+    type Result;
+
+    /// Compute the arithmetic mean
+    fn mean(self) -> Self::Result;
+}
+
+impl<I> Mean for I
+where
+    I: Iterator<Item = f32>,
+{
+    type Result = f32;
+
+    fn mean(self) -> f32 {
+        let mut num_items = 0;
+        let mut sum = 0.0;
+
+        for item in self {
+            num_items += 1;
+            sum += item;
+        }
+
+        sum / (num_items as f32)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::Mean;
+
+    #[test]
+    fn mean_test() {
+        assert_eq!([1.0, 3.0, 5.0].iter().cloned().mean(), 3.0);
+    }
+}
@@ -0,0 +1,35 @@
+use super::variance::Variance;
+#[allow(unused_imports)]
+use crate::F32Ext;
+
+/// Compute standard deviation
+pub trait StdDev {
+    /// Result type
+    type Result;
+
+    /// Compute standard deviation
+    fn stddev(self) -> Self::Result;
+}
+
+impl<T> StdDev for &[T]
+where
+    T: Copy + Into<f32>,
+{
+    /// Result type
+    type Result = f32;
+
+    fn stddev(self) -> f32 {
+        self.variance().sqrt()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::StdDev;
+
+    #[test]
+    fn stddev_test() {
+        assert_eq!([1.0, 3.0, 5.0].as_ref().stddev(), 2.0);
+        assert_eq!([1.0, 3.0, 5.0].as_ref().stddev(), 2.0);
+    }
+}
@@ -0,0 +1,123 @@
+//! Iterate over input slices after culling statistical outliers.
+
+use super::mean::Mean;
+#[allow(unused_imports)]
+use crate::F32Ext;
+use core::{iter, slice};
+
+/// Default number of standard deviations at which a value should be
+/// considered an outlier.
+pub const DEFAULT_THRESHOLD: f32 = 2.0;
+
+/// Iterate over the given input after culling outliers.
+pub trait Trim {
+    /// Result type
+    type Result: Copy + Into<f32>;
+
+    /// Trim this collection (cull outliers) at the default threshold of 2 standard
+    /// deviations.
+    fn trim(&self) -> Iter<'_, Self::Result> {
+        self.trim_at(DEFAULT_THRESHOLD)
+    }
+
+    /// Trim this collection (cull outliers) at the specified number of standard deviations.
+    fn trim_at(&self, threshold: f32) -> Iter<'_, Self::Result>;
+}
+
+impl<N> Trim for &[N]
+where
+    N: Copy,
+    f32: From<N>,
+{
+    type Result = N;
+
+    fn trim_at(&self, threshold: f32) -> Iter<'_, Self::Result> {
+        Iter::new(self, threshold)
+    }
+}
+
+/// A "trimmed" iterator which culls outliers at a given number of standard
+/// deviations from the mean.
+pub struct Iter<'a, N: Copy> {
+    /// Iterator over the input values
+    input: iter::Cloned<slice::Iter<'a, N>>,
+
+    /// Arithmetic mean of the input values as a float
+    mean: f32,
+
+    /// Standard deviation of the input values
+    stddev: f32,
+
+    /// Number of standard deviations at which values are considered outliers
+    threshold: f32,
+}
+
+impl<'a, N> Iter<'a, N>
+where
+    N: Copy,
+    f32: From<N>,
+{
+    /// Create a new trimmed iterator over an input slice.
+    ///
+    /// Inputs will be considered outliers at the given `threshold` number
+    /// of standard deviations (e.g. `2.0`).
+    pub fn new(input: &'a [N], threshold: f32) -> Self {
+        let len = input.len() as f32;
+        let input = input.iter().cloned();
+        let input_f32 = input.clone().map(f32::from);
+
+        // TODO(tarcieri): eliminate duplication with mean/variance/stddev in super
+        let mean = input_f32.clone().mean();
+        let sum = input_f32.fold(0.0, |sum, n| {
+            let n = n - mean;
+            sum + n * n
+        });
+        let variance = sum / (len - 1.0);
+        let stddev = variance.sqrt();
+
+        Self {
+            input,
+            mean,
+            stddev,
+            threshold,
+        }
+    }
+}
+
+impl<'a, N> Iterator for Iter<'a, N>
+where
+    N: Copy,
+    f32: From<N>,
+{
+    type Item = N;
+
+    fn next(&mut self) -> Option<N> {
+        while let Some(n) = self.input.next() {
+            let distance = (f32::from(n) - self.mean).abs();
+
+            // TODO(tarcieri): better method for finding outliers? (e.g. MAD, IQD)
+            if (distance / self.stddev) < self.threshold {
+                return Some(n);
+            }
+        }
+
+        None
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::Trim;
+
+    #[test]
+    fn stddev_test() {
+        let input: &[f32] = &[1.0, 2.0, 3.0, 999.0, 5.0];
+        let mut trimmed = input.trim_at(1.0);
+
+        assert_eq!(trimmed.next().unwrap(), 1.0);
+        assert_eq!(trimmed.next().unwrap(), 2.0);
+        assert_eq!(trimmed.next().unwrap(), 3.0);
+        assert_eq!(trimmed.next().unwrap(), 5.0);
+        assert_eq!(trimmed.next(), None);
+    }
+}
@@ -0,0 +1,40 @@
+use super::mean::Mean;
+
+/// Statistical variance
+pub trait Variance {
+    /// Result type
+    type Result;
+
+    /// Compute the statistical variance
+    fn variance(self) -> Self::Result;
+}
+
+impl<T> Variance for &[T]
+where
+    T: Copy + Into<f32>,
+{
+    /// Result type
+    type Result = f32;
+
+    fn variance(self) -> f32 {
+        let mean = self.iter().map(|n| (*n).into()).mean();
+        let mut sum = 0.0;
+
+        for item in self {
+            let n = (*item).into() - mean;
+            sum += n * n;
+        }
+
+        sum / (self.len() as f32 - 1.0)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::Variance;
+
+    #[test]
+    fn variance_test() {
+        assert_eq!([1.0, 3.0, 5.0].as_ref().variance(), 4.0);
+    }
+}
@@ -0,0 +1,75 @@
+//! Algebraic vector types generic over a component type.
+//!
+//! All vectors types impl the [Vector] trait, and all vector components
+//! impl the [Component] trait.
+
+mod component;
+mod iter;
+mod vector2d;
+mod vector3d;
+
+pub use self::{
+    component::Component,
+    iter::Iter,
+    vector2d::{F32x2, I16x2, I32x2, I8x2, U16x2, U32x2, U8x2, Vector2d},
+    vector3d::{F32x3, I16x3, I32x3, I8x3, U16x3, U32x3, U8x3, Vector3d},
+};
+
+use core::{fmt::Debug, iter::FromIterator};
+
+#[allow(unused_imports)]
+use crate::F32Ext;
+
+/// Algebraic vector generic over a given [`Component`] type.
+pub trait Vector<C>: Copy + Debug + Default + FromIterator<C> + Send + Sync
+where
+    C: Component,
+{
+    /// Number of axes
+    const AXES: usize;
+
+    /// Get the component value for a particular index
+    fn get(self, index: usize) -> Option<C>;
+
+    /// Compute the dot product of two vectors
+    fn dot(self, rhs: Self) -> C;
+
+    /// Instantiate a vector from a slice of components.
+    ///
+    /// Panics if the slice is not the right size.
+    fn from_slice(slice: &[C]) -> Self {
+        Self::from_iter(slice.iter().cloned())
+    }
+
+    /// Iterate over the components of this vector
+    fn iter(&self) -> Iter<'_, Self, C> {
+        Iter::new(self)
+    }
+
+    /// Compute the distance between two vectors
+    fn distance(self, rhs: Self) -> f32
+    where
+        C: Into<f32>,
+    {
+        let differences = self
+            .iter()
+            .zip(rhs.iter())
+            .map(|(a, b)| a.into() - b.into());
+
+        differences.map(|n| n * n).sum::<f32>().sqrt()
+    }
+
+    /// Compute the magnitude of a vector
+    fn magnitude(self) -> f32
+    where
+        C: Into<f32>,
+    {
+        self.iter()
+            .map(|n| {
+                let n = n.into();
+                n * n
+            })
+            .sum::<f32>()
+            .sqrt()
+    }
+}
@@ -0,0 +1,43 @@
+//! Components of numeric vectors.
+
+use crate::F32;
+use core::{
+    fmt::Debug,
+    ops::{Add, Div, Mul, Sub},
+};
+
+/// Components of numeric vectors.
+///
+/// All components must be [`Copy`] + [`Sized`] types which support a minimal
+/// set of arithmetic operations ([`Add`], [`Sub`], [`Mul`], [`Div`]), as well as
+/// [`Default`], [`PartialEq`] and [`PartialOrd`].
+///
+/// This trait is impl'd for the following primitive types:
+///
+/// - [`i8`], [`i16`], [`i32`]
+/// - [`u8`], [`u16`], [`u32`]
+/// - [`f32`]
+pub trait Component:
+    Copy
+    + Debug
+    + Default
+    + PartialEq
+    + PartialOrd
+    + Send
+    + Sized
+    + Sync
+    + Add<Output = Self>
+    + Sub<Output = Self>
+    + Mul<Output = Self>
+    + Div<Output = Self>
+{
+}
+
+impl Component for i8 {}
+impl Component for i16 {}
+impl Component for i32 {}
+impl Component for u8 {}
+impl Component for u16 {}
+impl Component for u32 {}
+impl Component for f32 {}
+impl Component for F32 {}
@@ -0,0 +1,54 @@
+//! Iterator over the components of an algebraic vector
+
+use super::{Component, Vector};
+use core::marker::PhantomData;
+
+/// Iterator over the components of an algebraic vector
+#[derive(Clone, Debug)]
+pub struct Iter<'a, V, C>
+where
+    V: Vector<C>,
+    C: Component,
+{
+    /// Reference to the original vector
+    vector: &'a V,
+
+    /// Iteration position within the vector
+    position: usize,
+
+    /// Component type
+    component: PhantomData<C>,
+}
+
+impl<'a, V, C> Iter<'a, V, C>
+where
+    V: Vector<C>,
+    C: Component,
+{
+    /// Create a new iterator over the vector's components
+    pub(super) fn new(vector: &'a V) -> Self {
+        Self {
+            vector,
+            position: 0,
+            component: PhantomData,
+        }
+    }
+}
+
+impl<'a, V, C> Iterator for Iter<'a, V, C>
+where
+    V: Vector<C>,
+    C: Component,
+{
+    type Item = C;
+
+    fn next(&mut self) -> Option<C> {
+        let item = self.vector.get(self.position);
+
+        if item.is_some() {
+            self.position += 1;
+        }
+
+        item
+    }
+}
@@ -0,0 +1,222 @@
+//! 2-dimensional vector
+
+use super::{Component, Vector};
+use core::{
+    iter::FromIterator,
+    ops::{Add, AddAssign, Index, Mul, MulAssign, Sub, SubAssign},
+};
+
+/// 2-dimensional XY vector of `i8` values
+pub type I8x2 = Vector2d<i8>;
+
+/// 2-dimensional XY vector of `i16` values
+pub type I16x2 = Vector2d<i16>;
+
+/// 2-dimensional XY vector of `i32` values
+pub type I32x2 = Vector2d<i32>;
+
+/// 2-dimensional XY vector of `u8` values
+pub type U8x2 = Vector2d<u8>;
+
+/// 2-dimensional XY vector of `u16` values
+pub type U16x2 = Vector2d<u16>;
+
+/// 2-dimensional XY vector of `u32` values
+pub type U32x2 = Vector2d<u32>;
+
+/// 2-dimensional XY vector of `f32` values
+pub type F32x2 = Vector2d<f32>;
+
+/// 2-dimensional vector
+#[derive(Copy, Clone, Debug, Default, PartialEq)]
+pub struct Vector2d<C: Component> {
+    /// X component
+    pub x: C,
+
+    /// Y component
+    pub y: C,
+}
+
+impl<C> Vector2d<C>
+where
+    C: Component,
+{
+    /// Return a 2-element array containing the coordinates
+    // TODO(tarcieri): move this to the `Vector` trait leveraging const generics?
+    pub fn to_array(&self) -> [C; 2] {
+        [self.x, self.y]
+    }
+}
+
+impl<C> FromIterator<C> for Vector2d<C>
+where
+    C: Component,
+{
+    fn from_iter<T>(into_iter: T) -> Self
+    where
+        T: IntoIterator<Item = C>,
+    {
+        let mut iter = into_iter.into_iter();
+
+        let x = iter.next().expect("no x-axis component in slice");
+        let y = iter.next().expect("no y-axis component in slice");
+
+        assert!(
+            iter.next().is_none(),
+            "too many items for 2-dimensional vector"
+        );
+
+        Self { x, y }
+    }
+}
+
+impl<C> Vector<C> for Vector2d<C>
+where
+    C: Component,
+{
+    const AXES: usize = 2;
+
+    fn get(self, index: usize) -> Option<C> {
+        match index {
+            0 => Some(self.x),
+            1 => Some(self.y),
+            _ => None,
+        }
+    }
+
+    fn dot(self, rhs: Self) -> C {
+        (self.x * rhs.x) + (self.y * rhs.y)
+    }
+}
+
+impl<C> From<(C, C)> for Vector2d<C>
+where
+    C: Component,
+{
+    fn from(vector: (C, C)) -> Self {
+        Self {
+            x: vector.0,
+            y: vector.1,
+        }
+    }
+}
+
+impl<C> Index<usize> for Vector2d<C>
+where
+    C: Component,
+{
+    type Output = C;
+
+    fn index(&self, i: usize) -> &C {
+        match i {
+            0 => &self.x,
+            1 => &self.y,
+            _ => panic!("index out of range"),
+        }
+    }
+}
+
+impl<C> Add for Vector2d<C>
+where
+    C: Component,
+{
+    type Output = Self;
+
+    fn add(self, rhs: Self) -> Self {
+        Self {
+            x: self.x + rhs.x,
+            y: self.y + rhs.y,
+        }
+    }
+}
+
+impl<C> AddAssign for Vector2d<C>
+where
+    C: Component,
+{
+    fn add_assign(&mut self, other: Self) {
+        *self = *self + other;
+    }
+}
+
+impl<C> Sub for Vector2d<C>
+where
+    C: Component,
+{
+    type Output = Self;
+
+    fn sub(self, rhs: Self) -> Self {
+        Self {
+            x: self.x - rhs.x,
+            y: self.y - rhs.y,
+        }
+    }
+}
+
+impl<C> SubAssign for Vector2d<C>
+where
+    C: Component,
+{
+    fn sub_assign(&mut self, other: Self) {
+        *self = *self - other;
+    }
+}
+
+impl<C> Mul<C> for Vector2d<C>
+where
+    C: Component,
+{
+    type Output = Self;
+
+    fn mul(self, rhs: C) -> Self {
+        Self {
+            x: self.x * rhs,
+            y: self.y * rhs,
+        }
+    }
+}
+
+impl<C> MulAssign<C> for Vector2d<C>
+where
+    C: Component,
+{
+    fn mul_assign(&mut self, rhs: C) {
+        *self = *self * rhs;
+    }
+}
+
+impl From<I8x2> for F32x2 {
+    fn from(vector: I8x2) -> F32x2 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+        }
+    }
+}
+
+impl From<I16x2> for F32x2 {
+    fn from(vector: I16x2) -> F32x2 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+        }
+    }
+}
+
+impl From<U8x2> for F32x2 {
+    fn from(vector: U8x2) -> F32x2 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+        }
+    }
+}
+
+impl From<U16x2> for F32x2 {
+    fn from(vector: U16x2) -> F32x2 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+        }
+    }
+}
@@ -0,0 +1,264 @@
+//! 3-dimensional vector
+
+use super::{Component, Vector};
+use crate::F32;
+use core::{
+    iter::FromIterator,
+    ops::{Add, AddAssign, Index, Mul, MulAssign, Sub, SubAssign},
+};
+
+/// 3-dimensional XYZ vector of `i8` values
+pub type I8x3 = Vector3d<i8>;
+
+/// 3-dimensional XYZ vector of `i16` values
+pub type I16x3 = Vector3d<i16>;
+
+/// 3-dimensional XYZ vector of `i32` values
+pub type I32x3 = Vector3d<i32>;
+
+/// 3-dimensional XYZ vector of `u8` values
+pub type U8x3 = Vector3d<u8>;
+
+/// 3-dimensional XYZ vector of `u16` values
+pub type U16x3 = Vector3d<u16>;
+
+/// 3-dimensional XYZ vector of `u32` values
+pub type U32x3 = Vector3d<u32>;
+
+/// 3-dimensional XYZ vector of `f32` values
+pub type F32x3 = Vector3d<f32>;
+
+/// 3-dimensional vector
+#[derive(Copy, Clone, Debug, Default, PartialEq)]
+pub struct Vector3d<C: Component> {
+    /// X component
+    pub x: C,
+
+    /// Y component
+    pub y: C,
+
+    /// Z component
+    pub z: C,
+}
+
+impl<C> Vector3d<C>
+where
+    C: Component,
+{
+    /// Return a 3-element array containing the coordinates
+    // TODO(tarcieri): move this to the `Vector` trait leveraging const generics?
+    pub fn to_array(&self) -> [C; 3] {
+        [self.x, self.y, self.z]
+    }
+}
+
+impl<C> FromIterator<C> for Vector3d<C>
+where
+    C: Component,
+{
+    fn from_iter<T>(into_iter: T) -> Self
+    where
+        T: IntoIterator<Item = C>,
+    {
+        let mut iter = into_iter.into_iter();
+
+        let x = iter.next().expect("no x-axis component in slice");
+        let y = iter.next().expect("no y-axis component in slice");
+        let z = iter.next().expect("no z-axis component in slice");
+
+        assert!(
+            iter.next().is_none(),
+            "too many items for 3-dimensional vector"
+        );
+
+        Self { x, y, z }
+    }
+}
+
+impl<C> Vector<C> for Vector3d<C>
+where
+    C: Component,
+{
+    const AXES: usize = 3;
+
+    fn get(self, index: usize) -> Option<C> {
+        match index {
+            0 => Some(self.x),
+            1 => Some(self.y),
+            2 => Some(self.z),
+            _ => None,
+        }
+    }
+
+    fn dot(self, rhs: Self) -> C {
+        (self.x * rhs.x) + (self.y * rhs.y) + (self.z * rhs.z)
+    }
+}
+
+impl<C> From<(C, C, C)> for Vector3d<C>
+where
+    C: Component,
+{
+    fn from(vector: (C, C, C)) -> Self {
+        Self {
+            x: vector.0,
+            y: vector.1,
+            z: vector.2,
+        }
+    }
+}
+
+impl<C> Index<usize> for Vector3d<C>
+where
+    C: Component,
+{
+    type Output = C;
+
+    fn index(&self, i: usize) -> &C {
+        match i {
+            0 => &self.x,
+            1 => &self.y,
+            2 => &self.z,
+            _ => panic!("index out of range"),
+        }
+    }
+}
+
+impl<C> Add for Vector3d<C>
+where
+    C: Component,
+{
+    type Output = Self;
+
+    fn add(self, rhs: Self) -> Self {
+        Self {
+            x: self.x + rhs.x,
+            y: self.y + rhs.y,
+            z: self.z + rhs.z,
+        }
+    }
+}
+
+impl<C> AddAssign for Vector3d<C>
+where
+    C: Component,
+{
+    fn add_assign(&mut self, other: Self) {
+        *self = *self + other;
+    }
+}
+
+impl<C> Sub for Vector3d<C>
+where
+    C: Component,
+{
+    type Output = Self;
+
+    fn sub(self, rhs: Self) -> Self {
+        Self {
+            x: self.x - rhs.x,
+            y: self.y - rhs.y,
+            z: self.z - rhs.z,
+        }
+    }
+}
+
+impl<C> SubAssign for Vector3d<C>
+where
+    C: Component,
+{
+    fn sub_assign(&mut self, other: Self) {
+        *self = *self - other;
+    }
+}
+
+/// Compute the cross product of two vectors
+impl<C> Mul for Vector3d<C>
+where
+    C: Component,
+{
+    type Output = Self;
+
+    fn mul(self, rhs: Self) -> Self {
+        Self {
+            x: (self.y * rhs.z) - (self.z * rhs.y),
+            y: (self.z * rhs.x) - (self.x * rhs.z),
+            z: (self.x * rhs.y) - (self.y * rhs.x),
+        }
+    }
+}
+
+/// Multiply vector by a scalar component
+impl<C> Mul<C> for Vector3d<C>
+where
+    C: Component,
+{
+    type Output = Self;
+
+    fn mul(self, rhs: C) -> Self {
+        Self {
+            x: self.x * rhs,
+            y: self.y * rhs,
+            z: self.z * rhs,
+        }
+    }
+}
+
+impl<C> MulAssign<C> for Vector3d<C>
+where
+    C: Component,
+{
+    fn mul_assign(&mut self, rhs: C) {
+        *self = *self * rhs;
+    }
+}
+
+impl From<I8x3> for F32x3 {
+    fn from(vector: I8x3) -> F32x3 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+            z: vector.z.into(),
+        }
+    }
+}
+
+impl From<I16x3> for F32x3 {
+    fn from(vector: I16x3) -> F32x3 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+            z: vector.z.into(),
+        }
+    }
+}
+
+impl From<U8x3> for F32x3 {
+    fn from(vector: U8x3) -> F32x3 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+            z: vector.z.into(),
+        }
+    }
+}
+
+impl From<U16x3> for F32x3 {
+    fn from(vector: U16x3) -> F32x3 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+            z: vector.z.into(),
+        }
+    }
+}
+
+impl From<Vector3d<F32>> for F32x3 {
+    fn from(vector: Vector3d<F32>) -> F32x3 {
+        Self {
+            x: vector.x.into(),
+            y: vector.y.into(),
+            z: vector.z.into(),
+        }
+    }
+}