rigless 0.3.0

// Copyright (c) 2026 Rigless Labs.
// SPDX-License-Identifier: MIT OR Apache-2.0
//! The floating-point extension's arithmetic, DDI 0403E.b A2.5, on `rustc_apfloat`'s software
//! IEEE 754, so a result is the same bits on every host.
//!
//! Each operation takes and answers registers as bits, a word for a single and a doubleword for a
//! double, reads its controls from the `FPSCR` you pass and sets its cumulative exception bits
//! there, A2.5.3. `rustc_apfloat` does the rounding. What the architecture adds to IEEE 754 is
//! done here: flush-to-zero, the default NaN, which NaN an operation hands on, and underflow
//! judged before rounding, A2.5.4.

use std::cmp::Ordering;
use std::num::TryFromIntError;

use rustc_apfloat::ieee::{Double, Half, Quad, Single};
use rustc_apfloat::{ExpInt, Float, FloatConvert, Round, Status, StatusAnd};

use super::registers::{fpscr, get, is_set, mask};
use crate::soc::FieldSpec;

/// `FPSCR.N`, `Z`, `C` and `V`, A2.5.3.
pub(crate) const FPSCR_NZCV: u32 =
    mask(fpscr::N) | mask(fpscr::Z) | mask(fpscr::C) | mask(fpscr::V);
/// `FPSCR.AHP`, `DN`, `FZ` and `RMode`, which a new floating-point context takes from `FPDSCR`.
pub(crate) const FPSCR_DEFAULTS: u32 =
    mask(fpscr::AHP) | mask(fpscr::DN) | mask(fpscr::FZ) | mask(fpscr::RMODE);
/// Every bit of `FPSCR` that isn't reserved.
pub(crate) const FPSCR_WRITABLE: u32 = FPSCR_NZCV
    | FPSCR_DEFAULTS
    | mask(fpscr::IDC)
    | mask(fpscr::IXC)
    | mask(fpscr::UFC)
    | mask(fpscr::OFC)
    | mask(fpscr::DZC)
    | mask(fpscr::IOC);

/// A single's sign bit, and a double's in its upper word.
pub(crate) const SIGN: u32 = 0x8000_0000;

const HALF_SIGN: u32 = 0x8000;
const HALF_EXPONENT: u32 = 0x7C00;
const HALF_FRACTION: u32 = 0x03FF;
const HALF_QUIET: u32 = 0x0200;
/// `FPDefaultNaN(16)`.
const HALF_DEFAULT_NAN: u32 = 0x7E00;
/// What the alternative half-precision format gives for an infinity or a value too big for it.
const HALF_ALTERNATIVE_MAX: u32 = 0x7FFF;
/// The bits of a half's fraction, which is where its exponent field starts.
const HALF_FRACTION_BITS: u32 = 10;
/// The half-precision format's exponent bias, A2.5.4.
const HALF_BIAS: ExpInt = 15;
/// The least exponent the alternative half-precision format needs its top exponent field for.
const HALF_ALTERNATIVE_TOP: ExpInt = 15;
/// The least exponent the alternative format can't hold in any mode.
const HALF_ALTERNATIVE_BEYOND: ExpInt = 17;

/// A register's bits, whose width picks the format they're read in, A2.5.4: a word holds a
/// single and a doubleword a double.
pub(crate) trait Bits: Copy + Into<u128> + TryFrom<u128, Error = TryFromIntError> {
    /// The format.
    type Float: Float + FloatConvert<Self::Wide> + FloatConvert<Double> + FloatConvert<Half>;
    /// A format that holds every number of this one, scaled by up to 2^32, and the product of any
    /// two, exactly.
    type Wide: Float;
}

impl Bits for u32 {
    type Float = Single;
    type Wide = Double;
}

impl Bits for u64 {
    type Float = Double;
    type Wide = Quad;
}

/// A fixed-point format, A7.7.226: how wide it is, how many of its bits are fraction, and whether
/// it's signed.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub(crate) struct Fixed {
    pub(crate) size: u32,
    pub(crate) fraction_bits: u32,
    pub(crate) signed: bool,
}

/// `FPNeg()`, A2.5.4: the sign alone, whatever the operand is.
pub(crate) fn negate<B: Bits>(operand: B) -> B {
    bits(-float(operand))
}

/// `FPAdd()`, A2.5.4.
pub(crate) fn add<B: Bits>(op1: B, op2: B, fpscr: &mut u32) -> B {
    binary(Binary::Add, op1, op2, fpscr)
}

/// `FPSub()`.
pub(crate) fn subtract<B: Bits>(op1: B, op2: B, fpscr: &mut u32) -> B {
    binary(Binary::Subtract, op1, op2, fpscr)
}

/// `FPMul()`.
pub(crate) fn multiply<B: Bits>(op1: B, op2: B, fpscr: &mut u32) -> B {
    binary(Binary::Multiply, op1, op2, fpscr)
}

/// `FPDiv()`.
pub(crate) fn divide<B: Bits>(op1: B, op2: B, fpscr: &mut u32) -> B {
    binary(Binary::Divide, op1, op2, fpscr)
}

/// `FPMulAdd()`: `addend + op1 * op2` with one rounding.
pub(crate) fn multiply_add<B: Bits>(addend: B, op1: B, op2: B, fpscr: &mut u32) -> B {
    let addend = flush(float(addend), fpscr);
    let (op1, op2) = (flush(float(op1), fpscr), flush(float(op2), fpscr));
    let product_invalid =
        (op1.is_infinite() && op2.is_zero()) || (op1.is_zero() && op2.is_infinite());
    let nan = process_nans(&[addend, op1, op2], fpscr);
    if addend.is_nan() && !addend.is_signaling() && product_invalid {
        raise(fpscr, fpscr::IOC);
        return bits(B::Float::qnan(None));
    }
    if let Some(nan) = nan {
        return bits(nan);
    }
    let mode = rounding(*fpscr);
    let flushing = is_set(fpscr::FZ, *fpscr);
    let fused = |mode| op1.mul_add_r(op2, addend, mode);
    let wide = || {
        let product = widen::<_, B::Wide>(op1)
            .mul_r(widen(op2), Round::TowardZero)
            .value;
        product.add_r(widen(addend), Round::TowardZero).value
    };
    bits(round(fpscr, mode, flushing, fused, wide))
}

/// `FPSqrt()`.
///
/// `rustc_apfloat` has no square root, so the root is the integer one of the significand,
/// scaled to keep two or more bits past the format's precision, with what it leaves over as a
/// sticky bit below them. Rounding that to the format rounds the exact root, which is always
/// normal.
pub(crate) fn square_root<B: Bits>(operand: B, fpscr: &mut u32) -> B {
    const SCALE: ExpInt = 64;
    let operand = flush(float(operand), fpscr);
    if operand.is_nan() {
        return bits(process_nan(operand, fpscr));
    }
    if operand.is_zero() || (operand.is_infinite() && !operand.is_negative()) {
        return bits(operand);
    }
    if operand.is_negative() {
        raise(fpscr, fpscr::IOC);
        return bits(B::Float::qnan(None));
    }
    let precision = ExpInt::try_from(B::Float::PRECISION).expect("a format's precision");
    let mut exponent = operand.ilogb() - (precision - 1);
    let mut significand = operand
        .scalbn(-exponent)
        .to_u128_r(B::Float::PRECISION, Round::TowardZero, &mut false)
        .value;
    if exponent % 2 != 0 {
        significand <<= 1;
        exponent -= 1;
    }
    let scaled = significand << SCALE;
    let root = scaled.isqrt();
    let sticky = u128::from(root * root != scaled);
    let StatusAnd { status, value } = B::Float::from_u128_r(root << 1 | sticky, rounding(*fpscr));
    if status.intersects(Status::INEXACT) {
        raise(fpscr, fpscr::IXC);
    }
    bits(value.scalbn((exponent - SCALE) / 2 - 1))
}

/// `FPCompare()`, answering N, Z, C and V where `FPSCR` holds them. `quiet_nan_exc` makes a quiet
/// NaN invalid too, as VCMPE does.
pub(crate) fn compare<B: Bits>(op1: B, op2: B, quiet_nan_exc: bool, fpscr: &mut u32) -> u32 {
    let (op1, op2) = (flush(float(op1), fpscr), flush(float(op2), fpscr));
    if op1.is_nan() || op2.is_nan() {
        if op1.is_signaling() || op2.is_signaling() || quiet_nan_exc {
            raise(fpscr, fpscr::IOC);
        }
        return mask(fpscr::C) | mask(fpscr::V);
    }
    match op1.partial_cmp(&op2) {
        Some(Ordering::Equal) => mask(fpscr::Z) | mask(fpscr::C),
        Some(Ordering::Less) => mask(fpscr::N),
        _ => mask(fpscr::C),
    }
}

