alice-det-math 0.5.0

//! `f32` kernels: Cephes `sinf` / `cosf` / `expf`, musl `logf`, plus
//! `powf` / `powi` / `cbrt` / `hypot`, and the `f32` entry points of the
//! double-precision fdlibm kernels in [`crate::double`].

use crate::double::{
    acos64, asin64, exp64, k_cos64, k_sin64, ln64, log10_64, log2_64, reduce_pio2_64,
};
use crate::ops::sqrt;

/// π/2 split into three single-precision pieces (Cephes `DP1..3` × 2) so
/// `x - k·π/2` is computed with ~24 extra bits for `|k| < 2^24`.
/// First piece of π/2 (Cephes `DP1`×2).
pub const PIO2_1: f32 = 1.570_312_5;
/// Second piece of π/2.
pub const PIO2_2: f32 = 4.837_512_969_970_703_125e-4;
/// Third piece of π/2.
pub const PIO2_3: f32 = 7.549_789_948_768_648e-8;
/// Fourth piece of π/2: `π/2 − (PIO2_1 + PIO2_2 + PIO2_3)` rounded to `f32`.
///
/// With it, and with `k·PIO2_3` formed exactly (Dekker two-product, halves
/// split with [`SPLIT_MASK`]), the reduced argument near a zero of `sin` /
/// `cos` keeps its relative accuracy out to `|x| ≤ 8192` (the three-piece
/// reduction lost up to 14 ulp at `|x| ≤ 100`).
pub const PIO2_4: f32 = -1.715_124_510_005_881_9e-15;
/// Bits kept by the high half of a Dekker split (12 of the 24 significand bits).
pub const SPLIT_MASK: u32 = 0xFFFF_F000;
/// `2/π`, the range-reduction multiplier.
pub const FRAC_2_PI: f32 = core::f32::consts::FRAC_2_PI;

/// Cephes `sinf` polynomial coefficients on `[-π/4, π/4]`.
pub const SIN_P: [f32; 3] = [-1.951_529_589_1e-4, 8.332_160_873_6e-3, -1.666_665_461_1e-1];
/// Cephes `cosf` polynomial coefficients on `[-π/4, π/4]`.
pub const COS_P: [f32; 3] = [
    2.443_315_711_809_948e-5,
    -1.388_731_625_493_765e-3,
    4.166_664_568_298_827e-2,
];

/// Sine polynomial on `[-π/4, π/4]` (Cephes `sinf`).
#[inline(always)]
fn sin_poly(r: f32) -> f32 {
    let z = r * r;
    let p = ((SIN_P[0] * z + SIN_P[1]) * z + SIN_P[2]) * z * r;
    r + p
}

/// Cosine polynomial on `[-π/4, π/4]` (Cephes `cosf`).
#[inline(always)]
fn cos_poly(r: f32) -> f32 {
    let z = r * r;
    let p = ((COS_P[0] * z + COS_P[1]) * z + COS_P[2]) * z * z;
    1.0 - 0.5 * z + p
}

/// Nearest integer of `t` for range reduction: `trunc(t + copysign(0.5, t))`.
///
/// Basic operations only (deterministic), 3 ops instead of the ~20 of the
/// exact `round`; it differs from `round` only when `t` is within an ulp of
/// an exact tie, where either neighbour is a valid reduction (`|r|` lands on
/// the `π/4` boundary the polynomials are fitted to).
#[inline(always)]
pub(crate) fn nearest_int(t: f32) -> f32 {
    // `t < 0` (not the sign bit: `-0.0` and NaN take +0.5, which changes
    // nothing) — the SIMD version is the same compare + blend
    let half = if t < 0.0 { -0.5 } else { 0.5 };
    // `as i32` saturates and maps NaN to 0, like the SIMD `trunc_int`
    ((t + half) as i32) as f32
}

