rustymath_core/

recipe.rs

//! Formulas that rustymath's kernels compute and basedpl's JIT compiles, each written once over the operations it uses (`Ops`): tolerant closeness, comparison (each item, or through the range of floats equal to one number), floor, ceiling and the residue; Mathish division, sign, the reciprocal and `π`'s two forms; and the powers that `Exponent` routes off the general route, with the square root they take for halves. `Floats` runs a recipe on floats, as the kernels do. A compiler's `Ops` records the same operations as steps, so both routes round alike, step for step.
use crate::{Exponent, Float, float::sealed::Core, tolerant::Band};
/// `Ops`' integer methods on bit patterns held in `u64`, `bits` wide, which `NativeOps` and `IntegerOps` share.
macro_rules! u64_integers {
    ($bits:expr) => {
        #[inline(always)]
        fn integer(&self, x: i64) -> u64 {
            x as u64
        }
        #[inline(always)]
        fn integer_op(&self, op: Integer, x: u64, y: u64) -> u64 {
            crate::recipe::integer_op(op, x, y, $bits)
        }
        #[inline(always)]
        fn integer_compare(&self, op: Cmp, x: u64, y: u64) -> bool {
            crate::recipe::integer_compare(op, x, y)
        }
        #[inline(always)]
        fn integer_select(&self, test: bool, x: u64, y: u64) -> u64 {
            if test { x } else { y }
        }
    };
}
/// The comparison and Boolean `Ops` methods that `NativeOps` and `IntegerOps` share, on plain values and `bool` masks.
macro_rules! plain_logic {
    ($t:ty) => {
        #[inline(always)]
        fn compare(&self, c: Cmp, x: $t, y: $t) -> bool {
            match c {
                Cmp::Less => x < y,
                Cmp::LessEqual => x <= y,
                Cmp::Equal => x == y,
            }
        }
        #[inline(always)]
        fn and(&self, x: bool, y: bool) -> bool {
            x & y
        }
        #[inline(always)]
        fn or(&self, x: bool, y: bool) -> bool {
            x | y
        }
        #[inline(always)]
        fn not(&self, x: bool) -> bool {
            !x
        }
        #[inline(always)]
        fn select(&self, test: bool, x: $t, y: $t) -> $t {
            if test { x } else { y }
        }
    };
}
mod integer;
pub use crate::arith::Overflowing;
pub use integer::IntegerOps;

/// An elementwise operation of one value that a recipe uses, named by id: an `Ops` implementor computes each in one match (`Ops::unary`).
#[derive(Clone, Copy, Debug, PartialEq, Hash, Eq)]
pub enum Unary {
    Negate,
    Magnitude,
    Sqrt,
    Floor,
    Ln,
    Sin,
    Cos,
    Sinh,
    Cosh,
    /// √(1 - x²).
    Sqrt1mx2,
    /// √(1 + x²).
    Sqrt1px2,
    /// ψ(x).
    Digamma,
    /// e^x.
    Exp,
    /// ln |Γ(x)|.
    Lgamma,
    /// 1 - x: the negation of a truth value held as 0 or 1, as a Boolean is.
    Not,
}

/// An elementwise operation of two values that a recipe uses, named by id (`Ops::binary`).
#[derive(Clone, Copy, Debug, PartialEq, Hash, Eq)]
pub enum Binary {
    Add,
    Subtract,
    Multiply,
    Divide,
    /// Runtime power retains the same exponent dispatch as the native kernel.
    Power,
    /// `x` where `x > y`, otherwise `y`: `y` where either is NaN.
    Larger,
    /// `x` where `x < y`, otherwise `y`: `y` where either is NaN.
    Smaller,
    /// C(n, k), the binomial coefficient, with `n` first.
    Binom,
    /// ∂C(n, k)/∂n, with `n` first.
    BinomDn,
    /// ∂C(n, k)/∂k, with `n` first.
    BinomDk,
    /// ψ⁽ⁿ⁾(x), the polygamma function, with its whole order `n` first.
    Polygamma,
    /// The bitwise and of two integers, which only integer values take (`Ops::integral`).
    BitAnd,
}

/// A comparison that a recipe uses, which gives a truth value (`Ops::compare`).
#[derive(Clone, Copy, Debug, PartialEq, Hash, Eq)]
pub enum Cmp { Less, LessEqual, Equal }

#[derive(Clone, Copy, Debug, PartialEq, Hash, Eq)]
pub enum Integer { Add, Subtract, And, Or, Xor, ShiftLeft, ShiftRight, SignedShiftRight }

/// `op` on the bit patterns `x` and `y` of integers `bits` wide. Sums and differences wrap. Left and signed right shifts follow `shifted`. An unsigned right shift reads `x` as `bits` wide and gives 0 for a count of 64 or more. `NativeOps` and `IntegerOps` both compute through it.
#[inline(always)]
pub(crate) fn integer_op(op: Integer, x: u64, y: u64, bits: u32) -> u64 {
    match op {
        Integer::Add => x.wrapping_add(y), Integer::Subtract => x.wrapping_sub(y), Integer::And => x & y, Integer::Or => x | y, Integer::Xor => x ^ y,
        Integer::ShiftLeft => shifted(x, y as i64), Integer::ShiftRight => { if y >= 64 { 0 } else { (x & (u64::MAX >> (64 - bits))) >> y } }
        Integer::SignedShiftRight => shifted(x, (y as i64).saturating_neg())
    }
}