/// `FPMaxNum()`, A2.5.4: the larger, with a quiet NaN beside a number taken as minus infinity.
pub(crate) fn maximum_number<B: Bits>(op1: B, op2: B, fpscr: &mut u32) -> B {
    extremum_number(op1, op2, true, fpscr)
}

/// `FPMinNum()`: the smaller, with a quiet NaN beside a number taken as plus infinity.
pub(crate) fn minimum_number<B: Bits>(op1: B, op2: B, fpscr: &mut u32) -> B {
    extremum_number(op1, op2, false, fpscr)
}

/// `FPMaxNum()` or `FPMinNum()`, by `maximum`.
///
/// DDI 0403E.b calls `FPMax()` and `FPMin()` without spelling them out, so they're IEEE
/// 754-2008's maxNum and minNum, 5.3.1, on what `FPUnpack()` gives, with A2.5.4's NaN handling
/// for the NaNs left. Where 754 lets either zero win, the larger is +0 and the smaller -0. The
/// one chosen is a number already, so `FPRound()` gives it back unchanged.
fn extremum_number<B: Bits>(op1: B, op2: B, maximum: bool, fpscr: &mut u32) -> B {
    let (mut op1, mut op2) = (flush(float(op1), fpscr), flush(float(op2), fpscr));
    let quiet = |operand: B::Float| operand.is_nan() && !operand.is_signaling();
    let losing = if maximum {
        -B::Float::INFINITY
    } else {
        B::Float::INFINITY
    };
    if quiet(op1) && !quiet(op2) {
        op1 = losing;
    } else if !quiet(op1) && quiet(op2) {
        op2 = losing;
    }
    if let Some(nan) = process_nans(&[op1, op2], fpscr) {
        return bits(nan);
    }
    let first = match op1.partial_cmp(&op2) {
        Some(Ordering::Greater) => maximum,
        Some(Ordering::Less) => !maximum,
        _ => op1.is_negative() != maximum,
    };
    bits(if first { op1 } else { op2 })
}

/// `FPRoundInt()`, A2.5.4: `operand` rounded to a whole number in its own format by `mode`,
/// raising Inexact only when `exact` asks for it.
pub(crate) fn round_to_integral<B: Bits>(
    operand: B,
    mode: Round,
    exact: bool,
    fpscr: &mut u32,
) -> B {
    let operand = flush(float(operand), fpscr);
    if operand.is_nan() {
        return bits(process_nan(operand, fpscr));
    }
    let StatusAnd { status, value } = operand.round_to_integral(mode);
    if exact && status.intersects(Status::INEXACT) {
        raise(fpscr, fpscr::IXC);
    }
    bits(value)
}

/// `FPToFixed()` and `FPToFixedDirected()`: `operand` in `format`, sign- or zero-extended to 32
/// bits, rounded by `mode`.
pub(crate) fn to_fixed<B: Bits>(operand: B, format: Fixed, mode: Round, fpscr: &mut u32) -> u32 {
    let operand = flush(float(operand), fpscr);
    if operand.is_nan() {
        raise(fpscr, fpscr::IOC);
        return 0;
    }
    // The wide format holds the operand exactly, scaled by up to 2^32, so the one rounding is to
    // integer.
    let scaled = widen::<_, B::Wide>(operand).scalbn(exponent(format.fraction_bits));
    let size = format.size as usize;
    let StatusAnd { status, value } = if format.signed {
        scaled.to_i128_r(size, mode, &mut false)
    } else {
        scaled
            .to_u128_r(size, mode, &mut false)
            .map(|value| i128::try_from(value).expect("at most 32 bits"))
    };
    if status.intersects(Status::INVALID_OP) {
        raise(fpscr, fpscr::IOC);
    } else if status.intersects(Status::INEXACT) {
        raise(fpscr, fpscr::IXC);
    }
    let word = i64::try_from(value).expect("saturated to the format");
    u32::try_from(word & i64::from(u32::MAX)).expect("masked")
}

/// `FixedToFP()`: the low `format.size` bits of `operand` in the format `B` holds, rounded by
/// `mode`.
pub(crate) fn from_fixed<B: Bits>(operand: u32, format: Fixed, mode: Round, fpscr: &mut u32) -> B {
    let low = u64::from(operand) & ((1 << format.size) - 1);
    let integer = if format.signed && low >> (format.size - 1) != 0 {
        i128::from(low) - (1 << format.size)
    } else {
        i128::from(low)
    };
    // A nonzero integer of 32 bits over at most 2^32 is a normal number in either format, so
    // rounding the integer rounds the result, and the scaling after it is exact.
    let StatusAnd { status, value } = B::Float::from_i128_r(integer, mode);
    if status.intersects(Status::INEXACT) {
        raise(fpscr, fpscr::IXC);
    }
    bits(value.scalbn(-exponent(format.fraction_bits)))
}

/// `FPSingleToHalf()` and `FPDoubleToHalf()`, answering the half-precision value in the low 16
/// bits, in the format `FPSCR.AHP` picks.
pub(crate) fn to_half<B: Bits>(operand: B, fpscr: &mut u32) -> u32 {
    let operand = flush(float(operand), fpscr);
    let alternative = is_set(fpscr::AHP, *fpscr);
    let sign = if operand.is_negative() { HALF_SIGN } else { 0 };
    if operand.is_nan() {
        if operand.is_signaling() || alternative {
            raise(fpscr, fpscr::IOC);
        }
        let payload = operand.to_bits() >> (fraction_bits::<B::Float>() - HALF_FRACTION_BITS);
        let payload = u32::try_from(payload & u128::from(HALF_QUIET - 1)).expect("masked");
        return if alternative {
            sign
        } else if is_set(fpscr::DN, *fpscr) {
            HALF_DEFAULT_NAN
        } else {
            sign | HALF_DEFAULT_NAN | payload
        };
    }
    if operand.is_infinite() {
        if alternative {
            raise(fpscr, fpscr::IOC);
            return sign | HALF_ALTERNATIVE_MAX;
        }
        return sign | HALF_EXPONENT;
    }
    if operand.is_zero() {
        return sign;
    }
    let mode = rounding(*fpscr);
    if alternative && operand.ilogb() >= HALF_ALTERNATIVE_TOP {
        // Half the value rounds as the value does, one exponent down, and IEEE's overflow of it is
        // where the alternative format runs out, A2.5.4's FPRound(). A mode rounding towards zero
        // keeps IEEE's largest number without calling it overflow, so a value already out of
        // range before rounding is caught by its exponent.
        let StatusAnd { status, value } = operand.scalbn(-1).convert_r(mode, &mut false);
        if operand.ilogb() >= HALF_ALTERNATIVE_BEYOND || status.intersects(Status::OVERFLOW) {
            raise(fpscr, fpscr::IOC);
            return sign | HALF_ALTERNATIVE_MAX;
        }
        if status.intersects(Status::INEXACT) {
            raise(fpscr, fpscr::IXC);
        }
        let half: Half = value;
        return half_bits(half) + (HALF_FRACTION + 1);
    }
    let rounded = |mode| operand.convert_r(mode, &mut false);
    half_bits(round(fpscr, mode, false, rounded, || {
        widen::<_, Double>(operand)
    }))
}

/// `FPHalfToSingle()` and `FPHalfToDouble()` of the half-precision value in the low 16 bits of
/// `operand`, in the format `FPSCR.AHP` picks. It's always exact.
pub(crate) fn from_half<B: Bits>(operand: u32, fpscr: &mut u32) -> B {
    let negative = operand & HALF_SIGN != 0;
    let signed = |value: B::Float| if negative { -value } else { value };
    let field = (operand & HALF_EXPONENT) >> HALF_FRACTION_BITS;
    let fraction = operand & HALF_FRACTION;
    let ieee_top = field == HALF_EXPONENT >> HALF_FRACTION_BITS && !is_set(fpscr::AHP, *fpscr);
    if ieee_top && fraction == 0 {
        return bits(signed(B::Float::INFINITY));
    }
    if ieee_top {
        if fraction & HALF_QUIET == 0 {
            raise(fpscr, fpscr::IOC);
        }
        if is_set(fpscr::DN, *fpscr) {
            return bits(B::Float::qnan(None));
        }
        let payload = u128::from(fraction & (HALF_QUIET - 1))
            << (fraction_bits::<B::Float>() - HALF_FRACTION_BITS);
        let nan = B::Float::from_bits(B::Float::qnan(None).to_bits() | payload);
        return bits(signed(nan));
    }
    // A2.5.4: a half is 0.fraction times 2^-14 with a zero exponent field, and 1.fraction times
    // 2^(exponent - 15) otherwise, the alternative format's top exponent among them.
    let (significand, power) = if field == 0 {
        (fraction, 1 - HALF_BIAS)
    } else {
        let field = ExpInt::try_from(field).expect("five bits");
        (fraction | (HALF_FRACTION + 1), field - HALF_BIAS)
    };
    let value = B::Float::from_u128(u128::from(significand)).value;
    bits(signed(value.scalbn(power - exponent(HALF_FRACTION_BITS))))
}