/// High 12 significand bits of `a` (exact: the low bits are cleared, never
/// rounded, so the split cannot overflow the way `4097·a` does).
#[inline(always)]
pub(crate) fn split_hi(a: f32) -> f32 {
    f32::from_bits(a.to_bits() & SPLIT_MASK)
}

/// Dekker's exact product: `a·b = p + e` with `p = fl(a·b)`, using a split of
/// both factors into 12-bit halves (every partial product is exact in `f32`)
/// and no fused multiply-add.
#[inline(always)]
pub(crate) fn two_prod(a: f32, b: f32) -> (f32, f32) {
    let p = a * b;
    let (ah, bh) = (split_hi(a), split_hi(b));
    let (al, bl) = (a - ah, b - bh);
    (p, (((ah * bh - p) + ah * bl) + al * bh) + al * bl)
}

/// Reduce `x` to `(k mod 4, r)` with `x = k·π/2 + r`, `|r| ≤ π/4`.
///
/// `k·PIO2_1` and `k·PIO2_2` are exact for `|k| < 2^13`; `k·PIO2_3` is formed
/// exactly as `p3 + e3`, and `PIO2_4` carries the rest of π/2.
#[inline(always)]
fn reduce_pio2(x: f32) -> (i32, f32) {
    let kf = nearest_int(x * FRAC_2_PI);
    let (p3, e3) = two_prod(kf, PIO2_3);
    let r = ((((x - kf * PIO2_1) - kf * PIO2_2) - p3) - e3) - kf * PIO2_4;
    // `kf` is integral and |kf| < 2^31 for every finite f32 we accept; the
    // cast saturates for huge inputs, which only affects the (already
    // meaningless) quadrant of astronomically large arguments.
    let k = (kf as i32) & 3;
    (k, r)
}

/// A NaN produced *by arithmetic* (`inf - inf` in the range reduction of an
/// astronomically large argument) carries the platform's default NaN, whose
/// sign bit differs between `x86` (negative) and `AArch64` (positive). Replace
/// it with the canonical NaN; finite results pass through untouched.
#[inline(always)]
pub(crate) fn canon(v: f32) -> f32 {
    if v.is_nan() {
        f32::NAN
    } else {
        v
    }
}

/// Deterministic `sin(x)`.
///
/// `|x| > 2^24` is beyond the single-precision reduction: the result is
/// deterministic but meaningless (`±inf` or the canonical NaN), as with
/// Cephes `sinf` (which returns 0 there) — reduce the argument first.
#[inline]
#[must_use]
pub fn sin(x: f32) -> f32 {
    if !x.is_finite() {
        return f32::NAN;
    }
    if x == 0.0 {
        return x; // ±0 → ±0 (the polynomial's `r + p` would turn −0 into +0)
    }
    let (k, r) = reduce_pio2(x);
    canon(match k {
        0 => sin_poly(r),
        1 => cos_poly(r),
        2 => -sin_poly(r),
        _ => -cos_poly(r),
    })
}

/// Deterministic `cos(x)` (see [`sin`] for the `|x| > 2^24` caveat).
#[inline]
#[must_use]
pub fn cos(x: f32) -> f32 {
    if !x.is_finite() {
        return f32::NAN;
    }
    let (k, r) = reduce_pio2(x);
    canon(match k {
        0 => cos_poly(r),
        1 => -sin_poly(r),
        2 => -cos_poly(r),
        _ => sin_poly(r),
    })
}

/// Deterministic `(sin(x), cos(x))` from one range reduction; bit-identical
/// to calling [`sin`] and [`cos`] separately.
#[inline]
#[must_use]
pub fn sin_cos(x: f32) -> (f32, f32) {
    if !x.is_finite() {
        return (f32::NAN, f32::NAN);
    }
    if x == 0.0 {
        return (x, 1.0); // as `sin`: ±0 → ±0
    }
    let (k, r) = reduce_pio2(x);
    let s = sin_poly(r);
    let c = cos_poly(r);
    let (s, c) = match k {
        0 => (s, c),
        1 => (c, -s),
        2 => (-s, -c),
        _ => (-c, s),
    };
    (canon(s), canon(c))
}