/// Comparison `op` of the bit patterns `x` and `y`, read as signed 64-bit integers.
#[inline(always)]
pub(crate) fn integer_compare(op: Cmp, x: u64, y: u64) -> bool {
    match op { Cmp::Less => (x as i64) < y as i64, Cmp::LessEqual => (x as i64) <= y as i64, Cmp::Equal => x == y }
}

/// `x`'s bits shifted left by `count`, or right with the sign filled in for a negative count. A count of 64 or more saturates: to 0 leftward, and to the sign rightward. Items hold their sign extended to 64 bits, so a count past an item's own width saturates the same way once the result narrows.
#[inline(always)]
fn shifted(x: u64, count: i64) -> u64 { if count >= 64 { 0 } else if count >= 0 { x << count } else { ((x as i64) >> count.unsigned_abs().min(63)) as u64 } }

/// The named methods of `Ops`, each the operation of its id. Inline metadata keeps each constant operation id visible across crate boundaries, so native item loops can specialize and vectorize.
macro_rules! named {
    (unary: $($u:ident $uid:ident),+; binary: $($b:ident $bid:ident),+; compare: $($c:ident $cid:ident),+) => {
        $( #[inline] fn $u(&self, x: Self::V) -> Self::V { self.unary(Unary::$uid, x) } )+
        $( #[inline] fn $b(&self, x: Self::V, y: Self::V) -> Self::V { self.binary(Binary::$bid, x, y) } )+
        $( #[inline] fn $c(&self, x: Self::V, y: Self::V) -> Self::B { self.compare(Cmp::$cid, x, y) } )+
    };
}

/// The operations that a recipe uses, on values `V` and truth values `B`, at float type `F`. Each rounds as IEEE's operation, or rustymath's function of its name, does at `F`'s width. An implementor computes the operations by id, in `unary`, `binary` and `compare`, and recipes call them by their names.
pub trait Ops<F: Float> {
    /// A value.
    type V: Copy;
    /// A truth value.
    type B: Copy;
    /// Integer bits at the resolved working width.
    type I: Copy;
    fn width(&self) -> crate::buffer::FloatWidth;
    fn fused(&self) -> bool;
    /// A binary64 temporary for short mixed-width reductions.
    type W: Copy;
    fn widen(&self, x: Self::V) -> Self::W;
    fn wide_number(&self, x: f64) -> Self::W;
    fn wide_multiply_add(&self, x: Self::W, y: Self::W, z: Self::W) -> Self::W;
    fn narrow(&self, x: Self::W) -> Self::V;
    fn integer(&self, x: i64) -> Self::I;
    fn integer_op(&self, op: Integer, x: Self::I, y: Self::I) -> Self::I;
    fn integer_compare(&self, op: Cmp, x: Self::I, y: Self::I) -> Self::B;
    fn integer_select(&self, test: Self::B, x: Self::I, y: Self::I) -> Self::I;
    fn to_integer(&self, x: Self::V) -> Self::I;
    fn from_integer(&self, x: Self::I) -> Self::V;
    fn to_bits(&self, x: Self::V) -> Self::I;
    fn from_bits(&self, x: Self::I) -> Self::V;
    /// The number `x`.
    fn number(&self, x: F) -> Self::V;
    fn unary(&self, op: Unary, x: Self::V) -> Self::V;
    fn binary(&self, op: Binary, x: Self::V, y: Self::V) -> Self::V;
    fn compare(&self, c: Cmp, x: Self::V, y: Self::V) -> Self::B;
    /// `x·y + z`, rounded once where the target fuses multiply-add (`FUSED`), as rustymath's `mla` gives it.
    fn multiply_add(&self, x: Self::V, y: Self::V, z: Self::V) -> Self::V;
    fn overflow(&self, op: Binary, x: Self::V, y: Self::V) -> Self::B;
    fn and(&self, x: Self::B, y: Self::B) -> Self::B;
    fn or(&self, x: Self::B, y: Self::B) -> Self::B;
    fn not(&self, x: Self::B) -> Self::B;
    /// `x` where `m` holds, otherwise `y`. `∂`'s recorder runs a recipe on whole values, where `m` and each array case must have one layout: build them from the same arguments.
    fn select(&self, m: Self::B, x: Self::V, y: Self::V) -> Self::V;
    /// 1 where `m` holds, otherwise 0.
    #[inline]
    fn one_where(&self, m: Self::B) -> Self::V { self.select(m, self.number(F::one()), self.number(F::zero())) }
    /// `v` where it's one number that every item shares, known before any item is read, such as the 3 of `x^3`. A recipe can pick a cheaper formula by it. `None` where `v` varies by item.
    #[inline]
    fn constant(&self, _v: Self::V) -> Option<F> { None }
    /// Whether values are integers, as a compiler's values at an integer type are. A recipe for an integer object asks before taking an integer-only step (`Binary::BitAnd`). Floats never are.
    #[inline]
    fn integral(&self) -> bool { false }
    /// Captured tolerance; recording backends return its runtime setting slot.
    #[inline]
    fn tolerance(&self, fallback: F) -> Self::V { self.number(fallback) }
    /// The sign of `x`, as [`sign`] gives it.
    #[inline]
    fn sign(&self, x: Self::V) -> Self::V { sign(self, x) }
    /// `1 ÷ x`, as [`reciprocal`] gives it.
    #[inline]
    fn reciprocal(&self, x: Self::V) -> Self::V { reciprocal(self, x) }
    /// `x` raised to `k`, a constant that `Exponent::of` routes off the general route, as [`power`] gives it.
    #[inline]
    fn power(&self, k: F, x: Self::V) -> Self::V { power(self, Exponent::of(k.wide(), self.fused()), x).expect("an exponent off the general route") }
    /// The real part of `x^y` through the existing power formula.
    #[inline]
    fn real_power(&self, x: Self::V, y: Self::V) -> Self::V { crate::pow::real_formula(self, x, y) }
    /// The partial in `x` of `x^y = r`; `halves` permits constant half-integer powers.
    #[inline]
    fn power_slope(&self, x: Self::V, y: Self::V, r: Self::V, halves: bool) -> Self::V { crate::derivative::power_slope(self, x, y, r, halves) }
    named! {
        unary: negate Negate, magnitude Magnitude, sqrt Sqrt, floor Floor, ln Ln, sin Sin, cos Cos, sinh Sinh, cosh Cosh, sqrt1mx2 Sqrt1mx2, sqrt1px2 Sqrt1px2, digamma Digamma, exp Exp, lgamma Lgamma, complement Not;
        binary: add Add, subtract Subtract, multiply Multiply, divide Divide, larger Larger, smaller Smaller, binom Binom, binom_dn BinomDn, binom_dk BinomDk, polygamma Polygamma, bit_and BitAnd;
        compare: less Less, less_equal LessEqual, equal Equal
    }
}

/// Conservative range interpretation of the shared predicates.
pub struct RangeOps(pub crate::buffer::DType);
impl RangeOps {
    fn compare(op: Cmp, x: crate::domain::Range, y: crate::domain::Range) -> Option<bool> {
        if x.is_empty() || y.is_empty() { return None; }
        match op {
            Cmp::Less if x.hi < y.lo => Some(true), Cmp::Less if x.lo >= y.hi => Some(false), Cmp::LessEqual if x.hi <= y.lo => Some(true),
            Cmp::LessEqual if x.lo > y.hi => Some(false), Cmp::Equal if x.constant().is_some() && x == y => Some(true),
            Cmp::Equal if x.hi < y.lo || y.hi < x.lo => Some(false), _ => None
        }
    }
    fn selected(test: Option<bool>, x: crate::domain::Range, y: crate::domain::Range) -> crate::domain::Range {
        match test { Some(true) => x, Some(false) => y, None => x.join(y) }
    }
}
impl<F: Float> Ops<F> for RangeOps {
    type V = crate::domain::Range;
    type B = Option<bool>;
    type I = crate::domain::Range;
    type W = crate::domain::Range;
    fn width(&self) -> crate::buffer::FloatWidth { self.0.numpy_float() }
    fn fused(&self) -> bool { crate::FUSED }
    fn integral(&self) -> bool { self.0.is_integral() }
    fn widen(&self, _: Self::V) -> Self::W { crate::domain::Range::ANY }
    fn wide_number(&self, x: f64) -> Self::W { crate::domain::Range::of(Some(x)) }
    fn wide_multiply_add(&self, _: Self::W, _: Self::W, _: Self::W) -> Self::W { crate::domain::Range::ANY }
    fn narrow(&self, _: Self::W) -> Self::V { crate::domain::Range::ANY }
    fn integer(&self, x: i64) -> Self::I { crate::domain::Range::of((x.unsigned_abs() <= 1 << 53).then_some(x as f64)) }
    fn integer_op(&self, _: Integer, _: Self::I, _: Self::I) -> Self::I { crate::domain::Range::ANY }
    fn integer_compare(&self, op: Cmp, x: Self::I, y: Self::I) -> Self::B { Self::compare(op, x, y) }
    fn integer_select(&self, test: Self::B, x: Self::I, y: Self::I) -> Self::I { Self::selected(test, x, y) }
    fn to_integer(&self, _: Self::V) -> Self::I { crate::domain::Range::ANY }
    fn from_integer(&self, _: Self::I) -> Self::V { crate::domain::Range::ANY }
    fn to_bits(&self, _: Self::V) -> Self::I { crate::domain::Range::ANY }
    fn from_bits(&self, _: Self::I) -> Self::V { crate::domain::Range::ANY }
    fn number(&self, x: F) -> Self::V { crate::domain::Range::of(Some(x.wide())) }
    fn constant(&self, x: Self::V) -> Option<F> { x.constant().map(F::narrow) }
    fn unary(&self, op: Unary, x: Self::V) -> Self::V { if op == Unary::Magnitude { x.magnitude() } else { crate::domain::Range::ANY } }
    fn binary(&self, _: Binary, _: Self::V, _: Self::V) -> Self::V { crate::domain::Range::ANY }
    fn compare(&self, op: Cmp, x: Self::V, y: Self::V) -> Self::B { Self::compare(op, x, y) }
    fn multiply_add(&self, _: Self::V, _: Self::V, _: Self::V) -> Self::V { crate::domain::Range::ANY }
    /// No overflow (`Some(false)`) where the operation's range, as its formula's facts give it, lies within the computing dtype's values. Past 2⁵³ a bound and the result's ends round, so a result stays strictly inside such a bound.
    fn overflow(&self, op: Binary, x: Self::V, y: Self::V) -> Self::B {
        use crate::facts::Formula;
        let result = match op {
            Binary::Add => crate::Add.range([x, y]), Binary::Subtract => crate::Subtract.range([x, y]), Binary::Multiply => crate::Multiply.range([x, y]),
            _ => return None
        };
        let fits = crate::domain::Range::of_dtype(self.0);
        let exact = |bound: f64| bound.abs() < 2f64.powi(53);
        let above = if exact(fits.lo) { result.lo >= fits.lo } else { result.lo > fits.lo };
        let below = if exact(fits.hi) { result.hi <= fits.hi } else { result.hi < fits.hi };
        (above && below).then_some(false)
    }
    fn and(&self, x: Self::B, y: Self::B) -> Self::B {
        match (x, y) { (Some(false), _) | (_, Some(false)) => Some(false), (Some(true), Some(true)) => Some(true), _ => None }
    }
    fn or(&self, x: Self::B, y: Self::B) -> Self::B {
        match (x, y) { (Some(true), _) | (_, Some(true)) => Some(true), (Some(false), Some(false)) => Some(false), _ => None }
    }
    fn not(&self, x: Self::B) -> Self::B { x.map(|x| !x) }
    fn select(&self, test: Self::B, x: Self::V, y: Self::V) -> Self::V { Self::selected(test, x, y) }
}

/// Scalar instantiation of the same formulae used by slice kernels and graph lowering.
#[derive(Clone, Copy, Debug, Default)]
pub struct NativeOps<const FUSED: bool> { pub tolerance: Option<f64> }
impl<const FUSED: bool> NativeOps<FUSED> { pub const DEFAULT: Self = Self { tolerance: None }; }
pub type Floats = NativeOps<{ crate::FUSED }>;
#[allow(non_upper_case_globals)]
pub const Floats: Floats = NativeOps::DEFAULT;

impl<F: Float, const FUSED: bool> Ops<F> for NativeOps<FUSED> {
    type V = F;
    type B = bool;
    type I = u64;
    type W = f64;
    #[inline(always)]
    fn fused(&self) -> bool { FUSED }
    #[inline(always)]
    fn widen(&self, x: F) -> f64 { x.wide() }
    #[inline(always)]
    fn wide_number(&self, x: f64) -> f64 { x }
    #[inline(always)]
    fn wide_multiply_add(&self, x: f64, y: f64, z: f64) -> f64 { if FUSED { x.mul_add(y, z) } else { x * y + z } }
    #[inline(always)]
    fn narrow(&self, x: f64) -> F { F::narrow(x) }
    #[inline(always)]
    fn width(&self) -> crate::buffer::FloatWidth { if F::Working::WIDE { crate::buffer::FloatWidth::F64 } else { crate::buffer::FloatWidth::F32 } }
    u64_integers!(64);
    #[inline(always)]
    fn to_integer(&self, x: F) -> u64 { (x.wide() as i64) as u64 }
    #[inline(always)]
    fn from_integer(&self, x: u64) -> F { F::narrow((x as i64) as f64) }
    #[inline(always)]
    fn to_bits(&self, x: F) -> u64 { if F::Working::WIDE { x.wide().to_bits() } else { (x.wide() as f32).to_bits() as u64 } }
    #[inline(always)]
    fn from_bits(&self, x: u64) -> F { F::narrow(if F::Working::WIDE { f64::from_bits(x) } else { f32::from_bits(x as u32) as f64 }) }
    #[inline(always)]
    fn number(&self, x: F) -> F { x }
    #[inline(always)]
    fn unary(&self, op: Unary, x: F) -> F {
        match op {
            Unary::Negate => -x, Unary::Magnitude => x.abs(), Unary::Sqrt => x.sqrt(), Unary::Floor => x.floor(), Unary::Ln => crate::ln(x),
            Unary::Sin => crate::sin(x), Unary::Cos => crate::cos(x), Unary::Sinh => crate::sinh(x), Unary::Cosh => crate::cosh(x),
            Unary::Sqrt1mx2 => crate::sqrt1mx2(x), Unary::Sqrt1px2 => crate::sqrt1px2(x), Unary::Digamma => crate::digamma(x), Unary::Exp => crate::exp(x),
            Unary::Lgamma => crate::lgamma(x), Unary::Not => F::one() - x
        }
    }
    #[inline(always)]
    fn binary(&self, op: Binary, x: F, y: F) -> F {
        match op {
            Binary::Add => x + y, Binary::Subtract => x - y, Binary::Multiply => x * y, Binary::Divide => x / y, Binary::Larger => { if x > y { x } else { y } }
            Binary::Smaller => { if x < y { x } else { y } } Binary::Binom => crate::binom(x, y), Binary::BinomDn => crate::binom_dn(x, y),
            Binary::Power => F::from_working(crate::pow::pow::<F::Working, FUSED>(x.working(), y.working())), Binary::BinomDk => crate::binom_dk(x, y),
            Binary::Polygamma => crate::polygamma(x, y),
            Binary::BitAnd => unreachable!("floats take no bitwise and: a recipe asks for one only where `Ops::integral` holds")
        }
    }
    plain_logic!(F);
    #[inline(always)]
    fn multiply_add(&self, x: F, y: F, z: F) -> F { F::from_working(x.working().mla::<FUSED>(y.working(), z.working())) }
    #[inline(always)]
    fn overflow(&self, _op: Binary, _x: F, _y: F) -> bool { false }
    #[inline(always)]
    fn constant(&self, value: F) -> Option<F> { Some(value) }
    #[inline(always)]
    fn tolerance(&self, fallback: F) -> F { self.tolerance.map_or(fallback, F::narrow) }
}

/// Whether `x` and `y` lie within tolerance `t` of each other: `|x-y| < t·max(|x|, |y|)`. With a strict `<`, an infinity is close to no number without a finiteness test, because `inf < inf` is false. NaN is close to nothing.
#[inline(always)]
pub fn close<F: Float, O: Ops<F> + ?Sized>(o: &O, t: O::V, x: O::V, y: O::V) -> O::B {
    let gap = o.subtract(x, y);
    let gap = o.magnitude(gap);
    let (a, b) = (o.magnitude(x), o.magnitude(y));
    let size = o.larger(a, b);
    let bound = o.multiply(t, size);
    o.less(gap, bound)
}

/// Whether `x` equals `y`: `x == y`, or `close`, which gives whether they're close. Each of these compositions runs `close` after its plain comparison: in the kernels' loops, that order is what lets LLVM vectorize 16 items at a time rather than 8.
#[inline(always)]
pub fn equal<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V, y: O::V, close: impl FnOnce(&O) -> O::B) -> O::B {
    if o.integral() { return o.equal(x, y); }
    let same = o.equal(x, y);
    let close = close(o);
    o.or(same, close)
}

/// Whether `x` is below `y` and doesn't equal it, where `close` gives whether they're close. `x < y` already implies `x != y`, so only `close` is tested.
#[inline(always)]
pub fn less<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V, y: O::V, close: impl FnOnce(&O) -> O::B) -> O::B {
    if o.integral() { return o.less(x, y); }
    let below = o.less(x, y);
    let close = close(o);
    let apart = o.not(close);
    o.and(below, apart)
}

/// Whether `x` is below `y` or equals it, where `close` gives whether they're close. `x <= y` already holds where `x == y`, so only `close` is added.
#[inline(always)]
pub fn less_equal<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V, y: O::V, close: impl FnOnce(&O) -> O::B) -> O::B {
    if o.integral() { return o.less_equal(x, y); }
    let below = o.less_equal(x, y);
    let close = close(o);
    o.or(below, close)
}