/// `FPSingleToDouble()` and `FPDoubleToSingle()`, which A2.5.4 names among the operations
/// without spelling out. They're `FPHalfToDouble()` and `FPDoubleToHalf()` in IEEE format: a NaN
/// keeps its sign and the top of its fraction, and a number rounds as `FPRound()` does.
pub(crate) fn convert<A: Bits, B: Bits>(operand: A, fpscr: &mut u32) -> B
where
    A::Float: FloatConvert<B::Float>,
{
    let operand = flush(float(operand), fpscr);
    if operand.is_nan() {
        if operand.is_signaling() {
            raise(fpscr, fpscr::IOC);
        }
        if is_set(fpscr::DN, *fpscr) {
            return bits(B::Float::qnan(None));
        }
        let (from, to) = (fraction_bits::<A::Float>(), fraction_bits::<B::Float>());
        let fraction = operand.to_bits() & (quiet_bit::<A::Float>() - 1);
        let payload = if to > from {
            fraction << (to - from)
        } else {
            fraction >> (from - to)
        };
        let nan = B::Float::from_bits(B::Float::qnan(None).to_bits() | payload);
        return bits(if operand.is_negative() { -nan } else { nan });
    }
    let mode = rounding(*fpscr);
    let flushing = is_set(fpscr::FZ, *fpscr);
    let rounded = |mode| operand.convert_r(mode, &mut false);
    bits(round(fpscr, mode, flushing, rounded, || {
        widen::<_, Double>(operand)
    }))
}

/// `FPProcessNaNs()` and `FPProcessNaNs3()`, A2.5.4: the first signaling NaN, or failing that the
/// first quiet one, made quiet or the default NaN.
fn process_nans<F: Float>(operands: &[F], fpscr: &mut u32) -> Option<F> {
    let nan = operands
        .iter()
        .find(|operand| operand.is_signaling())
        .or_else(|| operands.iter().find(|operand| operand.is_nan()))?;
    Some(process_nan(*nan, fpscr))
}

/// `FPProcessNaN()`.
fn process_nan<F: Float>(operand: F, fpscr: &mut u32) -> F {
    if operand.is_signaling() {
        raise(fpscr, fpscr::IOC);
    }
    if is_set(fpscr::DN, *fpscr) {
        F::qnan(None)
    } else {
        F::from_bits(operand.to_bits() | quiet_bit::<F>())
    }
}

/// The operations of two operands `rustc_apfloat` rounds, in any of its formats.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Binary {
    Add,
    Subtract,
    Multiply,
    Divide,
}

impl Binary {
    fn apply<F: Float>(self, op1: F, op2: F, mode: Round) -> StatusAnd<F> {
        match self {
            Self::Add => op1.add_r(op2, mode),
            Self::Subtract => op1.sub_r(op2, mode),
            Self::Multiply => op1.mul_r(op2, mode),
            Self::Divide => op1.div_r(op2, mode),
        }
    }
}

/// An operation of two operands, with `FPUnpack()`'s flushing and the NaNs done first.
fn binary<B: Bits>(operation: Binary, op1: B, op2: B, fpscr: &mut u32) -> B {
    let (op1, op2) = (flush(float(op1), fpscr), flush(float(op2), fpscr));
    if let Some(nan) = process_nans(&[op1, op2], fpscr) {
        return bits(nan);
    }
    let mode = rounding(*fpscr);
    let flushing = is_set(fpscr::FZ, *fpscr);
    let rounded = |mode| operation.apply(op1, op2, mode);
    let wide = || {
        operation
            .apply(widen::<_, B::Wide>(op1), widen(op2), Round::TowardZero)
            .value
    };
    bits(round(fpscr, mode, flushing, rounded, wide))
}

/// `FPRound()`, A2.5.4, of what `operation` computes from operands that aren't NaNs. `wide` is
/// the same towards zero in a wider format, which lands on the same side as the exact value of
/// both the smallest normal number and the overflow threshold.
///
/// `rustc_apfloat` judges underflow after rounding where the architecture does before, and
/// doesn't call it overflow when a mode keeps the largest normal number. So when a result is
/// either of those, `wide` decides.
fn round<F: FloatConvert<W>, W: Float>(
    fpscr: &mut u32,
    mode: Round,
    flushing: bool,
    operation: impl Fn(Round) -> StatusAnd<F>,
    wide: impl Fn() -> W,
) -> F {
    let StatusAnd { status, value } = operation(mode);
    if value.is_nan() {
        raise(fpscr, fpscr::IOC);
        return F::qnan(None);
    }
    if status.intersects(Status::DIV_BY_ZERO) {
        raise(fpscr, fpscr::DZC);
    }
    let inexact = status.intersects(Status::INEXACT);
    let magnitude = value.abs();
    let tiny = value.is_denormal()
        || (inexact && value.is_zero())
        || (inexact
            && magnitude.bitwise_eq(F::smallest_normalized())
            && wide().abs() < widen(F::smallest_normalized()));
    if tiny && flushing {
        raise(fpscr, fpscr::UFC);
        return F::ZERO.copy_sign(value);
    }
    let beyond = W::from_u128(1).value.scalbn(F::MAX_EXP + 1);
    let overflow = status.intersects(Status::OVERFLOW)
        || (inexact && magnitude.bitwise_eq(F::largest()) && wide().abs() >= beyond);
    if overflow {
        raise(fpscr, fpscr::OFC);
    }
    if tiny && inexact {
        raise(fpscr, fpscr::UFC);
    }
    if inexact {
        raise(fpscr, fpscr::IXC);
    }
    value
}

/// `FPUnpack()`'s flush-to-zero: a denormal operand is a zero of its sign when `FPSCR.FZ` is set.
fn flush<F: Float>(operand: F, fpscr: &mut u32) -> F {
    if is_set(fpscr::FZ, *fpscr) && operand.is_denormal() {
        raise(fpscr, fpscr::IDC);
        return F::ZERO.copy_sign(operand);
    }
    operand
}

/// The rounding mode `FPSCR.RMode` selects.
pub(crate) fn rounding(fpscr: u32) -> Round {
    match get(fpscr::RMODE, fpscr) {
        0b00 => Round::NearestTiesToEven,
        0b01 => Round::TowardPositive,
        0b10 => Round::TowardNegative,
        _ => Round::TowardZero,
    }
}

fn raise(fpscr: &mut u32, exception: &FieldSpec) {
    *fpscr |= mask(exception);
}

/// The bit that makes a NaN quiet, the fraction's top, A2.5.4.
fn quiet_bit<F: Float>() -> u128 {
    F::qnan(None).to_bits() & !F::INFINITY.to_bits()
}

/// How many bits a format's fraction has.
fn fraction_bits<F: Float>() -> u32 {
    u32::try_from(F::PRECISION - 1).expect("a format's precision")
}

fn float<B: Bits>(operand: B) -> B::Float {
    B::Float::from_bits(operand.into())
}

fn bits<B: Bits>(value: B::Float) -> B {
    B::try_from(value.to_bits()).expect("the format's width")
}

fn widen<F: FloatConvert<W>, W: Float>(value: F) -> W {
    value.convert(&mut false).value
}

fn half_bits(value: Half) -> u32 {
    u32::try_from(value.to_bits()).expect("a half is 16 bits")
}

fn exponent(fraction_bits: u32) -> ExpInt {
    ExpInt::try_from(fraction_bits).expect("at most 32 fraction bits")
}
#[cfg(test)]
mod tests {
    use super::super::registers::set;
    use super::*;

    /// The masks built from the field table are A2.5.3's: the flags, and every bit there is.
    #[test]
    fn the_fpscr_masks_are_the_manuals() {
        assert_eq!(FPSCR_NZCV, 0xF000_0000);
        assert_eq!(FPSCR_WRITABLE, 0xF7C0_009F);
    }