/// `log2(e)`, the `exp` range-reduction multiplier.
pub const LOG2E: f32 = core::f32::consts::LOG2_E;
/// `ln 2` split: high part (exact in 12 bits).
pub const LN2_HI: f32 = 0.693_145_751_953_125;
/// `ln 2` split: low part.
pub const LN2_LO: f32 = 1.428_606_765_330_187e-6;
/// Cephes `expf` polynomial for `e^r`, `|r| ≤ ln2/2` (highest degree first).
pub const EXP_P: [f32; 5] = [
    1.987_569_150_0e-4,
    1.398_199_950_7e-3,
    8.333_451_907_3e-3,
    4.166_579_589_4e-2,
    1.666_666_545_9e-1,
];
/// `exp` overflows to `+inf` above this and underflows to `0` below [`EXP_LO`].
pub const EXP_HI: f32 = 88.722_84;
/// See [`EXP_HI`].
pub const EXP_LO: f32 = -103.972_08;

/// Build `2^k` for `-126 ≤ k ≤ 127` directly from the exponent bits.
#[inline(always)]
fn pow2i(k: i32) -> f32 {
    debug_assert!((-126..=127).contains(&k));
    f32::from_bits(((k + 127) as u32) << 23)
}

/// Deterministic `exp(x)`.
///
/// Overflows to `+inf` above ≈ 88.72, underflows to `0.0` below ≈ −103.97
/// (subnormal results are produced via a two-step scale, so the tail is
/// gradual, not a cliff).
#[inline]
#[must_use]
pub fn exp(x: f32) -> f32 {
    if x.is_nan() {
        return f32::NAN;
    }
    if x > EXP_HI {
        return f32::INFINITY;
    }
    if x < EXP_LO {
        return 0.0;
    }
    let kf = nearest_int(x * LOG2E);
    let k = kf as i32;
    let r = (x - kf * LN2_HI) - kf * LN2_LO;
    let p = (((((EXP_P[0] * r + EXP_P[1]) * r + EXP_P[2]) * r + EXP_P[3]) * r + EXP_P[4]) * r
        + 0.5)
        * r
        * r
        + r
        + 1.0;
    // Scale by 2^k. Split the scale when k leaves the normal exponent range
    // so subnormal results and the overflow edge are handled without UB.
    if k > 127 {
        p * pow2i(127) * pow2i(k - 127)
    } else if k < -126 {
        p * pow2i(-126) * pow2i(k + 126)
    } else {
        p * pow2i(k)
    }
}

/// musl `logf` coefficient.
pub const LG1: f32 = 0.666_666_626_93;
/// musl `logf` coefficient.
pub const LG2: f32 = 0.400_009_721_52;
/// musl `logf` coefficient.
pub const LG3: f32 = 0.284_987_866_88;
/// musl `logf` coefficient.
pub const LG4: f32 = 0.242_790_788_41;
/// Mantissa bit pattern of `sqrt(2)` used to fold `m` into `[√2/2, √2)`.
pub const SQRT2_BITS: u32 = 0x3fb5_04f3;