/// The floor of `y` within tolerance `t`. `y` equals the whole number above it when raising `y` by `t·|y|` reaches that number. Capping the raise just below a half keeps a whole number from reaching the next one where `t·|y|` exceeds a half. Where floats are 1 or more apart, the sum rounds back to `y`. A zero result is `0`, since `-0` plus a zero raise is `0`.
#[inline(always)]
pub fn floor<F: Float, O: Ops<F> + ?Sized>(o: &O, t: O::V, y: O::V) -> O::V {
    if o.integral() { return y; }
    let size = o.magnitude(y);
    let raise = o.multiply(t, size);
    let most = o.number(F::narrow(0.5).next_down());
    let raise = o.smaller(raise, most);
    let raised = o.add(y, raise);
    o.floor(raised)
}

#[inline(always)]
pub(crate) fn integer_floor_divide<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V, y: O::V) -> O::V {
    let q = o.divide(x, y);
    let r = o.subtract(x, o.multiply(q, y));
    let zero = o.number(F::zero());
    let lower = o.or(o.and(o.less(r, zero), o.less(zero, y)), o.and(o.less(zero, r), o.less(y, zero)));
    o.subtract(q, o.one_where(lower))
}

/// The exact floor of `y`. Adding 0 turns IEEE's `-0` into `0`.
#[inline(always)]
pub fn exact_floor<F: Float, O: Ops<F> + ?Sized>(o: &O, y: O::V) -> O::V {
    let whole = o.floor(y);
    let zero = o.number(F::zero());
    o.add(whole, zero)
}