    const INFINITY: u32 = 0x7F80_0000;
    /// `FPDefaultNaN(32)`, Table A2-2.
    const DEFAULT_NAN: u32 = 0x7FC0_0000;
    const QUIET: u32 = 0x0040_0000;
    const ONE: u32 = 0x3F80_0000;
    const TWO: u32 = 0x4000_0000;
    const THREE: u32 = 0x4040_0000;
    const MIN_NORMAL: u32 = 0x0080_0000;
    const MAX_NORMAL: u32 = 0x7F7F_FFFF;
    const SIGNALING: u32 = 0x7F80_0001;
    const QUIET_NAN: u32 = 0x7FC0_1234;
    const ROUND_UP: u32 = set(fpscr::RMODE, 0, 0b01);
    const ROUND_DOWN: u32 = set(fpscr::RMODE, 0, 0b10);
    const ROUND_ZERO: u32 = set(fpscr::RMODE, 0, 0b11);
    /// What `FPCompare()` gives an unordered pair and an equal one, NZCV `0011` and `0110`, A2.5.4,
    /// Table A2-4, read off the page so a wrong flag field shows.
    const UNORDERED: u32 = 0x3000_0000;
    const EQUAL: u32 = 0x6000_0000;

    fn with(fpscr: u32, operation: impl FnOnce(&mut u32) -> u32) -> (u32, u32) {
        let mut fpscr = fpscr;
        let result = operation(&mut fpscr);
        (result, fpscr)
    }

    #[test]
    fn the_default_nan_is_the_one_apfloat_makes() {
        assert_eq!(bits::<u32>(Single::qnan(None)), DEFAULT_NAN);
        assert_eq!(half_bits(Half::qnan(None)), HALF_DEFAULT_NAN);
    }

    /// A2.5.4's `FPAdd()`: exact sums raise nothing, and an exact zero is negative only when
    /// rounding towards minus infinity.
    #[test]
    fn add_rounds_and_signs_zero_as_the_pseudocode_does() {
        assert_eq!(with(0, |f| add(ONE, TWO, f)), (THREE, 0));
        assert_eq!(with(0, |f| add(ONE, negate(ONE), f)), (0, 0));
        assert_eq!(
            with(ROUND_DOWN, |f| add(ONE, negate(ONE), f)),
            (SIGN, ROUND_DOWN)
        );
        assert_eq!(with(0, |f| add(SIGN, SIGN, f)), (SIGN, 0));
        assert_eq!(
            with(0, |f| add(ONE, 0x3380_0000, f)),
            (ONE, mask(fpscr::IXC)),
            "1 + 2^-24 ties to even"
        );
        assert_eq!(
            with(ROUND_UP, |f| add(ONE, 0x3380_0000, f)),
            (0x3F80_0001, ROUND_UP | mask(fpscr::IXC))
        );
        assert_eq!(with(0, |f| subtract(THREE, ONE, f)), (TWO, 0));
    }

    /// Table A2-3: overflow gives infinity or the largest normal by mode, with Inexact.
    #[test]
    fn overflow_gives_the_default_result_of_each_mode() {
        let flags = mask(fpscr::OFC) | mask(fpscr::IXC);
        assert_eq!(with(0, |f| multiply(MAX_NORMAL, TWO, f)), (INFINITY, flags));
        assert_eq!(
            with(ROUND_ZERO, |f| multiply(MAX_NORMAL, TWO, f)),
            (MAX_NORMAL, ROUND_ZERO | flags)
        );
        assert_eq!(
            with(ROUND_DOWN, |f| multiply(MAX_NORMAL, TWO, f)),
            (MAX_NORMAL, ROUND_DOWN | flags)
        );
        assert_eq!(
            with(ROUND_DOWN, |f| multiply(negate(MAX_NORMAL), TWO, f)),
            (SIGN | INFINITY, ROUND_DOWN | flags)
        );
    }

    /// A2.5.4: invalid operations give the default NaN, and NaN operands hand on the first
    /// signaling one made quiet, then the first quiet one.
    #[test]
    fn nans_propagate_in_the_order_the_pseudocode_takes_them() {
        assert_eq!(
            with(0, |f| add(INFINITY, negate(INFINITY), f)),
            (DEFAULT_NAN, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| multiply(0, INFINITY, f)),
            (DEFAULT_NAN, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| divide(0, 0, f)),
            (DEFAULT_NAN, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| add(QUIET_NAN, SIGNALING, f)),
            (SIGNALING | QUIET, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| add(ONE, QUIET_NAN | SIGN, f)),
            (QUIET_NAN | SIGN, 0)
        );
        assert_eq!(with(0, |f| add(QUIET_NAN, 0x7FC0_0001, f)), (QUIET_NAN, 0));
        assert_eq!(
            with(mask(fpscr::DN), |f| add(QUIET_NAN, ONE, f)),
            (DEFAULT_NAN, mask(fpscr::DN))
        );
        assert_eq!(
            with(mask(fpscr::DN), |f| add(ONE, SIGNALING, f)),
            (DEFAULT_NAN, mask(fpscr::DN) | mask(fpscr::IOC))
        );
    }

    /// `FPDiv()`: a finite dividend over zero is Division by Zero, and an infinite one isn't.
    #[test]
    fn division_by_zero_is_raised_only_for_a_finite_dividend() {
        assert_eq!(
            with(0, |f| divide(ONE, SIGN, f)),
            (SIGN | INFINITY, mask(fpscr::DZC))
        );
        assert_eq!(with(0, |f| divide(INFINITY, 0, f)), (INFINITY, 0));
        assert_eq!(
            with(0, |f| divide(ONE, THREE, f)),
            (0x3EAA_AAAB, mask(fpscr::IXC))
        );
        assert_eq!(with(ROUND_ZERO, |f| divide(ONE, THREE, f)).0, 0x3EAA_AAAA);
    }

    /// A2.5.4's flush-to-zero: a denormal input is a zero with Input Denormal, a result tiny
    /// before rounding is a zero with Underflow and no Inexact.
    #[test]
    fn flush_to_zero_flushes_inputs_and_results() {
        let denormal = 0x0000_0001;
        assert_eq!(
            with(mask(fpscr::FZ), |f| add(denormal, ONE, f)),
            (ONE, mask(fpscr::FZ) | mask(fpscr::IDC))
        );
        assert_eq!(with(0, |f| add(denormal, denormal, f)), (2, 0));
        assert_eq!(
            with(mask(fpscr::FZ), |f| multiply(MIN_NORMAL, 0x3F00_0000, f)),
            (0, mask(fpscr::FZ) | mask(fpscr::UFC)),
            "2^-127 is tiny"
        );
        assert_eq!(
            with(mask(fpscr::FZ), |f| multiply(
                negate(MIN_NORMAL),
                0x3F00_0000,
                f
            )),
            (SIGN, mask(fpscr::FZ) | mask(fpscr::UFC)),
            "the flushed zero keeps the sign"
        );
        assert_eq!(
            with(mask(fpscr::FZ), |f| negate(denormal) ^ *f).0,
            0x8000_0001 ^ mask(fpscr::FZ)
        );
    }

    /// A2.5.4: Underflow is a result tiny before rounding and inexact. One that rounds up to the
    /// smallest normal number underflowed, and one that rounds down to it didn't.
    #[test]
    fn underflow_is_judged_before_rounding() {
        assert_eq!(
            with(0, |f| multiply(MIN_NORMAL, 0x3F00_0000, f)),
            (0x0040_0000, 0),
            "an exact denormal raises nothing"
        );
        let just_below = 0x3F7F_FFFF;
        assert_eq!(
            with(ROUND_UP, |f| multiply(MIN_NORMAL, just_below, f)),
            (MIN_NORMAL, ROUND_UP | mask(fpscr::UFC) | mask(fpscr::IXC))
        );
        assert_eq!(
            with(0, |f| multiply(0x0080_0001, just_below, f)),
            (MIN_NORMAL, mask(fpscr::IXC)),
            "the exact value is above the smallest normal"
        );
        assert_eq!(
            with(0, |f| multiply(0x0000_0003, 0x3F00_0000, f)),
            (2, mask(fpscr::UFC) | mask(fpscr::IXC)),
            "1.5 units ties to 2"
        );
    }

    /// `FPMulAdd()`: one rounding, the NaN order addend first, and a quiet NaN addend still
    /// invalid with an infinite product of zero.
    #[test]
    fn multiply_add_rounds_once() {
        let a = 0x3F80_0001;
        let b = 0x3F7F_FFFE;
        let product = multiply(a, b, &mut 0);
        assert_eq!(
            with(0, |f| multiply_add(negate(product), a, b, f)).0,
            0xA880_0000
        );
        assert_eq!(
            with(0, |f| multiply_add(ONE, TWO, THREE, f)),
            (0x40E0_0000, 0)
        );
        assert_eq!(
            with(0, |f| multiply_add(QUIET_NAN, 0, INFINITY, f)),
            (DEFAULT_NAN, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| multiply_add(QUIET_NAN, SIGNALING, ONE, f)),
            (SIGNALING | QUIET, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| multiply_add(ONE, QUIET_NAN, ONE, f)),
            (QUIET_NAN, 0)
        );
        assert_eq!(
            with(0, |f| multiply_add(0, SIGN, ONE, f)),
            (0, 0),
            "+0 + -0 is +0"
        );
        assert_eq!(with(0, |f| multiply_add(SIGN, SIGN, ONE, f)), (SIGN, 0));
        assert_eq!(
            with(0, |f| multiply_add(INFINITY, negate(INFINITY), ONE, f)),
            (DEFAULT_NAN, mask(fpscr::IOC))
        );
    }