/// Deterministic natural logarithm `ln(x)`.
///
/// `ln(0) = -inf`, `ln(x < 0) = NaN`, `ln(inf) = inf`.
#[inline]
#[must_use]
pub fn ln(x: f32) -> f32 {
    if x.is_nan() || x < 0.0 {
        return f32::NAN;
    }
    if x == 0.0 {
        return f32::NEG_INFINITY;
    }
    if x.is_infinite() {
        return f32::INFINITY;
    }
    // Normalise subnormals so the exponent extraction below is exact.
    let mut bits = x.to_bits();
    let mut k: i32 = 0;
    if bits < 0x0080_0000 {
        let scaled = x * 33_554_432.0; // 2^25
        bits = scaled.to_bits();
        k -= 25;
    }
    // Reduce mantissa to [sqrt(2)/2, sqrt(2)).
    k += ((bits >> 23) as i32) - 127;
    bits = (bits & 0x007f_ffff) | 0x3f80_0000; // m in [1, 2)
    if bits >= SQRT2_BITS {
        // m >= sqrt(2): halve it
        bits -= 0x0080_0000;
        k += 1;
    }
    let m = f32::from_bits(bits);
    let f = m - 1.0;
    let s = f / (2.0 + f);
    let z = s * s;
    let w = z * z;
    let t1 = w * (LG2 + w * LG4);
    let t2 = z * (LG1 + w * LG3);
    let r = t2 + t1;
    let hfsq = 0.5 * f * f;
    let kf = k as f32;
    kf * LN2_HI - ((hfsq - (s * (hfsq + r) + kf * LN2_LO)) - f)
}

/// Deterministic base-2 logarithm `log2(x)`.
///
/// [`log2_64`] rounded once to `f32` (the same pattern as
/// [`asin`] / [`acos`] / [`powf`]), so the result is within 1 ulp of correctly
/// rounded over the measured domain. A power of two returns its exponent
/// exactly. `log2(0) = -inf`, `log2(x < 0) = NaN`, `log2(inf) = inf`, and a
/// NaN input returns the canonical NaN.
#[inline]
#[must_use]
pub fn log2(x: f32) -> f32 {
    log2_64(f64::from(x)) as f32
}

/// Deterministic base-10 logarithm `log10(x)`.
///
/// [`log10_64`] rounded once to `f32`. `log10(1) = 0`
/// exactly; `log10(0) = -inf`, `log10(x < 0) = NaN`, `log10(inf) = inf`, and a
/// NaN input returns the canonical NaN.
#[inline]
#[must_use]
pub fn log10(x: f32) -> f32 {
    log10_64(f64::from(x)) as f32
}

/// Deterministic `x^y` for real `y`.
///
/// Evaluated as `exp64(y · ln64 x)` in double precision and rounded once to
/// `f32`, so the argument of the exponential carries ~29 spare bits and the
/// result is within 1 ulp of the correctly rounded value over the measured
/// domain. Special cases follow IEEE `pow` for the inputs the engine uses:
/// `x < 0 → NaN` (non-integer `y` is the common case here, so no
/// integer-exponent sign rule is attempted), `x == 0 → 0 / 1 / inf` for
/// `y > 0 / y == 0 / y < 0`, `y == 0 → 1` and `x == 1 → 1` for every `x` /
/// `y` (NaN included, as IEEE `pow`).
#[inline]
#[must_use]
pub fn powf(x: f32, y: f32) -> f32 {
    // IEEE 754: pow(x, ±0) = 1 and pow(+1, y) = 1 for every x / y, NaN
    // included (at x = 1 the exact product `y · ln x` would split an
    // infinite or near-`MAX` y and overflow to a platform-signed NaN)
    if y == 0.0 || x == 1.0 {
        return 1.0;
    }
    if x.is_nan() || y.is_nan() || x < 0.0 {
        return f32::NAN;
    }
    if x == 0.0 {
        return if y > 0.0 { 0.0 } else { f32::INFINITY };
    }
    exp64(f64::from(y) * ln64(f64::from(x))) as f32
}

/// Deterministic integer power by repeated squaring (`n` may be negative).
///
/// Unlike `f32::powi`, which lowers to a compiler-rt routine, the
/// multiplication order here is fixed by this source.
#[inline]
#[must_use]
pub fn powi(x: f32, n: i32) -> f32 {
    if n == 0 {
        return 1.0;
    }
    if x.is_nan() {
        // canonical NaN: propagating the input payload through the
        // multiplications below is platform-dependent
        return f32::NAN;
    }
    let mut base = if n < 0 { 1.0 / x } else { x };
    let mut e = n.unsigned_abs();
    let mut acc = 1.0f32;
    while e > 0 {
        if e & 1 == 1 {
            acc *= base;
        }
        e >>= 1;
        if e > 0 {
            base *= base;
        }
    }
    acc
}