/// The ceiling of `y` through `floor`, the floor at the same tolerance: `0 - floor(-y)`.
#[inline(always)]
pub fn ceil<F: Float, O: Ops<F> + ?Sized>(o: &O, y: O::V, floor: impl FnOnce(&O, O::V) -> O::V) -> O::V {
    if o.integral() { return y; }
    let negated = o.negate(y);
    let whole = floor(o, negated);
    let zero = o.number(F::zero());
    o.subtract(zero, whole)
}

/// `y` modulo `x`: `y - x·q`, for `q` the floor of `y÷x` that `floor` gives, rounded once where the target fuses multiply-add. Within a tolerance `t`, a product `x·q` that equals `y` leaves no remainder. `x` = 0 gives `y`. At no tolerance, where `y÷x` passes the float's precision, the rounded quotient has lost its low digits: `tolerant::Residue` recomputes those items.
#[inline(always)]
pub fn residue<F: Float, O: Ops<F> + ?Sized>(o: &O, t: Option<O::V>, x: O::V, y: O::V, floor: impl FnOnce(&O, O::V) -> O::V) -> O::V {
    let quotient = o.divide(y, x);
    let q = floor(o, quotient);
    let minus = o.negate(x);
    let r = o.multiply_add(minus, q, y);
    let r = match t {
        Some(t) => {
            let product = o.multiply(x, q);
            let whole = equal(o, product, y, |o| close(o, t, product, y));
            let zero = o.number(F::zero());
            o.select(whole, zero, r)
        }
        None => r,
    };
    let zero = o.number(F::zero());
    let none = o.equal(x, zero);
    o.select(none, y, r)
}