    /// `FPSqrt()`: exact roots, rounded ones in every mode, and the special cases.
    #[test]
    fn square_root_is_correctly_rounded() {
        assert_eq!(
            with(0, |f| square_root(0x4110_0000_u32, f)),
            (THREE, 0),
            "of 9"
        );
        assert_eq!(
            with(0, |f| square_root(TWO, f)),
            (0x3FB5_04F3, mask(fpscr::IXC))
        );
        assert_eq!(with(ROUND_ZERO, |f| square_root(TWO, f)).0, 0x3FB5_04F3);
        assert_eq!(with(ROUND_UP, |f| square_root(TWO, f)).0, 0x3FB5_04F4);
        assert_eq!(with(0, |f| square_root(0x4000_0001_u32, f)).0, 0x3FB5_04F4);
        assert_eq!(
            with(0, |f| square_root(1_u32, f)),
            (0x1A35_04F3, mask(fpscr::IXC)),
            "of 2^-149"
        );
        assert_eq!(
            with(0, |f| square_root(0x0020_0000_u32, f)),
            (0x1F80_0000, 0),
            "of 2^-128"
        );
        assert_eq!(
            with(0, |f| square_root(0x0040_0000_u32, f)),
            (0x1FB5_04F3, mask(fpscr::IXC)),
            "of 2^-127"
        );
        assert_eq!(with(0, |f| square_root(SIGN, f)), (SIGN, 0));
        assert_eq!(
            with(0, |f| square_root(negate(ONE), f)),
            (DEFAULT_NAN, mask(fpscr::IOC))
        );
        assert_eq!(with(0, |f| square_root(INFINITY, f)), (INFINITY, 0));
        assert_eq!(
            with(0, |f| square_root(SIGNALING, f)),
            (SIGNALING | QUIET, mask(fpscr::IOC))
        );
        assert_eq!(
            with(mask(fpscr::FZ), |f| square_root(1_u32, f)),
            (0, mask(fpscr::FZ) | mask(fpscr::IDC))
        );
    }

    /// `FPCompare()`: Table A2-4's flags, an unordered pair, and which NaNs are invalid.
    #[test]
    fn compare_gives_table_a2_4() {
        assert_eq!(with(0, |f| compare(ONE, ONE, false, f)), (EQUAL, 0));
        assert_eq!(with(0, |f| compare(0_u32, SIGN, false, f)), (EQUAL, 0));
        assert_eq!(
            with(0, |f| compare(ONE, TWO, false, f)),
            (mask(fpscr::N), 0)
        );
        assert_eq!(
            with(0, |f| compare(TWO, ONE, false, f)),
            (mask(fpscr::C), 0)
        );
        assert_eq!(
            with(0, |f| compare(QUIET_NAN, ONE, false, f)),
            (UNORDERED, 0)
        );
        assert_eq!(
            with(0, |f| compare(QUIET_NAN, ONE, true, f)),
            (UNORDERED, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| compare(ONE, SIGNALING, false, f)),
            (UNORDERED, mask(fpscr::IOC))
        );
        assert_eq!(
            with(mask(fpscr::FZ), |f| compare(1_u32, 0, false, f)),
            (EQUAL, mask(fpscr::FZ) | mask(fpscr::IDC))
        );
    }

    /// `FPToFixed()`: rounding, saturation with Invalid Operation, and a NaN as zero.
    #[test]
    fn to_fixed_rounds_and_saturates() {
        let int = |signed| Fixed {
            size: 32,
            fraction_bits: 0,
            signed,
        };
        let two_and_a_half: u32 = 0x4020_0000;
        assert_eq!(
            with(0, |f| to_fixed(
                two_and_a_half,
                int(true),
                Round::TowardZero,
                f
            )),
            (2, mask(fpscr::IXC))
        );
        assert_eq!(
            with(0, |f| to_fixed(two_and_a_half, int(true), rounding(*f), f)),
            (2, mask(fpscr::IXC))
        );
        assert_eq!(
            with(ROUND_UP, |f| to_fixed(
                two_and_a_half,
                int(true),
                rounding(*f),
                f
            ))
            .0,
            3
        );
        assert_eq!(
            with(0, |f| to_fixed(
                negate(two_and_a_half),
                int(true),
                Round::TowardZero,
                f
            )),
            (0xFFFF_FFFE, mask(fpscr::IXC))
        );
        assert_eq!(
            with(0, |f| to_fixed(
                negate(ONE),
                int(false),
                Round::TowardZero,
                f
            )),
            (0, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| to_fixed(
                0xBF00_0000_u32,
                int(false),
                Round::TowardZero,
                f
            )),
            (0, mask(fpscr::IXC)),
            "-0.5 rounds to zero before saturating"
        );
        assert_eq!(
            with(0, |f| to_fixed(
                0x4F00_0000_u32,
                int(true),
                Round::TowardZero,
                f
            )),
            (0x7FFF_FFFF, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| to_fixed(
                negate(INFINITY),
                int(true),
                Round::TowardZero,
                f
            )),
            (0x8000_0000, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| to_fixed(
                0x4F00_0000_u32,
                int(false),
                Round::TowardZero,
                f
            )),
            (0x8000_0000, 0)
        );
        assert_eq!(
            with(0, |f| to_fixed(QUIET_NAN, int(true), Round::TowardZero, f)),
            (0, mask(fpscr::IOC))
        );
        let q15 = Fixed {
            size: 16,
            fraction_bits: 15,
            signed: true,
        };
        assert_eq!(
            with(0, |f| to_fixed(
                negate(0x3F00_0000_u32),
                q15,
                Round::TowardZero,
                f
            )),
            (0xFFFF_C000, 0),
            "-0.5 in Q15, sign-extended"
        );
        assert_eq!(
            with(0, |f| to_fixed(ONE, q15, Round::TowardZero, f)),
            (0x7FFF, mask(fpscr::IOC))
        );
    }

    /// `FixedToFP()`: the low bits as signed or unsigned, scaled, rounded once.
    #[test]
    fn from_fixed_rounds_once() {
        let int = |signed| Fixed {
            size: 32,
            fraction_bits: 0,
            signed,
        };
        assert_eq!(
            with(0, |f| from_fixed(3, int(true), rounding(*f), f)),
            (THREE, 0)
        );
        assert_eq!(
            with(0, |f| from_fixed(u32::MAX, int(true), rounding(*f), f)),
            (negate(ONE), 0)
        );
        assert_eq!(
            with(0, |f| from_fixed(u32::MAX, int(false), rounding(*f), f)),
            (0x4F80_0000, mask(fpscr::IXC))
        );
        assert_eq!(
            with(ROUND_ZERO, |f| from_fixed(
                u32::MAX,
                int(false),
                rounding(*f),
                f
            )),
            (0x4F7F_FFFF, ROUND_ZERO | mask(fpscr::IXC))
        );
        assert_eq!(
            with(ROUND_ZERO, |f| from_fixed(
                u32::MAX,
                int(false),
                Round::NearestTiesToEven,
                f
            ))
            .0,
            0x4F80_0000,
            "the fixed-point form rounds to nearest whatever FPSCR says"
        );
        let q15 = Fixed {
            size: 16,
            fraction_bits: 15,
            signed: true,
        };
        assert_eq!(
            with(0, |f| from_fixed(
                0xABCD_C000,
                q15,
                Round::NearestTiesToEven,
                f
            )),
            (negate(0x3F00_0000_u32), 0)
        );
        assert_eq!(
            with(0, |f| from_fixed(0, q15, Round::NearestTiesToEven, f)),
            (0, 0)
        );
    }