/// Deterministic cube root.
///
/// Bit-hack initial estimate (`bits / 3 + 0x2a51_37a0`) followed by three
/// Newton steps, then the sign is restored. `±0`, `±inf` and `NaN` pass
/// through.
#[inline]
#[must_use]
pub fn cbrt(x: f32) -> f32 {
    if x.is_nan() {
        return f32::NAN;
    }
    if x == 0.0 || x.is_infinite() {
        return x;
    }
    let ax = x.abs();
    // Subnormals: scale up by 2^24 (an exact cube-friendly power, 2^(3·8)),
    // take the root, scale back by 2^8. Above 2^96 scale *down* by 2^24 the
    // same way: the Newton step forms `y³`, which overflows near `f32::MAX`
    // and would leave `inf / inf` = a platform-signed NaN. Power-of-two
    // scaling is exact and the bit-hack estimate scales with it, so every
    // result that was finite before is bit-identical.
    let (ax, post) = if ax.to_bits() < 0x0080_0000 {
        (ax * 16_777_216.0, 1.0 / 256.0)
    } else if ax.to_bits() >= 0x6f80_0000 {
        (ax * (1.0 / 16_777_216.0), 256.0)
    } else {
        (ax, 1.0)
    };
    let mut y = f32::from_bits(ax.to_bits() / 3 + 0x2a51_37a0);
    for _ in 0..3 {
        let y2 = y * y;
        y -= (y2 * y - ax) / (3.0 * y2);
    }
    let y = y * post;
    if x < 0.0 {
        -y
    } else {
        y
    }
}

/// Deterministic `sqrt(x² + y²)` with power-of-two scaling so intermediate
/// squares neither overflow nor flush to zero.
#[inline]
#[must_use]
pub fn hypot(x: f32, y: f32) -> f32 {
    let ax = x.abs();
    let ay = y.abs();
    if ax.is_infinite() || ay.is_infinite() {
        return f32::INFINITY;
    }
    if ax.is_nan() || ay.is_nan() {
        return f32::NAN;
    }
    let m = if ax > ay { ax } else { ay };
    if m == 0.0 {
        return 0.0;
    }
    // Scaling by powers of two is exact, so the result is bit-identical to
    // the unscaled formula whenever that formula does not over/underflow.
    let (scale, unscale) = if m > 1.0e18 {
        (
            1.0 / 18_446_744_073_709_551_616.0,
            18_446_744_073_709_551_616.0,
        ) // 2^-64, 2^64
    } else if m < 1.0e-18 {
        (
            18_446_744_073_709_551_616.0,
            1.0 / 18_446_744_073_709_551_616.0,
        )
    } else {
        (1.0, 1.0)
    };
    let sx = ax * scale;
    let sy = ay * scale;
    sqrt(sx * sx + sy * sy) * unscale
}

// ---------------------------------------------------------------------------
// f32 inverse trigonometric / tan / tanh
//
// Evaluated in double precision with fdlibm's algorithms (basic IEEE
// operations only) and rounded once to `f32`, so the `f32` result is within
// 1 ulp of correctly rounded over the measured domains (see tests).
// ---------------------------------------------------------------------------

// fdlibm s_atanf.c: atan(0.5) / atan(1) / atan(1.5) / atan(inf) high + low parts
const ATANHI_F: [f32; 4] = [
    4.636_476_040e-01,
    7.853_981_256_5e-01,
    9.827_936_887_7e-01,
    1.570_796_251_3e+00,
];
const ATANLO_F: [f32; 4] = [
    5.012_158_244e-09,
    3.774_894_707_9e-08,
    3.447_321_717e-08,
    7.549_789_415_9e-08,
];
const AT_F: [f32; 5] = [
    3.333_332_836_6e-01,
    -1.999_915_838_2e-01,
    1.425_363_570_5e-01,
    -1.064_801_737_7e-01,
    6.168_760_731_8e-02,
];