/// The square root of `y`, with the values that `pow` gives `y^0.5` at -0 and -∞: 0 and ∞.
#[inline(always)]
pub fn root2<F: Float, O: Ops<F> + ?Sized>(o: &O, y: O::V) -> O::V {
    let root = o.sqrt(y);
    let zero = o.number(F::zero());
    let root = o.add(root, zero);
    let low = o.number(F::neg_infinity());
    let infinite = o.equal(y, low);
    let flipped = o.negate(y);
    o.select(infinite, flipped, root)
}

/// The real part of √y: √y where y is positive, 0 at or below 0, and NaN at NaN, which `larger` keeps.
#[inline(always)]
pub fn real_sqrt<F: Float, O: Ops<F> + ?Sized>(o: &O, y: O::V) -> O::V {
    let zero = o.number(F::zero());
    let clamped = o.larger(zero, y);
    o.sqrt(clamped)
}

/// (sin πx, cos πx), exact at whole and half x: x = n/2 + r with n the whole number nearest 2x and |r| up to ¼, where n mod 4 picks the quadrant.
#[inline(always)]
pub fn sincospi<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V) -> (O::V, O::V) {
    let (half, quarter) = (o.number(F::narrow(0.5)), o.number(F::narrow(0.25)));
    let n = o.floor(o.add(o.add(x, x), half));
    let a = o.multiply(o.subtract(x, o.multiply(n, half)), o.number(F::narrow(core::f64::consts::PI)));
    let (s, k) = (o.sin(a), o.cos(a));
    let q = o.subtract(n, o.multiply(o.number(F::narrow(4.0)), o.floor(o.multiply(n, quarter))));
    let [q0, q1, q2] = [0.0, 1.0, 2.0].map(|m| o.equal(q, o.number(F::narrow(m))));
    // 0 - v negates v exactly, and gives +0, never ¯0, at an exact half step.
    let zero = o.number(F::zero());
    let (ns, nk) = (o.subtract(zero, s), o.subtract(zero, k));
    (o.select(q0, s, o.select(q1, k, o.select(q2, ns, nk))), o.select(q0, k, o.select(q1, ns, o.select(q2, nk, s))))
}