    /// `FPSingleToHalf()` and `FPHalfToSingle()` in both formats, A2.5.4.
    #[test]
    fn half_precision_converts_in_either_format() {
        assert_eq!(with(0, |f| to_half(ONE, f)), (0x3C00, 0));
        assert_eq!(with(0, |f| from_half(0x3C00, f)), (ONE, 0));
        assert_eq!(
            with(0, |f| to_half(0x477F_E000_u32, f)),
            (0x7BFF, 0),
            "65504"
        );
        assert_eq!(
            with(0, |f| to_half(0x477F_F000_u32, f)),
            (0x7C00, mask(fpscr::OFC) | mask(fpscr::IXC)),
            "65520 ties to 65536"
        );
        assert_eq!(
            with(ROUND_ZERO, |f| to_half(0x477F_F000_u32, f)),
            (0x7BFF, ROUND_ZERO | mask(fpscr::IXC)),
            "towards zero it's the largest, and in range"
        );
        assert_eq!(
            with(ROUND_ZERO, |f| to_half(0x4780_0000_u32, f)),
            (0x7BFF, ROUND_ZERO | mask(fpscr::OFC) | mask(fpscr::IXC))
        );
        assert_eq!(with(0, |f| to_half(0x3380_0000_u32, f)), (1, 0), "2^-24");
        assert_eq!(
            with(0, |f| to_half(0x3300_0000_u32, f)),
            (0, mask(fpscr::UFC) | mask(fpscr::IXC))
        );
        assert_eq!(with(0, |f| to_half(QUIET_NAN, f)), (0x7E00, 0));
        assert_eq!(
            with(0, |f| to_half(SIGNALING | SIGN, f)),
            (0xFE00, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| from_half(0x7D00, f)),
            (0x7FE0_0000, mask(fpscr::IOC))
        );
        assert_eq!(with(0, |f| from_half(0x0001, f)), (0x3380_0000, 0));
        assert_eq!(with(0, |f| from_half(0xFC00, f)), (SIGN | INFINITY, 0));

        let ahp = mask(fpscr::AHP);
        assert_eq!(
            with(ahp, |f| to_half(0x4780_0000_u32, f)),
            (0x7C00, ahp),
            "65536"
        );
        assert_eq!(with(ahp, |f| from_half(0x7C00, f)), (0x4780_0000, ahp));
        assert_eq!(
            with(ahp, |f| to_half(0x47FF_E000_u32, f)),
            (0x7FFF, ahp),
            "131008"
        );
        assert_eq!(
            with(ahp, |f| to_half(0x4800_0000_u32, f)),
            (HALF_ALTERNATIVE_MAX, ahp | mask(fpscr::IOC)),
            "2^17 is too big and is Invalid, not Overflow"
        );
        assert_eq!(
            with(ahp, |f| to_half(INFINITY | SIGN, f)),
            (0xFFFF, ahp | mask(fpscr::IOC))
        );
        assert_eq!(
            with(ahp, |f| to_half(QUIET_NAN | SIGN, f)),
            (0x8000, ahp | mask(fpscr::IOC))
        );
        assert_eq!(
            with(ahp, |f| to_half(0x47FF_E800_u32, f)),
            (0x7FFF, ahp | mask(fpscr::IXC)),
            "131024"
        );
        assert_eq!(
            with(ahp, |f| to_half(0x47FF_F000_u32, f)),
            (0x7FFF, ahp | mask(fpscr::IOC)),
            "131040 ties to 2^17"
        );
    }

    /// A2.5.4's `FPRound()` in the alternative format: 2^17 or more is Invalid with no Inexact in
    /// every mode, even one that rounds it towards zero, and just below 2^17 is the largest number.
    #[test]
    fn the_alternative_format_runs_out_at_2_17_in_every_mode() {
        let ahp = mask(fpscr::AHP);
        for (name, mode, operand, expected) in [
            ("2^17 towards zero", ROUND_ZERO, 0x4800_0000_u32, 0x7FFF),
            ("2^18 towards minus", ROUND_DOWN, 0x4880_0000, 0x7FFF),
            ("-2^17 towards plus", ROUND_UP, 0xC800_0000, 0xFFFF),
            ("-2^17 towards zero", ROUND_ZERO, 0xC800_0000, 0xFFFF),
            ("2^17 towards plus", ROUND_UP, 0x4800_0000, 0x7FFF),
            (
                "just above 131008 rounds up to 2^17",
                ROUND_UP,
                0x47FF_E001,
                0x7FFF,
            ),
        ] {
            assert_eq!(
                with(ahp | mode, |f| to_half(operand, f)),
                (expected, ahp | mode | mask(fpscr::IOC)),
                "{name}"
            );
        }
        assert_eq!(
            with(ahp | ROUND_ZERO, |f| to_half(0x47FF_FFFF_u32, f)),
            (0x7FFF, ahp | ROUND_ZERO | mask(fpscr::IXC)),
            "just below 2^17 towards zero is in range"
        );
    }

    const D_SIGN: u64 = 0x8000_0000_0000_0000;
    const D_ONE: u64 = 0x3FF0_0000_0000_0000;
    const D_TWO: u64 = 0x4000_0000_0000_0000;
    const D_THREE: u64 = 0x4008_0000_0000_0000;
    const D_HALF: u64 = 0x3FE0_0000_0000_0000;
    const D_MIN_NORMAL: u64 = 0x0010_0000_0000_0000;
    const D_MAX_NORMAL: u64 = 0x7FEF_FFFF_FFFF_FFFF;
    const D_INFINITY: u64 = 0x7FF0_0000_0000_0000;
    /// `FPDefaultNaN(64)`, Table A2-2.
    const D_DEFAULT_NAN: u64 = 0x7FF8_0000_0000_0000;
    const D_QUIET: u64 = 0x0008_0000_0000_0000;
    const D_SIGNALING: u64 = 0x7FF0_0000_0000_0001;
    const D_QUIET_NAN: u64 = 0x7FF8_0000_0000_1234;
    /// 1 - 2^-53, the double just below 1.
    const D_JUST_BELOW_ONE: u64 = 0x3FEF_FFFF_FFFF_FFFF;
    const D_TWO_AND_A_HALF: u64 = 0x4004_0000_0000_0000;
    /// 1e10, past every 32-bit integer.
    const D_TEN_BILLION: u64 = 0x4202_A05F_2000_0000;

    fn with64(fpscr: u32, operation: impl FnOnce(&mut u32) -> u64) -> (u64, u32) {
        let mut fpscr = fpscr;
        let result = operation(&mut fpscr);
        (result, fpscr)
    }

    #[test]
    fn the_double_default_nan_is_the_one_apfloat_makes() {
        assert_eq!(bits::<u64>(Double::qnan(None)), D_DEFAULT_NAN);
    }

    /// A2.5.4's pseudocode takes N as 32 or 64, so a double rounds, signs its zeros, overflows,
    /// underflows before rounding and hands on its NaNs as a single does.
    #[test]
    fn double_precision_rounds_and_raises_as_single_precision_does() {
        assert_eq!(with64(0, |f| add(D_ONE, D_TWO, f)), (D_THREE, 0));
        let tiny = 0x3CA0_0000_0000_0000;
        assert_eq!(
            with64(0, |f| add(D_ONE, tiny, f)),
            (D_ONE, mask(fpscr::IXC)),
            "1 + 2^-53 ties to even"
        );
        assert_eq!(
            with64(ROUND_UP, |f| add(D_ONE, tiny, f)),
            (D_ONE + 1, ROUND_UP | mask(fpscr::IXC))
        );
        assert_eq!(
            with64(ROUND_DOWN, |f| subtract(D_ONE, D_ONE, f)),
            (D_SIGN, ROUND_DOWN)
        );
        assert_eq!(
            with64(0, |f| multiply(D_MAX_NORMAL, D_TWO, f)),
            (D_INFINITY, mask(fpscr::OFC) | mask(fpscr::IXC))
        );
        assert_eq!(
            with64(ROUND_ZERO, |f| multiply(D_MAX_NORMAL, D_TWO, f)),
            (
                D_MAX_NORMAL,
                ROUND_ZERO | mask(fpscr::OFC) | mask(fpscr::IXC)
            )
        );
        assert_eq!(
            with64(ROUND_UP, |f| multiply(D_MIN_NORMAL, D_JUST_BELOW_ONE, f)),
            (D_MIN_NORMAL, ROUND_UP | mask(fpscr::UFC) | mask(fpscr::IXC)),
            "tiny before rounding up to the smallest normal"
        );
        assert_eq!(
            with64(0, |f| multiply(D_MIN_NORMAL + 1, D_JUST_BELOW_ONE, f)),
            (D_MIN_NORMAL, mask(fpscr::IXC)),
            "the exact value is above the smallest normal"
        );
        assert_eq!(
            with64(0, |f| divide(D_ONE, D_SIGN, f)),
            (D_SIGN | D_INFINITY, mask(fpscr::DZC))
        );
        assert_eq!(
            with64(0, |f| divide(D_ONE, D_THREE, f)),
            (0x3FD5_5555_5555_5555, mask(fpscr::IXC))
        );
    }