/// Deterministic `atan(x)` (fdlibm `s_atanf.c`, single precision, ≤ 1 ulp).
///
/// Until 0.2.0 this was `atan64(x as f64) as f32` (the double-precision
/// kernel rounded once); the `f32` kernel is ~3× faster and within 1 ulp.
#[inline]
#[must_use]
pub fn atan(x: f32) -> f32 {
    let ix = x.to_bits() & 0x7fff_ffff;
    if ix >= 0x4c80_0000 {
        // |x| >= 2^26
        if x.is_nan() {
            return f32::NAN;
        }
        return if x > 0.0 {
            ATANHI_F[3] + ATANLO_F[3]
        } else {
            -ATANHI_F[3] - ATANLO_F[3]
        };
    }
    let (id, t): (i32, f32) = if ix < 0x3ee0_0000 {
        // |x| < 0.4375
        if ix < 0x3980_0000 {
            // |x| < 2^-12
            return x;
        }
        (-1, x)
    } else {
        let ax = x.abs();
        if ix < 0x3f98_0000 {
            // |x| < 1.1875
            if ix < 0x3f30_0000 {
                // |x| < 0.6875
                (0, (2.0 * ax - 1.0) / (2.0 + ax))
            } else {
                (1, (ax - 1.0) / (ax + 1.0))
            }
        } else if ix < 0x401c_0000 {
            // |x| < 2.4375
            (2, (ax - 1.5) / (1.0 + 1.5 * ax))
        } else {
            (3, -1.0 / ax)
        }
    };
    let z = t * t;
    let w = z * z;
    let s1 = z * (AT_F[0] + w * (AT_F[2] + w * AT_F[4]));
    let s2 = w * (AT_F[1] + w * AT_F[3]);
    if id < 0 {
        return t - t * (s1 + s2);
    }
    let i = id as usize;
    let r = ATANHI_F[i] - ((t * (s1 + s2) - ATANLO_F[i]) - t);
    if x < 0.0 {
        -r
    } else {
        r
    }
}

const PI_F: f32 = core::f32::consts::PI;
const PI_LO_F: f32 = -8.742_277_657_3e-08;
const PIO2_HI_F: f32 = core::f32::consts::FRAC_PI_2;
const PIO4_HI_F: f32 = core::f32::consts::FRAC_PI_4;

/// Deterministic `atan2(y, x)` in `(-π, π]`, IEEE special cases as fdlibm
/// (`e_atan2f.c`, single precision, ≤ 1 ulp).
///
/// Until 0.2.0 this was the double-precision kernel rounded once.
#[inline]
#[must_use]
pub fn atan2(y: f32, x: f32) -> f32 {
    if x.is_nan() || y.is_nan() {
        return f32::NAN;
    }
    if x == 1.0 {
        return atan(y);
    }
    // m: bit0 = sign(y), bit1 = sign(x)
    let m = u32::from(y.is_sign_negative()) | (u32::from(x.is_sign_negative()) << 1);
    if y == 0.0 {
        return match m {
            0 | 1 => y, // atan(±0, +anything) = ±0
            2 => PI_F,  // atan(+0, -anything) = π
            _ => -PI_F, // atan(-0, -anything) = -π
        };
    }
    if x == 0.0 {
        return if y < 0.0 { -PIO2_HI_F } else { PIO2_HI_F };
    }
    if x.is_infinite() {
        if y.is_infinite() {
            return match m {
                0 => PIO4_HI_F,
                1 => -PIO4_HI_F,
                2 => 3.0 * PIO4_HI_F,
                _ => -3.0 * PIO4_HI_F,
            };
        }
        return match m {
            0 => 0.0,
            1 => -0.0,
            2 => PI_F,
            _ => -PI_F,
        };
    }
    if y.is_infinite() {
        return if y < 0.0 { -PIO2_HI_F } else { PIO2_HI_F };
    }
    let ex = ((x.to_bits() >> 23) & 0xff) as i32;
    let ey = ((y.to_bits() >> 23) & 0xff) as i32;
    let k = ey - ex;
    let z = if k > 26 {
        PIO2_HI_F + 0.5 * PI_LO_F
    } else if x < 0.0 && k < -26 {
        0.0
    } else {
        atan((y / x).abs())
    };
    match m {
        0 => z,
        1 => -z,
        2 => PI_F - (z - PI_LO_F),
        _ => (z - PI_LO_F) - PI_F,
    }
}