/// Whether `v` lies in `band` against the range `lo` to `hi` of floats equal to one number (`Tolerant::range`), by one or two plain comparisons.
#[inline(always)]
pub fn in_band<F: Float, O: Ops<F> + ?Sized>(o: &O, band: Band, lo: O::V, hi: O::V, v: O::V) -> O::B {
    match band {
        Band::Inside | Band::Outside => {
            let (above, below) = (o.less_equal(lo, v), o.less_equal(v, hi));
            let inside = o.and(above, below);
            if matches!(band, Band::Outside) { o.not(inside) } else { inside }
        }
        Band::Below => o.less(v, lo),
        Band::NotAbove => o.less_equal(v, hi),
        Band::Above => o.less(hi, v),
        Band::NotBelow => o.less_equal(lo, v),
    }
}

/// `1 ÷ x`.
#[inline(always)]
pub fn reciprocal<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V) -> O::V { let one = o.number(F::one()); o.divide(one, x) }

/// `π × x`, with π at `F`'s width.
#[inline(always)]
pub fn pi_times<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V) -> O::V { let pi = o.number(F::narrow(std::f64::consts::PI)); o.multiply(x, pi) }

/// `π × (x ÷ y)`, an angle of `x` over `y` turns of half a circle.
#[inline(always)]
pub fn pi_ratio<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V, y: O::V) -> O::V { let ratio = o.divide(x, y); pi_times(o, ratio) }