    /// A2.5.4's NaN handling and flush-to-zero for doubles: the bit that quiets is bit 51, and FZ
    /// flushes double denormals in and out.
    #[test]
    fn double_nans_and_denormals_are_handled_as_single_ones_are() {
        assert_eq!(
            with64(0, |f| add(D_QUIET_NAN, D_SIGNALING, f)),
            (D_SIGNALING | D_QUIET, mask(fpscr::IOC))
        );
        assert_eq!(with64(0, |f| add(D_ONE, D_QUIET_NAN, f)), (D_QUIET_NAN, 0));
        assert_eq!(
            with64(mask(fpscr::DN), |f| add(D_QUIET_NAN, D_ONE, f)),
            (D_DEFAULT_NAN, mask(fpscr::DN))
        );
        assert_eq!(
            with64(0, |f| multiply(0, D_INFINITY, f)),
            (D_DEFAULT_NAN, mask(fpscr::IOC))
        );
        assert_eq!(
            with64(mask(fpscr::FZ), |f| add(1, D_ONE, f)),
            (D_ONE, mask(fpscr::FZ) | mask(fpscr::IDC))
        );
        assert_eq!(with64(0, |f| add(1, 1, f)), (2, 0), "denormals add exactly");
        assert_eq!(
            with64(mask(fpscr::FZ), |f| multiply(
                negate(D_MIN_NORMAL),
                D_HALF,
                f
            )),
            (D_SIGN, mask(fpscr::FZ) | mask(fpscr::UFC))
        );
        assert_eq!(
            negate(1_u64),
            D_SIGN | 1,
            "FPNeg() of a denormal flushes nothing"
        );
    }

    /// `FPMulAdd()` and `FPSqrt()` round once in double precision too. (1 + 2^-52)(1 - 2^-52) is
    /// 1 - 2^-104, which a product rounds to 1, and the root of 2 is above its double.
    #[test]
    fn double_fused_multiply_add_and_square_root_round_once() {
        let (a, b) = (D_ONE + 1, 0x3FEF_FFFF_FFFF_FFFE);
        let product = multiply(a, b, &mut 0);
        assert_eq!(product, D_ONE);
        assert_eq!(
            with64(0, |f| multiply_add(negate(product), a, b, f)),
            (0xB970_0000_0000_0000, 0),
            "-2^-104"
        );
        let nine = 0x4022_0000_0000_0000;
        assert_eq!(with64(0, |f| square_root(nine, f)), (D_THREE, 0));
        assert_eq!(
            with64(0, |f| square_root(D_TWO, f)),
            (0x3FF6_A09E_667F_3BCD, mask(fpscr::IXC))
        );
        assert_eq!(
            with64(ROUND_ZERO, |f| square_root(D_TWO, f)).0,
            0x3FF6_A09E_667F_3BCC
        );
        assert_eq!(
            with64(0, |f| square_root(1, f)),
            (0x1E60_0000_0000_0000, 0),
            "of 2^-1074 is 2^-537"
        );
        assert_eq!(
            with64(0, |f| square_root(negate(D_ONE), f)),
            (D_DEFAULT_NAN, mask(fpscr::IOC))
        );
        assert_eq!(
            with64(mask(fpscr::FZ), |f| square_root(1, f)),
            (0, mask(fpscr::FZ) | mask(fpscr::IDC))
        );
    }

    /// `FPCompare()` on doubles, Table A2-4.
    #[test]
    fn double_compare_gives_table_a2_4() {
        assert_eq!(
            with(0, |f| compare(D_ONE, D_TWO, false, f)),
            (mask(fpscr::N), 0)
        );
        assert_eq!(with(0, |f| compare(0, D_SIGN, false, f)), (EQUAL, 0));
        assert_eq!(
            with(0, |f| compare(D_QUIET_NAN, D_ONE, true, f)),
            (UNORDERED, mask(fpscr::IOC))
        );
        assert_eq!(
            with(mask(fpscr::FZ), |f| compare(1_u64, 0, false, f)),
            (EQUAL, mask(fpscr::FZ) | mask(fpscr::IDC))
        );
    }

    /// `FPToFixed()` and `FPToFixedDirected()` from a double, in every mode, and `FixedToFP()` to
    /// one, which a 32-bit integer always is exactly.
    #[test]
    fn doubles_convert_to_and_from_fixed_point() {
        let int = |signed| Fixed {
            size: 32,
            fraction_bits: 0,
            signed,
        };
        for (mode, expected) in [
            (Round::TowardZero, 2),
            (Round::TowardPositive, 3),
            (Round::NearestTiesToEven, 2),
            (Round::NearestTiesToAway, 3),
            (Round::TowardNegative, 2),
        ] {
            assert_eq!(
                with(0, |f| to_fixed(D_TWO_AND_A_HALF, int(true), mode, f)),
                (expected, mask(fpscr::IXC)),
                "2.5 by {mode:?}"
            );
        }
        assert_eq!(
            with(0, |f| to_fixed(
                negate(D_TWO_AND_A_HALF),
                int(true),
                Round::NearestTiesToAway,
                f
            )),
            (0xFFFF_FFFD, mask(fpscr::IXC)),
            "-2.5 away from zero"
        );
        assert_eq!(
            with(0, |f| to_fixed(
                D_TEN_BILLION,
                int(true),
                Round::TowardZero,
                f
            )),
            (0x7FFF_FFFF, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, |f| to_fixed(
                negate(D_HALF),
                int(false),
                Round::TowardZero,
                f
            )),
            (0, mask(fpscr::IXC))
        );
        assert_eq!(
            with(0, |f| to_fixed(
                negate(D_HALF),
                int(false),
                Round::TowardNegative,
                f
            )),
            (0, mask(fpscr::IOC)),
            "-0.5 rounds to -1, which saturates"
        );
        assert_eq!(
            with(0, |f| to_fixed(1_u64, int(true), Round::TowardPositive, f)),
            (1, mask(fpscr::IXC))
        );
        assert_eq!(
            with(mask(fpscr::FZ), |f| to_fixed(
                1_u64,
                int(true),
                Round::TowardPositive,
                f
            )),
            (0, mask(fpscr::FZ) | mask(fpscr::IDC))
        );
        assert_eq!(
            with64(0, |f| from_fixed(
                u32::MAX,
                int(false),
                Round::NearestTiesToEven,
                f
            )),
            (0x41EF_FFFF_FFE0_0000, 0),
            "2^32 - 1 is exact"
        );
        assert_eq!(
            with64(0, |f| from_fixed(u32::MAX, int(true), Round::TowardZero, f)),
            (negate(D_ONE), 0)
        );
        let fraction = Fixed {
            size: 32,
            fraction_bits: 32,
            signed: false,
        };
        assert_eq!(
            with64(0, |f| from_fixed(
                u32::MAX,
                fraction,
                Round::NearestTiesToEven,
                f
            )),
            (0x3FEF_FFFF_FFE0_0000, 0),
            "1 - 2^-32 is exact"
        );
    }

    /// `FPDoubleToHalf()` rounds once, with no single in between, and `FPHalfToDouble()` is exact,
    /// NaNs keeping the top of their fraction, A2.5.4.
    #[test]
    fn doubles_convert_to_and_from_half_precision() {
        assert_eq!(with(0, |f| to_half(D_ONE, f)), (0x3C00, 0));
        assert_eq!(
            with(0, |f| to_half(0x3FF0_0200_0000_1000_u64, f)),
            (0x3C01, mask(fpscr::IXC)),
            "1 + 2^-11 + 2^-40 rounds up, where a single on the way would tie to even"
        );
        assert_eq!(
            with(0, |f| to_half(0x40EF_FE00_0000_0000_u64, f)),
            (0x7C00, mask(fpscr::OFC) | mask(fpscr::IXC)),
            "65520"
        );
        assert_eq!(
            with(mask(fpscr::AHP), |f| to_half(0x4100_0000_0000_0000_u64, f)),
            (HALF_ALTERNATIVE_MAX, mask(fpscr::AHP) | mask(fpscr::IOC)),
            "2^17"
        );
        assert_eq!(
            with(0, |f| to_half(0x7FF8_0400_0000_0000_u64, f)),
            (0x7E01, 0),
            "bits 50 to 42 of the fraction"
        );
        assert_eq!(with64(0, |f| from_half(0x3C00, f)), (D_ONE, 0));
        assert_eq!(
            with64(0, |f| from_half(0x0001, f)),
            (0x3E70_0000_0000_0000, 0),
            "2^-24"
        );
        assert_eq!(
            with64(0, |f| from_half(0xFC00, f)),
            (D_SIGN | D_INFINITY, 0)
        );
        assert_eq!(
            with64(0, |f| from_half(0x7D00, f)),
            (0x7FFC_0000_0000_0000, mask(fpscr::IOC)),
            "quieted, with the fraction's bits 8 to 0 at 50 to 42"
        );
        assert_eq!(
            with64(mask(fpscr::AHP), |f| from_half(0x7C00, f)),
            (0x40F0_0000_0000_0000, mask(fpscr::AHP)),
            "the alternative format's 2^16"
        );
    }