/// Deterministic `asin(x)`; `NaN` outside `[-1, 1]`.
#[inline]
#[must_use]
pub fn asin(x: f32) -> f32 {
    asin64(f64::from(x)) as f32
}

/// Deterministic `acos(x)`; `NaN` outside `[-1, 1]`.
#[inline]
#[must_use]
pub fn acos(x: f32) -> f32 {
    acos64(f64::from(x)) as f32
}

/// Deterministic `tan(x)` (`sin/cos` of the double-precision kernels).
///
/// Argument reduction is exact for `|x| < 2^20·π/2`; beyond that the result
/// is still deterministic but loses accuracy (as `sin` / `cos` do).
#[inline]
#[must_use]
pub fn tan(x: f32) -> f32 {
    if !x.is_finite() {
        return f32::NAN;
    }
    if x == 0.0 {
        return x; // ±0 → ±0, as `sin`
    }
    let (k, r) = reduce_pio2_64(f64::from(x));
    let (s, c) = (k_sin64(r), k_cos64(r));
    // tan(x + kπ/2): even k → sin/cos, odd k → -cos/sin
    let t = if k & 1 == 0 { s / c } else { -c / s };
    canon(t as f32)
}

/// Deterministic `tanh(x)`.
///
/// `|x| < 2^-14 → x` (the cubic term is below half an `f32` ulp), otherwise
/// `(e^{2|x|} − 1) / (e^{2|x|} + 1)` in double precision via [`exp64`].
#[inline]
#[must_use]
pub fn tanh(x: f32) -> f32 {
    if x.is_nan() {
        return f32::NAN;
    }
    let ax = x.abs();
    if ax < 6.103_515_625e-5 {
        return x;
    }
    if ax > 10.0 {
        // 1 - 2e^{-20} rounds to 1.0 in f32
        return if x < 0.0 { -1.0 } else { 1.0 };
    }
    let e = exp64(2.0 * f64::from(ax));
    let t = ((e - 1.0) / (e + 1.0)) as f32;
    if x < 0.0 {
        -t
    } else {
        t
    }
}

#[cfg(test)]
mod two_prod_tests {
    use super::{two_prod, PIO2_3};

    /// oracle: a 24×24-bit product is exact in `f64`, so `p + e` must equal it
    fn exact(a: f32) {
        let (p, e) = two_prod(a, PIO2_3);
        let want = f64::from(a) * f64::from(PIO2_3);
        assert_eq!(f64::from(p) + f64::from(e), want, "k = {a}");
        assert_eq!(p, a * PIO2_3, "k = {a}");
    }

    #[test]
    fn k_times_pio2_3_is_split_exactly() {
        let mut n = 0;
        for k in -(1 << 14)..=(1 << 14) {
            exact(k as f32);
            n += 1;
        }
        // a stride through 2^24, and the largest finite multiples
        for k in ((1 << 14)..=(1 << 24)).step_by(4093) {
            exact(k as f32);
            exact(-(k as f32));
            n += 2;
        }
        for k in [f32::MAX, -f32::MAX, 2.0e38, 1.0e35, 8.3e34] {
            let (p, e) = two_prod(k, PIO2_3);
            assert!(p.is_finite() && e.is_finite(), "k = {k}: {p} {e}");
            n += 1;
        }
        assert!(n > 1 << 15);
    }
}