/// `x ÷ y`, or `x` where both are zero, as Mathish division gives `0÷0`.
#[inline(always)]
pub fn mathish_divide<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V, y: O::V) -> O::V {
    let quotient = o.divide(x, y);
    let zero = o.number(F::zero());
    let (a, b) = (o.equal(x, zero), o.equal(y, zero));
    let both = o.and(a, b);
    o.select(both, x, quotient)
}

/// The sign of `x`: `-1`, `1`, or `x` itself where `x` is a zero or NaN.
#[inline(always)]
pub fn sign<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V) -> O::V {
    let (zero, one, minus) = (o.number(F::zero()), o.number(F::one()), o.number(-F::one()));
    let (above, below) = (o.less(zero, x), o.less(x, zero));
    o.select(above, one, o.select(below, minus, x))
}

/// `x` raised to an exponent admitted by `e`. Large and negative powers use [`whole`]; small positive powers use [`multiplied`]. Halves first take `root2`. `None` for `Exponent::General`.
#[inline(always)]
pub fn power<F: Float, O: Ops<F> + ?Sized>(o: &O, e: Exponent, x: O::V) -> Option<O::V> {
    let (m, negative, x) = match e {
        Exponent::Whole { m, negative } => (m, negative, x), Exponent::Halves { m, negative } => (m, negative, root2(o, x)), Exponent::General => return None
    };
    Some(if !o.integral() && (m > crate::pow::MULTIPLIED || (negative && m > 1)) { whole(o, m, negative, x) } else { multiplied(o, m, negative, x) })
}

/// Positive powers up to `MULTIPLIED`, squaring from m's lowest set bit, and the reciprocal for m = 1. A rounding at the step that gives x^k carries into x^m about m/k times. A result that the float holds exactly comes out exactly.
#[inline(always)]
fn multiplied<F: Float, O: Ops<F> + ?Sized>(o: &O, m: u32, negative: bool, x: O::V) -> O::V {
    let (mut base, mut k) = (x, m);
    while k & 1 == 0 {
        base = o.multiply(base, base);
        k >>= 1;
    }
    let mut r = base;
    while k > 1 {
        k >>= 1;
        base = o.multiply(base, base);
        if k & 1 == 1 { r = o.multiply(r, base); }
    }
    if negative { reciprocal(o, r) } else { r }
}

/// A power that [`whole`] has part computed: `h + l` in double words, or `h` alone where `l` is `None`.
pub(crate) type Word<V> = (V, Option<V>);

/// Powers admitted for fused multiply-add. Negative powers start with the reciprocal and its rounding error, avoiding an overflowing positive power whose reciprocal is still representable. Raises [`base`] by squaring and multiplying from m's highest set bit, with a [`step`] for each lower bit. A result that the float holds exactly comes out exactly.
#[inline(always)]
pub(crate) fn whole<F: Float, O: Ops<F> + ?Sized>(o: &O, m: u32, negative: bool, x: O::V) -> O::V {
    let (b, mut error) = base(o, negative, x);
    let mut w = (b, None);
    for i in (0..u32::BITS - 1 - m.leading_zeros()).rev() { w = step(o, w, b, m, i, &mut error); }
    w.0
}

/// The number that [`whole`] raises to m: x, or q = 1/x with its relative error ε = 1 - x·q where `negative`, since 1/x is q·(1 + ε) to within ε².
#[inline(always)]
pub(crate) fn base<F: Float, O: Ops<F> + ?Sized>(o: &O, negative: bool, x: O::V) -> (O::V, Option<O::V>) {
    if !negative { return (x, None); }
    let q = reciprocal(o, x);
    let one = o.number(F::one());
    (q, Some(o.multiply_add(o.negate(x), q, one)))
}

/// The steps of bit position `i` of [`whole`]'s m: `w` squared, then times the base `b` where that bit is set.
#[inline(always)]
pub(crate) fn step<F: Float, O: Ops<F> + ?Sized>(o: &O, w: Word<O::V>, b: O::V, m: u32, i: u32, error: &mut Option<O::V>) -> Word<O::V> {
    let w = times(o, w, None, m, i, error);
    if (m >> i) & 1 == 1 { times(o, w, Some(b), m, i, error) } else { w }
}

