use super::{Env, Flags, MinMax, Propagate, Round, Spec, Tininess};
const EXTRA: u32 = 3;
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(super) enum Class {
Zero,
Finite,
Inf,
Nan,
}
#[derive(Debug, Clone, Copy)]
pub(super) struct Parts {
pub sign: bool,
pub exp: i32,
pub frac: u64,
pub class: Class,
pub snan: bool,
}
impl Parts {
pub(super) const fn zero(sign: bool) -> Parts {
Parts {
sign,
exp: 0,
frac: 0,
class: Class::Zero,
snan: false,
}
}
fn normalized(mut self) -> Parts {
if self.class == Class::Finite && self.frac != 0 {
let n = self.frac.leading_zeros();
self.frac <<= n;
self.exp -= n as i32;
}
self
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(super) enum Rounded {
Zero,
Finite { exp: i32, frac: u64 },
Inf,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(super) enum Outcome {
Num(bool, Rounded),
DefaultNan,
Nan { sign: bool, payload: u64 },
}
pub(super) fn shr_sticky(v: u128, n: u32) -> u128 {
if n == 0 {
return v;
}
if n >= 128 {
return u128::from(v != 0);
}
let lost = v & ((1u128 << n) - 1);
(v >> n) | u128::from(lost != 0)
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub(super) struct Wide {
hi: u128,
lo: u128,
}
impl Wide {
const ZERO: Wide = Wide { hi: 0, lo: 0 };
fn from_u128(v: u128) -> Wide {
Wide { hi: 0, lo: v }
}
fn shl(self, n: u32) -> Wide {
if n == 0 {
return self;
}
if n >= 256 {
return Wide::ZERO;
}
if n >= 128 {
Wide {
hi: self.lo << (n - 128),
lo: 0,
}
} else {
Wide {
hi: (self.hi << n) | (self.lo >> (128 - n)),
lo: self.lo << n,
}
}
}
fn shr_sticky(self, n: u32) -> Wide {
if n == 0 {
return self;
}
if n >= 256 {
return Wide::from_u128(u128::from(self != Wide::ZERO));
}
if n >= 128 {
let lost = self.lo != 0 || (n > 128 && self.hi & ((1u128 << (n - 128)) - 1) != 0);
Wide {
hi: 0,
lo: (self.hi >> (n - 128)) | u128::from(lost),
}
} else {
let lost = self.lo & ((1u128 << n) - 1) != 0;
Wide {
hi: self.hi >> n,
lo: (self.hi << (128 - n)) | (self.lo >> n) | u128::from(lost),
}
}
}
fn add(self, other: Wide) -> Wide {
let (lo, carry) = self.lo.overflowing_add(other.lo);
Wide {
hi: self
.hi
.wrapping_add(other.hi)
.wrapping_add(u128::from(carry)),
lo,
}
}
fn sub(self, other: Wide) -> Wide {
let (lo, borrow) = self.lo.overflowing_sub(other.lo);
Wide {
hi: self
.hi
.wrapping_sub(other.hi)
.wrapping_sub(u128::from(borrow)),
lo,
}
}
fn narrow(self) -> (u128, u32) {
if self.hi == 0 {
return (self.lo, 0);
}
let shift = 128 - self.hi.leading_zeros();
(self.shr_sticky(shift).lo, shift)
}
}
pub(super) fn nan_result(ops: &[Parts], env: Env) -> Option<(Outcome, Flags)> {
let mut any = false;
let mut signaling = false;
for p in ops {
if p.class == Class::Nan {
any = true;
signaling |= p.snan;
}
}
if !any {
return None;
}
let flags = if signaling {
Flags::INVALID
} else {
Flags::NONE
};
let pick = |p: &Parts| Outcome::Nan {
sign: p.sign,
payload: p.frac,
};
let first = |want_snan: bool| {
ops.iter()
.find(|p| p.class == Class::Nan && p.snan == want_snan)
};
let out = match env.nan.propagate {
Propagate::Default => Outcome::DefaultNan,
Propagate::FirstNan => pick(ops.iter().find(|p| p.class == Class::Nan)?),
Propagate::SignalingFirst => pick(first(true).or_else(|| first(false))?),
Propagate::LargerSignificand => {
let quiet = first(false).is_some();
pick(
ops.iter()
.rev()
.filter(|p| p.class == Class::Nan && p.snan != quiet)
.max_by_key(|p| p.frac)?,
)
}
};
Some((out, flags))
}
pub(super) fn round_exact(
sign: bool,
exp: i32,
sig: u128,
spec: Spec,
env: Env,
) -> (Outcome, Flags) {
if sig == 0 {
return (Outcome::Num(sign, Rounded::Zero), Flags::NONE);
}
let p = spec.precision;
let target = p + EXTRA - 1;
let msb = 127 - sig.leading_zeros();
let (mut sig, mut exp) = if msb > target {
let n = msb - target;
(shr_sticky(sig, n), exp + n as i32)
} else {
let n = target - msb;
(sig << n, exp - n as i32)
};
let lead = exp + target as i32;
let tiny_before = lead < spec.emin();
if exp + (EXTRA as i32) < spec.min_ulp {
let n = (spec.min_ulp - (exp + EXTRA as i32)) as u32;
sig = shr_sticky(sig, n);
exp += n as i32;
}
let rem = sig & ((1 << EXTRA) - 1);
let half = 1u128 << (EXTRA - 1);
let mut keep = sig >> EXTRA;
let inexact = rem != 0;
let up = match env.round {
Round::TiesEven => rem > half || (rem == half && keep & 1 != 0),
Round::TiesAway => rem >= half,
Round::TowardZero => false,
Round::TowardNegative => inexact && sign,
Round::TowardPositive => inexact && !sign,
};
let mut ulp = exp + EXTRA as i32;
if up {
keep += 1;
if keep >> p != 0 {
keep >>= 1;
ulp += 1;
}
}
if keep == 0 {
return (
Outcome::Num(sign, Rounded::Zero),
Flags::INEXACT | Flags::UNDERFLOW,
);
}
let lead_final = ulp + (127 - keep.leading_zeros()) as i32;
if lead_final > spec.emax {
let to_inf = match env.round {
Round::TiesEven | Round::TiesAway => true,
Round::TowardZero => false,
Round::TowardNegative => sign,
Round::TowardPositive => !sign,
};
let out = if to_inf {
Rounded::Inf
} else {
Rounded::Finite {
exp: spec.emax - (p as i32 - 1),
frac: u64::MAX >> (64 - p),
}
};
return (Outcome::Num(sign, out), Flags::OVERFLOW | Flags::INEXACT);
}
let tiny = match env.tininess {
Tininess::AfterRounding => lead_final < spec.emin(),
Tininess::BeforeRounding => tiny_before,
};
if tiny && env.flush_outputs {
return (
Outcome::Num(sign, Rounded::Zero),
Flags::INEXACT | Flags::UNDERFLOW,
);
}
let mut flags = if inexact { Flags::INEXACT } else { Flags::NONE };
if tiny && inexact {
flags |= Flags::UNDERFLOW;
}
(
Outcome::Num(
sign,
Rounded::Finite {
exp: ulp,
frac: keep as u64,
},
),
flags,
)
}
fn cancelled_sign(env: Env) -> bool {
env.round == Round::TowardNegative
}
pub(super) fn add(a: Parts, b: Parts, negate: bool, spec: Spec, env: Env) -> (Outcome, Flags) {
let mut b = b;
if negate {
b.sign = !b.sign;
}
if let Some(r) = nan_result(&[a, b], env) {
return r;
}
match (a.class, b.class) {
(Class::Inf, Class::Inf) => {
if a.sign == b.sign {
(Outcome::Num(a.sign, Rounded::Inf), Flags::NONE)
} else {
(Outcome::DefaultNan, Flags::INVALID)
}
}
(Class::Inf, _) => (Outcome::Num(a.sign, Rounded::Inf), Flags::NONE),
(_, Class::Inf) => (Outcome::Num(b.sign, Rounded::Inf), Flags::NONE),
(Class::Zero, Class::Zero) => {
let sign = if a.sign == b.sign {
a.sign
} else {
cancelled_sign(env)
};
(Outcome::Num(sign, Rounded::Zero), Flags::NONE)
}
(Class::Zero, _) => round_exact(b.sign, b.exp, u128::from(b.frac), spec, env),
(_, Class::Zero) => round_exact(a.sign, a.exp, u128::from(a.frac), spec, env),
_ => {
let (a, b) = (a.normalized(), b.normalized());
let (hi, lo) = if a.exp >= b.exp { (a, b) } else { (b, a) };
const SLACK: u32 = 62;
let sig_hi = u128::from(hi.frac) << SLACK;
let sig_lo = shr_sticky(u128::from(lo.frac) << SLACK, (hi.exp - lo.exp) as u32);
let (sign, sig) = if hi.sign == lo.sign {
(hi.sign, sig_hi + sig_lo)
} else if sig_hi >= sig_lo {
(hi.sign, sig_hi - sig_lo)
} else {
(lo.sign, sig_lo - sig_hi)
};
if sig == 0 {
return (
Outcome::Num(cancelled_sign(env), Rounded::Zero),
Flags::NONE,
);
}
round_exact(sign, hi.exp - SLACK as i32, sig, spec, env)
}
}
}
pub(super) fn mul(a: Parts, b: Parts, spec: Spec, env: Env) -> (Outcome, Flags) {
if let Some(r) = nan_result(&[a, b], env) {
return r;
}
let sign = a.sign ^ b.sign;
match (a.class, b.class) {
(Class::Inf, Class::Zero) | (Class::Zero, Class::Inf) => {
(Outcome::DefaultNan, Flags::INVALID)
}
(Class::Inf, _) | (_, Class::Inf) => (Outcome::Num(sign, Rounded::Inf), Flags::NONE),
(Class::Zero, _) | (_, Class::Zero) => (Outcome::Num(sign, Rounded::Zero), Flags::NONE),
_ => {
let sig = u128::from(a.frac) * u128::from(b.frac);
round_exact(sign, a.exp + b.exp, sig, spec, env)
}
}
}
pub(super) fn div(a: Parts, b: Parts, spec: Spec, env: Env) -> (Outcome, Flags) {
if let Some(r) = nan_result(&[a, b], env) {
return r;
}
let sign = a.sign ^ b.sign;
match (a.class, b.class) {
(Class::Inf, Class::Inf) | (Class::Zero, Class::Zero) => {
(Outcome::DefaultNan, Flags::INVALID)
}
(Class::Inf, _) => (Outcome::Num(sign, Rounded::Inf), Flags::NONE),
(_, Class::Inf) => (Outcome::Num(sign, Rounded::Zero), Flags::NONE),
(_, Class::Zero) => (Outcome::Num(sign, Rounded::Inf), Flags::DIV_BY_ZERO),
(Class::Zero, _) => (Outcome::Num(sign, Rounded::Zero), Flags::NONE),
_ => {
let (a, b) = (a.normalized(), b.normalized());
let den = u128::from(b.frac);
let n1 = u128::from(a.frac) << 63;
let (q1, r1) = (n1 / den, n1 % den);
let n2 = r1 << 8;
let (q2, r2) = (n2 / den, n2 % den);
let mut q = (q1 << 8) | q2;
if r2 != 0 {
q |= 1;
}
round_exact(sign, a.exp - b.exp - 71, q, spec, env)
}
}
}
fn isqrt(radicand: u128, extra: u32) -> (u128, bool) {
if radicand == 0 {
return (0, false);
}
let mut rem: u128 = 0;
let mut root: u128 = 0;
let mut shift = (127 - radicand.leading_zeros()) & !1;
loop {
rem = (rem << 2) | ((radicand >> shift) & 3);
root <<= 1;
let trial = (root << 1) | 1;
if rem >= trial {
rem -= trial;
root |= 1;
}
if shift == 0 {
break;
}
shift -= 2;
}
for _ in 0..extra {
rem <<= 2;
root <<= 1;
let trial = (root << 1) | 1;
if rem >= trial {
rem -= trial;
root |= 1;
}
}
(root, rem != 0)
}
pub(super) fn sqrt(a: Parts, spec: Spec, env: Env) -> (Outcome, Flags) {
if let Some(r) = nan_result(&[a], env) {
return r;
}
match a.class {
Class::Zero => (Outcome::Num(a.sign, Rounded::Zero), Flags::NONE),
_ if a.sign => (Outcome::DefaultNan, Flags::INVALID),
Class::Inf => (Outcome::Num(false, Rounded::Inf), Flags::NONE),
_ => {
let a = a.normalized();
let mut frac = u128::from(a.frac);
let mut exp = a.exp;
if exp & 1 != 0 {
frac <<= 1;
exp -= 1;
}
const EXTRA_PAIRS: u32 = 40;
let (root, rest) = isqrt(frac, EXTRA_PAIRS);
let sig = if rest { root | 1 } else { root };
round_exact(false, exp / 2 - EXTRA_PAIRS as i32, sig, spec, env)
}
}
}
pub(super) fn fma(a: Parts, b: Parts, c: Parts, spec: Spec, env: Env) -> (Outcome, Flags) {
if matches!(
(a.class, b.class),
(Class::Inf, Class::Zero) | (Class::Zero, Class::Inf)
) {
return (Outcome::DefaultNan, Flags::INVALID);
}
if let Some(r) = nan_result(&[a, b, c], env) {
return r;
}
let psign = a.sign ^ b.sign;
if a.class == Class::Inf || b.class == Class::Inf {
if c.class == Class::Inf && c.sign != psign {
return (Outcome::DefaultNan, Flags::INVALID);
}
return (Outcome::Num(psign, Rounded::Inf), Flags::NONE);
}
if c.class == Class::Inf {
return (Outcome::Num(c.sign, Rounded::Inf), Flags::NONE);
}
if a.class == Class::Zero || b.class == Class::Zero {
if c.class == Class::Zero {
let sign = if psign == c.sign {
psign
} else {
cancelled_sign(env)
};
return (Outcome::Num(sign, Rounded::Zero), Flags::NONE);
}
return round_exact(c.sign, c.exp, u128::from(c.frac), spec, env);
}
if c.class == Class::Zero {
return mul(a, b, spec, env);
}
let (a, b, c) = (a.normalized(), b.normalized(), c.normalized());
let prod = u128::from(a.frac) * u128::from(b.frac);
let pexp = a.exp + b.exp;
let base = core::cmp::max(pexp + 128, c.exp + 64) - 250;
let place = |v: u128, exp: i32| {
let shift = exp - base;
if shift >= 0 {
Wide::from_u128(v).shl(shift as u32)
} else {
Wide::from_u128(v).shr_sticky((-shift) as u32)
}
};
let psig = place(prod, pexp);
let csig = place(u128::from(c.frac), c.exp);
let (sign, sum) = if psign == c.sign {
(psign, psig.add(csig))
} else if psig >= csig {
(psign, psig.sub(csig))
} else {
(c.sign, csig.sub(psig))
};
if sum == Wide::ZERO {
return (
Outcome::Num(cancelled_sign(env), Rounded::Zero),
Flags::NONE,
);
}
let (sig, shift) = sum.narrow();
round_exact(sign, base + shift as i32, sig, spec, env)
}
pub(super) fn compare(a: Parts, b: Parts) -> core::cmp::Ordering {
use core::cmp::Ordering;
let magnitude = |p: &Parts| match p.class {
Class::Zero => (i32::MIN, 0u64),
Class::Inf => (i32::MAX, u64::MAX),
_ => {
let n = p.normalized();
(n.exp, n.frac)
}
};
if a.class == Class::Zero && b.class == Class::Zero {
return Ordering::Equal;
}
match (a.sign, b.sign) {
(false, true) => Ordering::Greater,
(true, false) => Ordering::Less,
(false, false) => magnitude(&a).cmp(&magnitude(&b)),
(true, true) => magnitude(&b).cmp(&magnitude(&a)),
}
}
pub(super) fn min_max(
a: Parts,
b: Parts,
want_min: bool,
env: Env,
) -> (Result<usize, Outcome>, Flags) {
use core::cmp::Ordering;
let nan_a = a.class == Class::Nan;
let nan_b = b.class == Class::Nan;
let flags = if a.snan || b.snan {
Flags::INVALID
} else {
Flags::NONE
};
match env.min_max {
MinMax::NonNan => match (nan_a, nan_b) {
(true, true) => return (Err(Outcome::DefaultNan), flags),
(true, false) => return (Ok(1), flags),
(false, true) => return (Ok(0), flags),
(false, false) => {}
},
MinMax::SecondOperand => {
if nan_a || nan_b {
return (Ok(1), flags);
}
let ord = compare(a, b);
let first_wins = if want_min {
ord == Ordering::Less
} else {
ord == Ordering::Greater
};
return (Ok(usize::from(!first_wins)), flags);
}
MinMax::PropagateNan => {
if let Some((out, f)) = nan_result(&[a, b], env) {
return (Err(out), f);
}
}
}
let ord = if a.class == Class::Zero && b.class == Class::Zero {
match (a.sign, b.sign) {
(true, false) => Ordering::Less,
(false, true) => Ordering::Greater,
_ => Ordering::Equal,
}
} else {
compare(a, b)
};
let take_a = if want_min {
ord != Ordering::Greater
} else {
ord != Ordering::Less
};
(Ok(usize::from(!take_a)), flags)
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(super) enum IntValue {
Nan,
Inf(bool),
Value {
sign: bool,
magnitude: u128,
inexact: bool,
},
}
pub(super) fn to_integer(p: Parts, env: Env) -> IntValue {
match p.class {
Class::Nan => IntValue::Nan,
Class::Inf => IntValue::Inf(p.sign),
Class::Zero => IntValue::Value {
sign: p.sign,
magnitude: 0,
inexact: false,
},
Class::Finite => {
let shift = p.exp + EXTRA as i32;
let scaled: u128 = if shift >= 0 {
if shift > 63 {
return IntValue::Value {
sign: p.sign,
magnitude: u128::MAX,
inexact: false,
};
}
u128::from(p.frac) << shift
} else {
shr_sticky(u128::from(p.frac), (-shift) as u32)
};
let rem = scaled & ((1 << EXTRA) - 1);
let half = 1u128 << (EXTRA - 1);
let mut magnitude = scaled >> EXTRA;
let inexact = rem != 0;
let up = match env.round {
Round::TiesEven => rem > half || (rem == half && magnitude & 1 != 0),
Round::TiesAway => rem >= half,
Round::TowardZero => false,
Round::TowardNegative => inexact && p.sign,
Round::TowardPositive => inexact && !p.sign,
};
if up {
magnitude += 1;
}
IntValue::Value {
sign: p.sign,
magnitude,
inexact,
}
}
}
}