    /// `FPSingleToDouble()` is exact, and `FPDoubleToSingle()` rounds as `FPRound()` does, with
    /// FZ and its exceptions. A NaN keeps its sign and the top of its fraction.
    #[test]
    fn singles_and_doubles_convert_as_fpround_has_it() {
        assert_eq!(
            with64(0, |f| convert::<u32, u64>(0x3EAA_AAAB, f)),
            (0x3FD5_5555_6000_0000, 0)
        );
        assert_eq!(
            with(0, |f| convert::<u64, u32>(0x3FD5_5555_5555_5555, f)),
            (0x3EAA_AAAB, mask(fpscr::IXC)),
            "a third"
        );
        let two_to_minus_130 = 0x37D0_0000_0000_0000;
        assert_eq!(
            with(0, |f| convert::<u64, u32>(two_to_minus_130, f)),
            (0x0008_0000, 0),
            "an exact denormal"
        );
        assert_eq!(
            with(mask(fpscr::FZ), |f| convert::<u64, u32>(
                two_to_minus_130,
                f
            )),
            (0, mask(fpscr::FZ) | mask(fpscr::UFC))
        );
        assert_eq!(
            with(0, |f| convert::<u64, u32>(D_MAX_NORMAL, f)),
            (INFINITY, mask(fpscr::OFC) | mask(fpscr::IXC))
        );
        assert_eq!(
            with(ROUND_ZERO, |f| convert::<u64, u32>(D_MAX_NORMAL, f)),
            (MAX_NORMAL, ROUND_ZERO | mask(fpscr::OFC) | mask(fpscr::IXC))
        );
        assert_eq!(
            with(0, |f| convert::<u64, u32>(0xFFF4_0000_0000_0000, f)),
            (0xFFE0_0000, mask(fpscr::IOC))
        );
        assert_eq!(
            with64(0, |f| convert::<u32, u64>(SIGNALING, f)),
            (0x7FF8_0000_2000_0000, mask(fpscr::IOC))
        );
        assert_eq!(
            with64(mask(fpscr::DN), |f| convert::<u32, u64>(QUIET_NAN, f)),
            (D_DEFAULT_NAN, mask(fpscr::DN))
        );
        assert_eq!(
            with64(0, |f| convert::<u32, u64>(1, f)),
            (0x36A0_0000_0000_0000, 0),
            "2^-149"
        );
        assert_eq!(
            with64(mask(fpscr::FZ), |f| convert::<u32, u64>(1, f)),
            (0, mask(fpscr::FZ) | mask(fpscr::IDC))
        );
    }

    /// `FPRoundInt()`, A2.5.4, in every mode and with ties away: a zero result keeps the operand's
    /// sign, and only an exact rounding raises Inexact.
    #[test]
    fn round_to_integral_rounds_by_its_mode() {
        let (two_and_a_half, minus_two_and_a_half) = (0x4020_0000, 0xC020_0000);
        let (half, minus_half) = (0x3F00_0000, 0xBF00_0000);
        let away = Round::NearestTiesToAway;
        let even = Round::NearestTiesToEven;
        let up = Round::TowardPositive;
        let down = Round::TowardNegative;
        let zero = Round::TowardZero;
        #[rustfmt::skip]
        let cases = [
            ("2.5 away", two_and_a_half, away, false, (THREE, 0)),
            ("-2.5 away", minus_two_and_a_half, away, false, (0xC040_0000, 0)),
            ("-0.4 away", 0xBECC_CCCD, away, false, (SIGN, 0)),
            ("2.5 even", two_and_a_half, even, false, (TWO, 0)),
            ("3.5 even", 0x4060_0000, even, false, (0x4080_0000, 0)),
            ("-0.5 even", minus_half, even, false, (SIGN, 0)),
            ("2.1 up", 0x4006_6666, up, false, (THREE, 0)),
            ("-2.5 up", minus_two_and_a_half, up, false, (0xC000_0000, 0)),
            ("-0.5 up", minus_half, up, false, (SIGN, 0)),
            ("2.5 down", two_and_a_half, down, false, (TWO, 0)),
            ("-2.5 down", minus_two_and_a_half, down, false, (0xC040_0000, 0)),
            ("0.5 down", half, down, false, (0, 0)),
            ("-2.5 to zero", minus_two_and_a_half, zero, false, (0xC000_0000, 0)),
            ("2.9 to zero", 0x4039_999A, zero, false, (TWO, 0)),
            ("2.5 exact", two_and_a_half, even, true, (TWO, mask(fpscr::IXC))),
            ("2 exact", TWO, even, true, (TWO, 0)),
            ("2^30", 0x4E80_0000, even, true, (0x4E80_0000, 0)),
            ("infinity", INFINITY | SIGN, up, true, (INFINITY | SIGN, 0)),
            ("-0", SIGN, up, true, (SIGN, 0)),
            ("a signaling NaN", SIGNALING, even, false, (SIGNALING | QUIET, mask(fpscr::IOC))),
            ("2^-149 up", 1, up, false, (ONE, 0)),
        ];
        for (name, operand, mode, exact, expected) in cases {
            assert_eq!(
                with(0, |f| round_to_integral(operand, mode, exact, f)),
                expected,
                "{name}"
            );
        }
        assert_eq!(
            with(mask(fpscr::FZ), |f| round_to_integral(1_u32, up, false, f)),
            (0, mask(fpscr::FZ) | mask(fpscr::IDC)),
            "a flushed denormal is a zero"
        );
        let below_2_52 = 0x432F_FFFF_FFFF_FFFF;
        assert_eq!(
            with64(0, |f| round_to_integral(below_2_52, even, true, f)),
            (0x4330_0000_0000_0000, mask(fpscr::IXC)),
            "2^52 - 0.5 ties to the even 2^52"
        );
        assert_eq!(
            with64(0, |f| round_to_integral(below_2_52, zero, false, f)),
            (0x432F_FFFF_FFFF_FFFE, 0)
        );
        assert_eq!(
            with64(0, |f| round_to_integral(
                negate(D_TWO_AND_A_HALF),
                down,
                false,
                f
            )),
            (negate(D_THREE), 0)
        );
    }

    /// `FPMaxNum()` and `FPMinNum()`, A2.5.4: a quiet NaN beside a number loses, a signaling one
    /// is Invalid, +0 is the larger zero, and FZ flushes the operands.
    #[test]
    fn maximum_and_minimum_number_take_a_number_over_a_quiet_nan() {
        let max = |op1, op2| move |f: &mut u32| maximum_number(op1, op2, f);
        let min = |op1, op2| move |f: &mut u32| minimum_number(op1, op2, f);
        assert_eq!(with(0, max(ONE, TWO)), (TWO, 0));
        assert_eq!(with(0, max(TWO, ONE)), (TWO, 0));
        assert_eq!(with(0, min(TWO, ONE)), (ONE, 0));
        assert_eq!(with(0, max(0, SIGN)), (0, 0));
        assert_eq!(with(0, max(SIGN, 0)), (0, 0));
        assert_eq!(with(0, min(0, SIGN)), (SIGN, 0));
        assert_eq!(with(0, min(SIGN, 0)), (SIGN, 0));
        assert_eq!(with(0, max(QUIET_NAN, ONE)), (ONE, 0));
        assert_eq!(with(0, min(ONE, QUIET_NAN)), (ONE, 0));
        assert_eq!(
            with(0, max(QUIET_NAN, INFINITY | SIGN)),
            (INFINITY | SIGN, 0)
        );
        assert_eq!(with(0, max(QUIET_NAN, 0x7FC0_0001)), (QUIET_NAN, 0));
        assert_eq!(
            with(0, max(SIGNALING, ONE)),
            (SIGNALING | QUIET, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, max(ONE, SIGNALING)),
            (SIGNALING | QUIET, mask(fpscr::IOC))
        );
        assert_eq!(
            with(0, min(QUIET_NAN, SIGNALING)),
            (SIGNALING | QUIET, mask(fpscr::IOC))
        );
        assert_eq!(
            with(mask(fpscr::DN), max(QUIET_NAN, QUIET_NAN)),
            (DEFAULT_NAN, mask(fpscr::DN))
        );
        assert_eq!(with(0, max(1, 0)), (1, 0), "a denormal is a number");
        assert_eq!(
            with(mask(fpscr::FZ), max(SIGN | 1, 0)),
            (0, mask(fpscr::FZ) | mask(fpscr::IDC)),
            "flushed to -0, which loses to +0"
        );
        assert_eq!(with64(0, |f| maximum_number(D_ONE, D_TWO, f)), (D_TWO, 0));
        assert_eq!(
            with64(0, |f| minimum_number(D_QUIET_NAN, D_ONE, f)),
            (D_ONE, 0)
        );
        assert_eq!(
            with64(0, |f| maximum_number(D_SIGNALING, D_ONE, f)),
            (D_SIGNALING | D_QUIET, mask(fpscr::IOC))
        );
    }
}