/// One step of [`whole`] at bit position `i`: `w` squared, or times `b` for `by` = `Some(b)`. A rounding at the step that gives x^k carries into x^m about m/k times. Each step at a position above 0 runs in double words, taking its product's rounding error exactly from a fused multiply-add. The first step after them rounds the double word to one word. That step also adds the base's error `error`, times m and halved once for each square from that step on. This turns a power of q into a power of 1/x. The steps after it run in single words. The roundings left keep x^m within 4 ulps.
#[inline(always)]
fn times<F: Float, O: Ops<F> + ?Sized>(o: &O, (h, l): Word<O::V>, by: Option<O::V>, m: u32, i: u32, error: &mut Option<O::V>) -> Word<O::V> {
    let y = by.unwrap_or(h);
    // The coefficient of l in the product: 2h for a square, b for a product with b.
    let twice = |o: &O| by.unwrap_or_else(|| o.add(h, h));
    if i > 0 {
        let p = o.multiply(h, y);
        let e = o.multiply_add(h, y, o.negate(p));
        return (p, Some(match l { Some(l) => o.multiply_add(twice(o), l, e), None => e }));
    }
    let squares = i + by.is_none() as u32;
    let fold = error.take().map(|e| o.multiply(e, o.number(F::narrow(m as f64 / (1u64 << squares) as f64))));
    let l = match (l, fold) {
        (None, None) => return (o.multiply(h, y), None), (Some(l), None) => l, (None, Some(c)) => o.multiply(h, c),
        (Some(l), Some(c)) => o.multiply_add(h, c, l)
    };
    let k = o.multiply(twice(o), l);
    // k is infinite or NaN only where h is, or where h·y overflows. h·y alone is the power there.
    let k = o.select(finite(o, k), k, o.number(F::zero()));
    (o.multiply_add(h, y, k), None)
}

/// Whether `v` is finite: `v - v` is 0 there, and NaN at an infinity or NaN.
#[inline(always)]
fn finite<F: Float, O: Ops<F> + ?Sized>(o: &O, v: O::V) -> O::B { o.equal(o.subtract(v, v), o.number(F::zero())) }

#[inline(always)]
pub fn horner<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V, coefficients: &[f64]) -> O::V {
    let mut result = o.number(F::narrow(coefficients[0]));
    for &coefficient in &coefficients[1..] { result = o.multiply_add(result, x, o.number(F::narrow(coefficient))); }
    result
}

#[inline(always)]
pub fn ln2<F: Float, O: Ops<F> + ?Sized>(o: &O) -> (O::V, O::V) {
    let (hi, lo) = if o.width() == crate::buffer::FloatWidth::F64 {
        (0.69314718055966295651160180568695068359375, 0.28235290563031577122588448175013436025525412068e-12)
    } else { (0.693145751953125, 1.428606765330187e-06) };
    (o.number(F::narrow(hi)), o.number(F::narrow(lo)))
}

#[inline(always)]
pub fn rounding_shift<F: Float, O: Ops<F> + ?Sized>(o: &O) -> O::V {
    o.number(F::narrow(if o.width() == crate::buffer::FloatWidth::F64 { 6755399441055744.0 } else { 12582912.0 }))
}

/// 2^n from z = rounding_shift + n, where n lies within the normal exponent range.
#[inline(always)]
pub fn pow2<F: Float, O: Ops<F> + ?Sized>(o: &O, z: O::V) -> O::V {
    let (mantissa, bias) = if o.width() == crate::buffer::FloatWidth::F64 { (52, 1023) } else { (23, 127) };
    let bits = o.integer_op(Integer::Add, o.to_bits(z), o.integer(bias));
    o.from_bits(o.integer_op(Integer::ShiftLeft, bits, o.integer(mantissa)))
}

/// y times 2^n, with one rounding when the result is subnormal.
#[inline(always)]
pub fn ldexp<F: Float, O: Ops<F> + ?Sized>(o: &O, y: O::V, n: O::V) -> O::V {
    let shift = rounding_shift(o);
    let half = o.multiply_add(n, o.number(F::narrow(0.5)), shift);
    let rest = o.add(o.subtract(n, o.subtract(half, shift)), shift);
    o.multiply(o.multiply(y, pow2(o, half)), pow2(o, rest))
}

/// Positive normal x = m times 2^e, with m in [0.75, 1.5].
#[inline(always)]
pub fn split<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V) -> (O::V, O::V) {
    let (mantissa, bias) = if o.width() == crate::buffer::FloatWidth::F64 { (52, 1023) } else { (23, 127) };
    let scaled = o.multiply(x, o.number(F::narrow(1.0 / 0.75)));
    let exponent = o.integer_op(Integer::ShiftRight, o.to_bits(scaled), o.integer(mantissa));
    let exponent = o.integer_op(Integer::And, exponent, o.integer(2 * bias + 1));
    let unbiased = o.integer_op(Integer::Subtract, exponent, o.integer(bias));
    let offset = o.integer_op(Integer::ShiftLeft, unbiased, o.integer(mantissa));
    (o.from_bits(o.integer_op(Integer::Subtract, o.to_bits(x), offset)), o.from_integer(unbiased))
}

/// Positive x = m times 2^e, including subnormals.
#[inline(always)]
pub fn frexp<F: Float, O: Ops<F> + ?Sized>(o: &O, x: O::V) -> (O::V, O::V) {
    let smallest = if o.width() == crate::buffer::FloatWidth::F64 { f64::MIN_POSITIVE } else { f32::MIN_POSITIVE as f64 };
    let small = o.less(x, o.number(F::narrow(smallest)));
    let scale = o.select(small, o.number(F::narrow(18446744073709551616.0)), o.number(F::one()));
    let (m, e) = split(o, o.multiply(x, scale));
    (m, o.subtract(e, o.select(small, o.number(F::narrow(64.0)), o.number(F::zero()))))
}