1 use crate::{Category, ExpInt, IEK_INF, IEK_NAN, IEK_ZERO};
2 use crate::{Float, FloatConvert, ParseError, Round, Status, StatusAnd};
4 use core::cmp::{self, Ordering};
5 use core::convert::TryFrom;
6 use core::fmt::{self, Write};
7 use core::marker::PhantomData;
10 use smallvec::{smallvec, SmallVec};
13 pub struct IeeeFloat<S> {
14 /// Absolute significand value (including the integer bit).
17 /// The signed unbiased exponent of the value.
20 /// What kind of floating point number this is.
23 /// Sign bit of the number.
26 marker: PhantomData<S>,
29 /// Fundamental unit of big integer arithmetic, but also
30 /// large to store the largest significands by itself.
32 const LIMB_BITS: usize = 128;
33 fn limbs_for_bits(bits: usize) -> usize {
34 (bits + LIMB_BITS - 1) / LIMB_BITS
37 /// Enum that represents what fraction of the LSB truncated bits of an fp number
40 /// This essentially combines the roles of guard and sticky bits.
42 #[derive(Copy, Clone, PartialEq, Eq, Debug)]
44 // Example of truncated bits:
45 ExactlyZero, // 000000
46 LessThanHalf, // 0xxxxx x's not all zero
47 ExactlyHalf, // 100000
48 MoreThanHalf, // 1xxxxx x's not all zero
51 /// Represents floating point arithmetic semantics.
52 pub trait Semantics: Sized {
53 /// Total number of bits in the in-memory format.
56 /// Number of bits in the significand. This includes the integer bit.
57 const PRECISION: usize;
59 /// The largest E such that 2<sup>E</sup> is representable; this matches the
60 /// definition of IEEE 754.
61 const MAX_EXP: ExpInt;
63 /// The smallest E such that 2<sup>E</sup> is a normalized number; this
64 /// matches the definition of IEEE 754.
65 const MIN_EXP: ExpInt = -Self::MAX_EXP + 1;
67 /// The significand bit that marks NaN as quiet.
68 const QNAN_BIT: usize = Self::PRECISION - 2;
70 /// The significand bitpattern to mark a NaN as quiet.
71 /// NOTE: for X87DoubleExtended we need to set two bits instead of 2.
72 const QNAN_SIGNIFICAND: Limb = 1 << Self::QNAN_BIT;
74 fn from_bits(bits: u128) -> IeeeFloat<Self> {
75 assert!(Self::BITS > Self::PRECISION);
77 let sign = bits & (1 << (Self::BITS - 1));
78 let exponent = (bits & !sign) >> (Self::PRECISION - 1);
79 let mut r = IeeeFloat {
80 sig: [bits & ((1 << (Self::PRECISION - 1)) - 1)],
81 // Convert the exponent from its bias representation to a signed integer.
82 exp: (exponent as ExpInt) - Self::MAX_EXP,
83 category: Category::Zero,
88 if r.exp == Self::MIN_EXP - 1 && r.sig == [0] {
89 // Exponent, significand meaningless.
90 r.category = Category::Zero;
91 } else if r.exp == Self::MAX_EXP + 1 && r.sig == [0] {
92 // Exponent, significand meaningless.
93 r.category = Category::Infinity;
94 } else if r.exp == Self::MAX_EXP + 1 && r.sig != [0] {
95 // Sign, exponent, significand meaningless.
96 r.category = Category::NaN;
98 r.category = Category::Normal;
99 if r.exp == Self::MIN_EXP - 1 {
101 r.exp = Self::MIN_EXP;
104 sig::set_bit(&mut r.sig, Self::PRECISION - 1);
111 fn to_bits(x: IeeeFloat<Self>) -> u128 {
112 assert!(Self::BITS > Self::PRECISION);
114 // Split integer bit from significand.
115 let integer_bit = sig::get_bit(&x.sig, Self::PRECISION - 1);
116 let mut significand = x.sig[0] & ((1 << (Self::PRECISION - 1)) - 1);
117 let exponent = match x.category {
118 Category::Normal => {
119 if x.exp == Self::MIN_EXP && !integer_bit {
127 // FIXME(eddyb) Maybe we should guarantee an invariant instead?
131 Category::Infinity => {
132 // FIXME(eddyb) Maybe we should guarantee an invariant instead?
136 Category::NaN => Self::MAX_EXP + 1,
139 // Convert the exponent from a signed integer to its bias representation.
140 let exponent = (exponent + Self::MAX_EXP) as u128;
142 ((x.sign as u128) << (Self::BITS - 1)) | (exponent << (Self::PRECISION - 1)) | significand
146 impl<S> Copy for IeeeFloat<S> {}
147 impl<S> Clone for IeeeFloat<S> {
148 fn clone(&self) -> Self {
153 macro_rules! ieee_semantics {
154 ($($name:ident = $sem:ident($bits:tt : $exp_bits:tt)),*) => {
156 $(pub type $name = IeeeFloat<$sem>;)*
157 $(impl Semantics for $sem {
158 const BITS: usize = $bits;
159 const PRECISION: usize = ($bits - 1 - $exp_bits) + 1;
160 const MAX_EXP: ExpInt = (1 << ($exp_bits - 1)) - 1;
167 Single = SingleS(32:8),
168 Double = DoubleS(64:11),
172 pub struct X87DoubleExtendedS;
173 pub type X87DoubleExtended = IeeeFloat<X87DoubleExtendedS>;
174 impl Semantics for X87DoubleExtendedS {
175 const BITS: usize = 80;
176 const PRECISION: usize = 64;
177 const MAX_EXP: ExpInt = (1 << (15 - 1)) - 1;
179 /// For x87 extended precision, we want to make a NaN, not a
180 /// pseudo-NaN. Maybe we should expose the ability to make
182 const QNAN_SIGNIFICAND: Limb = 0b11 << Self::QNAN_BIT;
184 /// Integer bit is explicit in this format. Intel hardware (387 and later)
185 /// does not support these bit patterns:
186 /// exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
187 /// exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
188 /// exponent = 0, integer bit 1 ("pseudodenormal")
189 /// exponent != 0 nor all 1's, integer bit 0 ("unnormal")
190 /// At the moment, the first two are treated as NaNs, the second two as Normal.
191 fn from_bits(bits: u128) -> IeeeFloat<Self> {
192 let sign = bits & (1 << (Self::BITS - 1));
193 let exponent = (bits & !sign) >> Self::PRECISION;
194 let mut r = IeeeFloat {
195 sig: [bits & ((1 << (Self::PRECISION - 1)) - 1)],
196 // Convert the exponent from its bias representation to a signed integer.
197 exp: (exponent as ExpInt) - Self::MAX_EXP,
198 category: Category::Zero,
203 if r.exp == Self::MIN_EXP - 1 && r.sig == [0] {
204 // Exponent, significand meaningless.
205 r.category = Category::Zero;
206 } else if r.exp == Self::MAX_EXP + 1 && r.sig == [1 << (Self::PRECISION - 1)] {
207 // Exponent, significand meaningless.
208 r.category = Category::Infinity;
209 } else if r.exp == Self::MAX_EXP + 1 && r.sig != [1 << (Self::PRECISION - 1)] {
210 // Sign, exponent, significand meaningless.
211 r.category = Category::NaN;
213 r.category = Category::Normal;
214 if r.exp == Self::MIN_EXP - 1 {
216 r.exp = Self::MIN_EXP;
223 fn to_bits(x: IeeeFloat<Self>) -> u128 {
224 // Get integer bit from significand.
225 let integer_bit = sig::get_bit(&x.sig, Self::PRECISION - 1);
226 let mut significand = x.sig[0] & ((1 << Self::PRECISION) - 1);
227 let exponent = match x.category {
228 Category::Normal => {
229 if x.exp == Self::MIN_EXP && !integer_bit {
237 // FIXME(eddyb) Maybe we should guarantee an invariant instead?
241 Category::Infinity => {
242 // FIXME(eddyb) Maybe we should guarantee an invariant instead?
243 significand = 1 << (Self::PRECISION - 1);
246 Category::NaN => Self::MAX_EXP + 1,
249 // Convert the exponent from a signed integer to its bias representation.
250 let exponent = (exponent + Self::MAX_EXP) as u128;
252 ((x.sign as u128) << (Self::BITS - 1)) | (exponent << Self::PRECISION) | significand
256 float_common_impls!(IeeeFloat<S>);
258 impl<S: Semantics> PartialEq for IeeeFloat<S> {
259 fn eq(&self, rhs: &Self) -> bool {
260 self.partial_cmp(rhs) == Some(Ordering::Equal)
264 impl<S: Semantics> PartialOrd for IeeeFloat<S> {
265 fn partial_cmp(&self, rhs: &Self) -> Option<Ordering> {
266 match (self.category, rhs.category) {
267 (Category::NaN, _) | (_, Category::NaN) => None,
269 (Category::Infinity, Category::Infinity) => Some((!self.sign).cmp(&(!rhs.sign))),
271 (Category::Zero, Category::Zero) => Some(Ordering::Equal),
273 (Category::Infinity, _) | (Category::Normal, Category::Zero) => {
274 Some((!self.sign).cmp(&self.sign))
277 (_, Category::Infinity) | (Category::Zero, Category::Normal) => {
278 Some(rhs.sign.cmp(&(!rhs.sign)))
281 (Category::Normal, Category::Normal) => {
282 // Two normal numbers. Do they have the same sign?
283 Some((!self.sign).cmp(&(!rhs.sign)).then_with(|| {
284 // Compare absolute values; invert result if negative.
285 let result = self.cmp_abs_normal(*rhs);
287 if self.sign { result.reverse() } else { result }
294 impl<S> Neg for IeeeFloat<S> {
296 fn neg(mut self) -> Self {
297 self.sign = !self.sign;
302 /// Prints this value as a decimal string.
304 /// \param precision The maximum number of digits of
305 /// precision to output. If there are fewer digits available,
306 /// zero padding will not be used unless the value is
307 /// integral and small enough to be expressed in
308 /// precision digits. 0 means to use the natural
309 /// precision of the number.
310 /// \param width The maximum number of zeros to
311 /// consider inserting before falling back to scientific
312 /// notation. 0 means to always use scientific notation.
314 /// \param alternate Indicate whether to remove the trailing zero in
315 /// fraction part or not. Also setting this parameter to true forces
316 /// producing of output more similar to default printf behavior.
317 /// Specifically the lower e is used as exponent delimiter and exponent
318 /// always contains no less than two digits.
320 /// Number precision width Result
321 /// ------ --------- ----- ------
322 /// 1.01E+4 5 2 10100
323 /// 1.01E+4 4 2 1.01E+4
324 /// 1.01E+4 5 1 1.01E+4
325 /// 1.01E-2 5 2 0.0101
326 /// 1.01E-2 4 2 0.0101
327 /// 1.01E-2 4 1 1.01E-2
328 impl<S: Semantics> fmt::Display for IeeeFloat<S> {
329 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
330 let width = f.width().unwrap_or(3);
331 let alternate = f.alternate();
333 match self.category {
334 Category::Infinity => {
336 return f.write_str("-Inf");
338 return f.write_str("+Inf");
342 Category::NaN => return f.write_str("NaN"),
352 if let Some(n) = f.precision() {
357 f.write_str("e+00")?;
359 f.write_str("0.0E+0")?;
367 Category::Normal => {}
374 // We use enough digits so the number can be round-tripped back to an
375 // APFloat. The formula comes from "How to Print Floating-Point Numbers
376 // Accurately" by Steele and White.
377 // FIXME: Using a formula based purely on the precision is conservative;
378 // we can print fewer digits depending on the actual value being printed.
380 // precision = 2 + floor(S::PRECISION / lg_2(10))
381 let precision = f.precision().unwrap_or(2 + S::PRECISION * 59 / 196);
383 // Decompose the number into an APInt and an exponent.
384 let mut exp = self.exp - (S::PRECISION as ExpInt - 1);
385 let mut sig = vec![self.sig[0]];
387 // Ignore trailing binary zeros.
388 let trailing_zeros = sig[0].trailing_zeros();
389 let _: Loss = sig::shift_right(&mut sig, &mut exp, trailing_zeros as usize);
391 // Change the exponent from 2^e to 10^e.
396 let shift = exp as usize;
397 sig.resize(limbs_for_bits(S::PRECISION + shift), 0);
398 sig::shift_left(&mut sig, &mut exp, shift);
401 let mut texp = -exp as usize;
403 // We transform this using the identity:
404 // (N)(2^-e) == (N)(5^e)(10^-e)
406 // Multiply significand by 5^e.
407 // N * 5^0101 == N * 5^(1*1) * 5^(0*2) * 5^(1*4) * 5^(0*8)
408 let mut sig_scratch = vec![];
410 let mut p5_scratch = vec![];
415 p5_scratch.resize(p5.len() * 2, 0);
417 sig::mul(&mut p5_scratch, &mut 0, &p5, &p5, p5.len() * 2 * LIMB_BITS);
418 while p5_scratch.last() == Some(&0) {
421 mem::swap(&mut p5, &mut p5_scratch);
424 sig_scratch.resize(sig.len() + p5.len(), 0);
425 let _: Loss = sig::mul(
430 (sig.len() + p5.len()) * LIMB_BITS,
432 while sig_scratch.last() == Some(&0) {
435 mem::swap(&mut sig, &mut sig_scratch);
442 let mut buffer = vec![];
444 // Ignore digits from the significand until it is no more
445 // precise than is required for the desired precision.
446 // 196/59 is a very slight overestimate of lg_2(10).
447 let required = (precision * 196 + 58) / 59;
448 let mut discard_digits = sig::omsb(&sig).saturating_sub(required) * 59 / 196;
449 let mut in_trail = true;
450 while !sig.is_empty() {
451 // Perform short division by 10 to extract the rightmost digit.
456 // Use 64-bit division and remainder, with 32-bit chunks from sig.
457 sig::each_chunk(&mut sig, 32, |chunk| {
458 let chunk = chunk as u32;
459 let combined = ((rem as u64) << 32) | (chunk as u64);
460 rem = (combined % 10) as u8;
461 (combined / 10) as u32 as Limb
464 // Reduce the sigificand to avoid wasting time dividing 0's.
465 while sig.last() == Some(&0) {
471 // Ignore digits we don't need.
472 if discard_digits > 0 {
478 // Drop trailing zeros.
479 if in_trail && digit == 0 {
483 buffer.push(b'0' + digit);
487 assert!(!buffer.is_empty(), "no characters in buffer!");
489 // Drop down to precision.
490 // FIXME: don't do more precise calculations above than are required.
491 if buffer.len() > precision {
492 // The most significant figures are the last ones in the buffer.
493 let mut first_sig = buffer.len() - precision;
496 // FIXME: this probably shouldn't use 'round half up'.
498 // Rounding down is just a truncation, except we also want to drop
499 // trailing zeros from the new result.
500 if buffer[first_sig - 1] < b'5' {
501 while first_sig < buffer.len() && buffer[first_sig] == b'0' {
505 // Rounding up requires a decimal add-with-carry. If we continue
506 // the carry, the newly-introduced zeros will just be truncated.
507 for x in &mut buffer[first_sig..] {
517 exp += first_sig as ExpInt;
518 buffer.drain(..first_sig);
520 // If we carried through, we have exactly one digit of precision.
521 if buffer.is_empty() {
526 let digits = buffer.len();
528 // Check whether we should use scientific notation.
529 let scientific = if width == 0 {
534 // But we shouldn't make the number look more precise than it is.
535 exp as usize > width || digits + exp as usize > precision
537 // Power of the most significant digit.
538 let msd = exp + (digits - 1) as ExpInt;
545 -msd as usize > width
549 // Scientific formatting is pretty straightforward.
551 exp += digits as ExpInt - 1;
553 f.write_char(buffer[digits - 1] as char)?;
555 let truncate_zero = !alternate;
556 if digits == 1 && truncate_zero {
559 for &d in buffer[..digits - 1].iter().rev() {
560 f.write_char(d as char)?;
563 // Fill with zeros up to precision.
564 if !truncate_zero && precision > digits - 1 {
565 for _ in 0..=precision - digits {
569 // For alternate we use lower 'e'.
570 f.write_char(if alternate { 'e' } else { 'E' })?;
572 // Exponent always at least two digits if we do not truncate zeros.
574 write!(f, "{:+}", exp)?;
576 write!(f, "{:+03}", exp)?;
582 // Non-scientific, positive exponents.
584 for &d in buffer.iter().rev() {
585 f.write_char(d as char)?;
593 // Non-scientific, negative exponents.
594 let unit_place = -exp as usize;
595 if unit_place < digits {
596 for &d in buffer[unit_place..].iter().rev() {
597 f.write_char(d as char)?;
600 for &d in buffer[..unit_place].iter().rev() {
601 f.write_char(d as char)?;
605 for _ in digits..unit_place {
608 for &d in buffer.iter().rev() {
609 f.write_char(d as char)?;
617 impl<S: Semantics> fmt::Debug for IeeeFloat<S> {
618 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
621 "{}({:?} | {}{:?} * 2^{})",
624 if self.sign { "-" } else { "+" },
631 impl<S: Semantics> Float for IeeeFloat<S> {
632 const BITS: usize = S::BITS;
633 const PRECISION: usize = S::PRECISION;
634 const MAX_EXP: ExpInt = S::MAX_EXP;
635 const MIN_EXP: ExpInt = S::MIN_EXP;
637 const ZERO: Self = IeeeFloat {
640 category: Category::Zero,
645 const INFINITY: Self = IeeeFloat {
648 category: Category::Infinity,
653 // FIXME(eddyb) remove when qnan becomes const fn.
654 const NAN: Self = IeeeFloat {
655 sig: [S::QNAN_SIGNIFICAND],
657 category: Category::NaN,
662 fn qnan(payload: Option<u128>) -> Self {
664 sig: [S::QNAN_SIGNIFICAND
665 | payload.map_or(0, |payload| {
666 // Zero out the excess bits of the significand.
667 payload & ((1 << S::QNAN_BIT) - 1)
670 category: Category::NaN,
676 fn snan(payload: Option<u128>) -> Self {
677 let mut snan = Self::qnan(payload);
679 // We always have to clear the QNaN bit to make it an SNaN.
680 sig::clear_bit(&mut snan.sig, S::QNAN_BIT);
682 // If there are no bits set in the payload, we have to set
683 // *something* to make it a NaN instead of an infinity;
684 // conventionally, this is the next bit down from the QNaN bit.
685 if snan.sig[0] & !S::QNAN_SIGNIFICAND == 0 {
686 sig::set_bit(&mut snan.sig, S::QNAN_BIT - 1);
692 fn largest() -> Self {
693 // We want (in interchange format):
695 // significand = 1..1
697 sig: [(1 << S::PRECISION) - 1],
699 category: Category::Normal,
705 // We want (in interchange format):
707 // significand = 0..01
708 const SMALLEST: Self = IeeeFloat {
711 category: Category::Normal,
716 fn smallest_normalized() -> Self {
717 // We want (in interchange format):
719 // significand = 10..0
721 sig: [1 << (S::PRECISION - 1)],
723 category: Category::Normal,
729 fn add_r(mut self, rhs: Self, round: Round) -> StatusAnd<Self> {
730 let status = match (self.category, rhs.category) {
731 (Category::Infinity, Category::Infinity) => {
732 // Differently signed infinities can only be validly
734 if self.sign != rhs.sign {
742 // Sign may depend on rounding mode; handled below.
743 (_, Category::Zero) | (Category::NaN, _) | (Category::Infinity, Category::Normal) => {
747 (Category::Zero, _) | (_, Category::NaN) | (_, Category::Infinity) => {
752 // This return code means it was not a simple case.
753 (Category::Normal, Category::Normal) => {
754 let loss = sig::add_or_sub(
763 self = unpack!(status=, self.normalize(round, loss));
765 // Can only be zero if we lost no fraction.
766 assert!(self.category != Category::Zero || loss == Loss::ExactlyZero);
772 // If two numbers add (exactly) to zero, IEEE 754 decrees it is a
773 // positive zero unless rounding to minus infinity, except that
774 // adding two like-signed zeroes gives that zero.
775 if self.category == Category::Zero
776 && (rhs.category != Category::Zero || self.sign != rhs.sign)
778 self.sign = round == Round::TowardNegative;
784 fn mul_r(mut self, rhs: Self, round: Round) -> StatusAnd<Self> {
785 self.sign ^= rhs.sign;
787 match (self.category, rhs.category) {
788 (Category::NaN, _) => {
793 (_, Category::NaN) => {
795 self.category = Category::NaN;
800 (Category::Zero, Category::Infinity) | (Category::Infinity, Category::Zero) => {
801 Status::INVALID_OP.and(Self::NAN)
804 (_, Category::Infinity) | (Category::Infinity, _) => {
805 self.category = Category::Infinity;
809 (Category::Zero, _) | (_, Category::Zero) => {
810 self.category = Category::Zero;
814 (Category::Normal, Category::Normal) => {
816 let mut wide_sig = [0; 2];
818 sig::mul(&mut wide_sig, &mut self.exp, &self.sig, &rhs.sig, S::PRECISION);
819 self.sig = [wide_sig[0]];
821 self = unpack!(status=, self.normalize(round, loss));
822 if loss != Loss::ExactlyZero {
823 status |= Status::INEXACT;
830 fn mul_add_r(mut self, multiplicand: Self, addend: Self, round: Round) -> StatusAnd<Self> {
831 // If and only if all arguments are normal do we need to do an
832 // extended-precision calculation.
833 if !self.is_finite_non_zero() || !multiplicand.is_finite_non_zero() || !addend.is_finite() {
835 self = unpack!(status=, self.mul_r(multiplicand, round));
837 // FS can only be Status::OK or Status::INVALID_OP. There is no more work
838 // to do in the latter case. The IEEE-754R standard says it is
839 // implementation-defined in this case whether, if ADDEND is a
840 // quiet NaN, we raise invalid op; this implementation does so.
842 // If we need to do the addition we can do so with normal
844 if status == Status::OK {
845 self = unpack!(status=, self.add_r(addend, round));
847 return status.and(self);
850 // Post-multiplication sign, before addition.
851 self.sign ^= multiplicand.sign;
853 // Allocate space for twice as many bits as the original significand, plus one
854 // extra bit for the addition to overflow into.
855 assert!(limbs_for_bits(S::PRECISION * 2 + 1) <= 2);
856 let mut wide_sig = sig::widening_mul(self.sig[0], multiplicand.sig[0]);
858 let mut loss = Loss::ExactlyZero;
859 let mut omsb = sig::omsb(&wide_sig);
860 self.exp += multiplicand.exp;
862 // Assume the operands involved in the multiplication are single-precision
863 // FP, and the two multiplicants are:
864 // lhs = a23 . a22 ... a0 * 2^e1
865 // rhs = b23 . b22 ... b0 * 2^e2
866 // the result of multiplication is:
867 // lhs = c48 c47 c46 . c45 ... c0 * 2^(e1+e2)
868 // Note that there are three significant bits at the left-hand side of the
869 // radix point: two for the multiplication, and an overflow bit for the
870 // addition (that will always be zero at this point). Move the radix point
871 // toward left by two bits, and adjust exponent accordingly.
874 if addend.is_non_zero() {
875 // Normalize our MSB to one below the top bit to allow for overflow.
876 let ext_precision = 2 * S::PRECISION + 1;
877 if omsb != ext_precision - 1 {
878 assert!(ext_precision > omsb);
879 sig::shift_left(&mut wide_sig, &mut self.exp, (ext_precision - 1) - omsb);
882 // The intermediate result of the multiplication has "2 * S::PRECISION"
883 // significant bit; adjust the addend to be consistent with mul result.
884 let mut ext_addend_sig = [addend.sig[0], 0];
886 // Extend the addend significand to ext_precision - 1. This guarantees
887 // that the high bit of the significand is zero (same as wide_sig),
888 // so the addition will overflow (if it does overflow at all) into the top bit.
889 sig::shift_left(&mut ext_addend_sig, &mut 0, ext_precision - 1 - S::PRECISION);
890 loss = sig::add_or_sub(
899 omsb = sig::omsb(&wide_sig);
902 // Convert the result having "2 * S::PRECISION" significant-bits back to the one
903 // having "S::PRECISION" significant-bits. First, move the radix point from
904 // position "2*S::PRECISION - 1" to "S::PRECISION - 1". The exponent need to be
905 // adjusted by "2*S::PRECISION - 1" - "S::PRECISION - 1" = "S::PRECISION".
906 self.exp -= S::PRECISION as ExpInt + 1;
908 // In case MSB resides at the left-hand side of radix point, shift the
909 // mantissa right by some amount to make sure the MSB reside right before
910 // the radix point (i.e., "MSB . rest-significant-bits").
911 if omsb > S::PRECISION {
912 let bits = omsb - S::PRECISION;
913 loss = sig::shift_right(&mut wide_sig, &mut self.exp, bits).combine(loss);
916 self.sig[0] = wide_sig[0];
919 self = unpack!(status=, self.normalize(round, loss));
920 if loss != Loss::ExactlyZero {
921 status |= Status::INEXACT;
924 // If two numbers add (exactly) to zero, IEEE 754 decrees it is a
925 // positive zero unless rounding to minus infinity, except that
926 // adding two like-signed zeroes gives that zero.
927 if self.category == Category::Zero
928 && !status.intersects(Status::UNDERFLOW)
929 && self.sign != addend.sign
931 self.sign = round == Round::TowardNegative;
937 fn div_r(mut self, rhs: Self, round: Round) -> StatusAnd<Self> {
938 self.sign ^= rhs.sign;
940 match (self.category, rhs.category) {
941 (Category::NaN, _) => {
946 (_, Category::NaN) => {
947 self.category = Category::NaN;
953 (Category::Infinity, Category::Infinity) | (Category::Zero, Category::Zero) => {
954 Status::INVALID_OP.and(Self::NAN)
957 (Category::Infinity, _) | (Category::Zero, _) => Status::OK.and(self),
959 (Category::Normal, Category::Infinity) => {
960 self.category = Category::Zero;
964 (Category::Normal, Category::Zero) => {
965 self.category = Category::Infinity;
966 Status::DIV_BY_ZERO.and(self)
969 (Category::Normal, Category::Normal) => {
971 let dividend = self.sig[0];
980 self = unpack!(status=, self.normalize(round, loss));
981 if loss != Loss::ExactlyZero {
982 status |= Status::INEXACT;
989 fn c_fmod(mut self, rhs: Self) -> StatusAnd<Self> {
990 match (self.category, rhs.category) {
992 | (Category::Zero, Category::Infinity)
993 | (Category::Zero, Category::Normal)
994 | (Category::Normal, Category::Infinity) => Status::OK.and(self),
996 (_, Category::NaN) => {
998 self.category = Category::NaN;
1000 Status::OK.and(self)
1003 (Category::Infinity, _) | (_, Category::Zero) => Status::INVALID_OP.and(Self::NAN),
1005 (Category::Normal, Category::Normal) => {
1006 while self.is_finite_non_zero()
1007 && rhs.is_finite_non_zero()
1008 && self.cmp_abs_normal(rhs) != Ordering::Less
1010 let mut v = rhs.scalbn(self.ilogb() - rhs.ilogb());
1011 if self.cmp_abs_normal(v) == Ordering::Less {
1017 self = unpack!(status=, self - v);
1018 assert_eq!(status, Status::OK);
1020 Status::OK.and(self)
1025 fn round_to_integral(self, round: Round) -> StatusAnd<Self> {
1026 // If the exponent is large enough, we know that this value is already
1027 // integral, and the arithmetic below would potentially cause it to saturate
1028 // to +/-Inf. Bail out early instead.
1029 if self.is_finite_non_zero() && self.exp + 1 >= S::PRECISION as ExpInt {
1030 return Status::OK.and(self);
1033 // The algorithm here is quite simple: we add 2^(p-1), where p is the
1034 // precision of our format, and then subtract it back off again. The choice
1035 // of rounding modes for the addition/subtraction determines the rounding mode
1036 // for our integral rounding as well.
1037 // NOTE: When the input value is negative, we do subtraction followed by
1038 // addition instead.
1039 assert!(S::PRECISION <= 128);
1041 let magic_const = unpack!(status=, Self::from_u128(1 << (S::PRECISION - 1)));
1042 let magic_const = magic_const.copy_sign(self);
1044 if status != Status::OK {
1045 return status.and(self);
1049 r = unpack!(status=, r.add_r(magic_const, round));
1050 if status != Status::OK && status != Status::INEXACT {
1051 return status.and(self);
1054 // Restore the input sign to handle 0.0/-0.0 cases correctly.
1055 r.sub_r(magic_const, round).map(|r| r.copy_sign(self))
1058 fn next_up(mut self) -> StatusAnd<Self> {
1059 // Compute nextUp(x), handling each float category separately.
1060 match self.category {
1061 Category::Infinity => {
1063 // nextUp(-inf) = -largest
1064 Status::OK.and(-Self::largest())
1066 // nextUp(+inf) = +inf
1067 Status::OK.and(self)
1071 // IEEE-754R 2008 6.2 Par 2: nextUp(sNaN) = qNaN. Set Invalid flag.
1072 // IEEE-754R 2008 6.2: nextUp(qNaN) = qNaN. Must be identity so we do not
1073 // change the payload.
1074 if self.is_signaling() {
1075 // For consistency, propagate the sign of the sNaN to the qNaN.
1076 Status::INVALID_OP.and(Self::NAN.copy_sign(self))
1078 Status::OK.and(self)
1082 // nextUp(pm 0) = +smallest
1083 Status::OK.and(Self::SMALLEST)
1085 Category::Normal => {
1086 // nextUp(-smallest) = -0
1087 if self.is_smallest() && self.sign {
1088 return Status::OK.and(-Self::ZERO);
1091 // nextUp(largest) == INFINITY
1092 if self.is_largest() && !self.sign {
1093 return Status::OK.and(Self::INFINITY);
1096 // Excluding the integral bit. This allows us to test for binade boundaries.
1097 let sig_mask = (1 << (S::PRECISION - 1)) - 1;
1099 // nextUp(normal) == normal + inc.
1101 // If we are negative, we need to decrement the significand.
1103 // We only cross a binade boundary that requires adjusting the exponent
1105 // 1. exponent != S::MIN_EXP. This implies we are not in the
1106 // smallest binade or are dealing with denormals.
1107 // 2. Our significand excluding the integral bit is all zeros.
1108 let crossing_binade_boundary =
1109 self.exp != S::MIN_EXP && self.sig[0] & sig_mask == 0;
1111 // Decrement the significand.
1113 // We always do this since:
1114 // 1. If we are dealing with a non-binade decrement, by definition we
1115 // just decrement the significand.
1116 // 2. If we are dealing with a normal -> normal binade decrement, since
1117 // we have an explicit integral bit the fact that all bits but the
1118 // integral bit are zero implies that subtracting one will yield a
1119 // significand with 0 integral bit and 1 in all other spots. Thus we
1120 // must just adjust the exponent and set the integral bit to 1.
1121 // 3. If we are dealing with a normal -> denormal binade decrement,
1122 // since we set the integral bit to 0 when we represent denormals, we
1123 // just decrement the significand.
1124 sig::decrement(&mut self.sig);
1126 if crossing_binade_boundary {
1127 // Our result is a normal number. Do the following:
1128 // 1. Set the integral bit to 1.
1129 // 2. Decrement the exponent.
1130 sig::set_bit(&mut self.sig, S::PRECISION - 1);
1134 // If we are positive, we need to increment the significand.
1136 // We only cross a binade boundary that requires adjusting the exponent if
1137 // the input is not a denormal and all of said input's significand bits
1138 // are set. If all of said conditions are true: clear the significand, set
1139 // the integral bit to 1, and increment the exponent. If we have a
1140 // denormal always increment since moving denormals and the numbers in the
1141 // smallest normal binade have the same exponent in our representation.
1142 let crossing_binade_boundary =
1143 !self.is_denormal() && self.sig[0] & sig_mask == sig_mask;
1145 if crossing_binade_boundary {
1147 sig::set_bit(&mut self.sig, S::PRECISION - 1);
1151 "We can not increment an exponent beyond the MAX_EXP \
1152 allowed by the given floating point semantics."
1156 sig::increment(&mut self.sig);
1159 Status::OK.and(self)
1164 fn from_bits(input: u128) -> Self {
1165 // Dispatch to semantics.
1169 fn from_u128_r(input: u128, round: Round) -> StatusAnd<Self> {
1172 exp: S::PRECISION as ExpInt - 1,
1173 category: Category::Normal,
1175 marker: PhantomData,
1177 .normalize(round, Loss::ExactlyZero)
1180 fn from_str_r(mut s: &str, mut round: Round) -> Result<StatusAnd<Self>, ParseError> {
1182 return Err(ParseError("Invalid string length"));
1185 // Handle special cases.
1187 "inf" | "INFINITY" => return Ok(Status::OK.and(Self::INFINITY)),
1188 "-inf" | "-INFINITY" => return Ok(Status::OK.and(-Self::INFINITY)),
1189 "nan" | "NaN" => return Ok(Status::OK.and(Self::NAN)),
1190 "-nan" | "-NaN" => return Ok(Status::OK.and(-Self::NAN)),
1194 // Handle a leading minus sign.
1195 let minus = s.starts_with('-');
1196 if minus || s.starts_with('+') {
1199 return Err(ParseError("String has no digits"));
1203 // Adjust the rounding mode for the absolute value below.
1208 let r = if s.starts_with("0x") || s.starts_with("0X") {
1211 return Err(ParseError("Invalid string"));
1213 Self::from_hexadecimal_string(s, round)?
1215 Self::from_decimal_string(s, round)?
1218 Ok(r.map(|r| if minus { -r } else { r }))
1221 fn to_bits(self) -> u128 {
1222 // Dispatch to semantics.
1226 fn to_u128_r(self, width: usize, round: Round, is_exact: &mut bool) -> StatusAnd<u128> {
1227 // The result of trying to convert a number too large.
1228 let overflow = if self.sign {
1229 // Negative numbers cannot be represented as unsigned.
1232 // Largest unsigned integer of the given width.
1238 match self.category {
1239 Category::NaN => Status::INVALID_OP.and(0),
1241 Category::Infinity => Status::INVALID_OP.and(overflow),
1244 // Negative zero can't be represented as an int.
1245 *is_exact = !self.sign;
1249 Category::Normal => {
1252 // Step 1: place our absolute value, with any fraction truncated, in
1254 let truncated_bits = if self.exp < 0 {
1255 // Our absolute value is less than one; truncate everything.
1256 // For exponent -1 the integer bit represents .5, look at that.
1257 // For smaller exponents leftmost truncated bit is 0.
1258 S::PRECISION - 1 + (-self.exp) as usize
1260 // We want the most significant (exponent + 1) bits; the rest are
1262 let bits = self.exp as usize + 1;
1264 // Hopelessly large in magnitude?
1266 return Status::INVALID_OP.and(overflow);
1269 if bits < S::PRECISION {
1270 // We truncate (S::PRECISION - bits) bits.
1271 r = self.sig[0] >> (S::PRECISION - bits);
1274 // We want at least as many bits as are available.
1275 r = self.sig[0] << (bits - S::PRECISION);
1280 // Step 2: work out any lost fraction, and increment the absolute
1281 // value if we would round away from zero.
1282 let mut loss = Loss::ExactlyZero;
1283 if truncated_bits > 0 {
1284 loss = Loss::through_truncation(&self.sig, truncated_bits);
1285 if loss != Loss::ExactlyZero
1286 && self.round_away_from_zero(round, loss, truncated_bits)
1288 r = r.wrapping_add(1);
1290 return Status::INVALID_OP.and(overflow); // Overflow.
1295 // Step 3: check if we fit in the destination.
1297 return Status::INVALID_OP.and(overflow);
1300 if loss == Loss::ExactlyZero {
1304 Status::INEXACT.and(r)
1310 fn cmp_abs_normal(self, rhs: Self) -> Ordering {
1311 assert!(self.is_finite_non_zero());
1312 assert!(rhs.is_finite_non_zero());
1314 // If exponents are equal, do an unsigned comparison of the significands.
1315 self.exp.cmp(&rhs.exp).then_with(|| sig::cmp(&self.sig, &rhs.sig))
1318 fn bitwise_eq(self, rhs: Self) -> bool {
1319 if self.category != rhs.category || self.sign != rhs.sign {
1323 if self.category == Category::Zero || self.category == Category::Infinity {
1327 if self.is_finite_non_zero() && self.exp != rhs.exp {
1334 fn is_negative(self) -> bool {
1338 fn is_denormal(self) -> bool {
1339 self.is_finite_non_zero()
1340 && self.exp == S::MIN_EXP
1341 && !sig::get_bit(&self.sig, S::PRECISION - 1)
1344 fn is_signaling(self) -> bool {
1345 // IEEE-754R 2008 6.2.1: A signaling NaN bit string should be encoded with the
1346 // first bit of the trailing significand being 0.
1347 self.is_nan() && !sig::get_bit(&self.sig, S::QNAN_BIT)
1350 fn category(self) -> Category {
1354 fn get_exact_inverse(self) -> Option<Self> {
1355 // Special floats and denormals have no exact inverse.
1356 if !self.is_finite_non_zero() {
1360 // Check that the number is a power of two by making sure that only the
1361 // integer bit is set in the significand.
1362 if self.sig != [1 << (S::PRECISION - 1)] {
1367 let mut reciprocal = Self::from_u128(1).value;
1369 reciprocal = unpack!(status=, reciprocal / self);
1370 if status != Status::OK {
1374 // Avoid multiplication with a denormal, it is not safe on all platforms and
1375 // may be slower than a normal division.
1376 if reciprocal.is_denormal() {
1380 assert!(reciprocal.is_finite_non_zero());
1381 assert_eq!(reciprocal.sig, [1 << (S::PRECISION - 1)]);
1386 fn ilogb(mut self) -> ExpInt {
1393 if self.is_infinite() {
1396 if !self.is_denormal() {
1400 let sig_bits = (S::PRECISION - 1) as ExpInt;
1401 self.exp += sig_bits;
1402 self = self.normalize(Round::NearestTiesToEven, Loss::ExactlyZero).value;
1406 fn scalbn_r(mut self, exp: ExpInt, round: Round) -> Self {
1407 // If exp is wildly out-of-scale, simply adding it to self.exp will
1408 // overflow; clamp it to a safe range before adding, but ensure that the range
1409 // is large enough that the clamp does not change the result. The range we
1410 // need to support is the difference between the largest possible exponent and
1411 // the normalized exponent of half the smallest denormal.
1413 let sig_bits = (S::PRECISION - 1) as i32;
1414 let max_change = S::MAX_EXP as i32 - (S::MIN_EXP as i32 - sig_bits) + 1;
1416 // Clamp to one past the range ends to let normalize handle overflow.
1417 let exp_change = cmp::min(cmp::max(exp as i32, -max_change - 1), max_change);
1418 self.exp = self.exp.saturating_add(exp_change as ExpInt);
1419 self = self.normalize(round, Loss::ExactlyZero).value;
1421 sig::set_bit(&mut self.sig, S::QNAN_BIT);
1426 fn frexp_r(mut self, exp: &mut ExpInt, round: Round) -> Self {
1427 *exp = self.ilogb();
1429 // Quiet signalling nans.
1430 if *exp == IEK_NAN {
1431 sig::set_bit(&mut self.sig, S::QNAN_BIT);
1435 if *exp == IEK_INF {
1439 // 1 is added because frexp is defined to return a normalized fraction in
1440 // +/-[0.5, 1.0), rather than the usual +/-[1.0, 2.0).
1441 if *exp == IEK_ZERO {
1446 self.scalbn_r(-*exp, round)
1450 impl<S: Semantics, T: Semantics> FloatConvert<IeeeFloat<T>> for IeeeFloat<S> {
1451 fn convert_r(self, round: Round, loses_info: &mut bool) -> StatusAnd<IeeeFloat<T>> {
1452 let mut r = IeeeFloat {
1455 category: self.category,
1457 marker: PhantomData,
1460 // x86 has some unusual NaNs which cannot be represented in any other
1461 // format; note them here.
1462 fn is_x87_double_extended<S: Semantics>() -> bool {
1463 S::QNAN_SIGNIFICAND == X87DoubleExtendedS::QNAN_SIGNIFICAND
1465 let x87_special_nan = is_x87_double_extended::<S>()
1466 && !is_x87_double_extended::<T>()
1467 && r.category == Category::NaN
1468 && (r.sig[0] & S::QNAN_SIGNIFICAND) != S::QNAN_SIGNIFICAND;
1470 // If this is a truncation of a denormal number, and the target semantics
1471 // has larger exponent range than the source semantics (this can happen
1472 // when truncating from PowerPC double-double to double format), the
1473 // right shift could lose result mantissa bits. Adjust exponent instead
1474 // of performing excessive shift.
1475 let mut shift = T::PRECISION as ExpInt - S::PRECISION as ExpInt;
1476 if shift < 0 && r.is_finite_non_zero() {
1477 let mut exp_change = sig::omsb(&r.sig) as ExpInt - S::PRECISION as ExpInt;
1478 if r.exp + exp_change < T::MIN_EXP {
1479 exp_change = T::MIN_EXP - r.exp;
1481 if exp_change < shift {
1485 shift -= exp_change;
1486 r.exp += exp_change;
1490 // If this is a truncation, perform the shift.
1491 let loss = if shift < 0 && (r.is_finite_non_zero() || r.category == Category::NaN) {
1492 sig::shift_right(&mut r.sig, &mut 0, -shift as usize)
1497 // If this is an extension, perform the shift.
1498 if shift > 0 && (r.is_finite_non_zero() || r.category == Category::NaN) {
1499 sig::shift_left(&mut r.sig, &mut 0, shift as usize);
1503 if r.is_finite_non_zero() {
1504 r = unpack!(status=, r.normalize(round, loss));
1505 *loses_info = status != Status::OK;
1506 } else if r.category == Category::NaN {
1507 *loses_info = loss != Loss::ExactlyZero || x87_special_nan;
1509 // For x87 extended precision, we want to make a NaN, not a special NaN if
1510 // the input wasn't special either.
1511 if !x87_special_nan && is_x87_double_extended::<T>() {
1512 sig::set_bit(&mut r.sig, T::PRECISION - 1);
1515 // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1516 // does not give you back the same bits. This is dubious, and we
1517 // don't currently do it. You're really supposed to get
1518 // an invalid operation signal at runtime, but nobody does that.
1519 status = Status::OK;
1521 *loses_info = false;
1522 status = Status::OK;
1529 impl<S: Semantics> IeeeFloat<S> {
1530 /// Handle positive overflow. We either return infinity or
1531 /// the largest finite number. For negative overflow,
1532 /// negate the `round` argument before calling.
1533 fn overflow_result(round: Round) -> StatusAnd<Self> {
1536 Round::NearestTiesToEven | Round::NearestTiesToAway | Round::TowardPositive => {
1537 (Status::OVERFLOW | Status::INEXACT).and(Self::INFINITY)
1539 // Otherwise we become the largest finite number.
1540 Round::TowardNegative | Round::TowardZero => Status::INEXACT.and(Self::largest()),
1544 /// Returns `true` if, when truncating the current number, with `bit` the
1545 /// new LSB, with the given lost fraction and rounding mode, the result
1546 /// would need to be rounded away from zero (i.e., by increasing the
1547 /// signficand). This routine must work for `Category::Zero` of both signs, and
1548 /// `Category::Normal` numbers.
1549 fn round_away_from_zero(&self, round: Round, loss: Loss, bit: usize) -> bool {
1550 // NaNs and infinities should not have lost fractions.
1551 assert!(self.is_finite_non_zero() || self.is_zero());
1553 // Current callers never pass this so we don't handle it.
1554 assert_ne!(loss, Loss::ExactlyZero);
1557 Round::NearestTiesToAway => loss == Loss::ExactlyHalf || loss == Loss::MoreThanHalf,
1558 Round::NearestTiesToEven => {
1559 if loss == Loss::MoreThanHalf {
1563 // Our zeros don't have a significand to test.
1564 if loss == Loss::ExactlyHalf && self.category != Category::Zero {
1565 return sig::get_bit(&self.sig, bit);
1570 Round::TowardZero => false,
1571 Round::TowardPositive => !self.sign,
1572 Round::TowardNegative => self.sign,
1576 fn normalize(mut self, round: Round, mut loss: Loss) -> StatusAnd<Self> {
1577 if !self.is_finite_non_zero() {
1578 return Status::OK.and(self);
1581 // Before rounding normalize the exponent of Category::Normal numbers.
1582 let mut omsb = sig::omsb(&self.sig);
1585 // OMSB is numbered from 1. We want to place it in the integer
1586 // bit numbered PRECISION if possible, with a compensating change in
1588 let mut final_exp = self.exp.saturating_add(omsb as ExpInt - S::PRECISION as ExpInt);
1590 // If the resulting exponent is too high, overflow according to
1591 // the rounding mode.
1592 if final_exp > S::MAX_EXP {
1593 let round = if self.sign { -round } else { round };
1594 return Self::overflow_result(round).map(|r| r.copy_sign(self));
1597 // Subnormal numbers have exponent MIN_EXP, and their MSB
1598 // is forced based on that.
1599 if final_exp < S::MIN_EXP {
1600 final_exp = S::MIN_EXP;
1603 // Shifting left is easy as we don't lose precision.
1604 if final_exp < self.exp {
1605 assert_eq!(loss, Loss::ExactlyZero);
1607 let exp_change = (self.exp - final_exp) as usize;
1608 sig::shift_left(&mut self.sig, &mut self.exp, exp_change);
1610 return Status::OK.and(self);
1613 // Shift right and capture any new lost fraction.
1614 if final_exp > self.exp {
1615 let exp_change = (final_exp - self.exp) as usize;
1616 loss = sig::shift_right(&mut self.sig, &mut self.exp, exp_change).combine(loss);
1618 // Keep OMSB up-to-date.
1619 omsb = omsb.saturating_sub(exp_change);
1623 // Now round the number according to round given the lost
1626 // As specified in IEEE 754, since we do not trap we do not report
1627 // underflow for exact results.
1628 if loss == Loss::ExactlyZero {
1629 // Canonicalize zeros.
1631 self.category = Category::Zero;
1634 return Status::OK.and(self);
1637 // Increment the significand if we're rounding away from zero.
1638 if self.round_away_from_zero(round, loss, 0) {
1640 self.exp = S::MIN_EXP;
1643 // We should never overflow.
1644 assert_eq!(sig::increment(&mut self.sig), 0);
1645 omsb = sig::omsb(&self.sig);
1647 // Did the significand increment overflow?
1648 if omsb == S::PRECISION + 1 {
1649 // Renormalize by incrementing the exponent and shifting our
1650 // significand right one. However if we already have the
1651 // maximum exponent we overflow to infinity.
1652 if self.exp == S::MAX_EXP {
1653 self.category = Category::Infinity;
1655 return (Status::OVERFLOW | Status::INEXACT).and(self);
1658 let _: Loss = sig::shift_right(&mut self.sig, &mut self.exp, 1);
1660 return Status::INEXACT.and(self);
1664 // The normal case - we were and are not denormal, and any
1665 // significand increment above didn't overflow.
1666 if omsb == S::PRECISION {
1667 return Status::INEXACT.and(self);
1670 // We have a non-zero denormal.
1671 assert!(omsb < S::PRECISION);
1673 // Canonicalize zeros.
1675 self.category = Category::Zero;
1678 // The Category::Zero case is a denormal that underflowed to zero.
1679 (Status::UNDERFLOW | Status::INEXACT).and(self)
1682 fn from_hexadecimal_string(s: &str, round: Round) -> Result<StatusAnd<Self>, ParseError> {
1683 let mut r = IeeeFloat {
1686 category: Category::Normal,
1688 marker: PhantomData,
1691 let mut any_digits = false;
1692 let mut has_exp = false;
1693 let mut bit_pos = LIMB_BITS as isize;
1694 let mut loss = None;
1696 // Without leading or trailing zeros, irrespective of the dot.
1697 let mut first_sig_digit = None;
1698 let mut dot = s.len();
1700 for (p, c) in s.char_indices() {
1701 // Skip leading zeros and any (hexa)decimal point.
1704 return Err(ParseError("String contains multiple dots"));
1707 } else if let Some(hex_value) = c.to_digit(16) {
1710 if first_sig_digit.is_none() {
1714 first_sig_digit = Some(p);
1717 // Store the number while we have space.
1720 r.sig[0] |= (hex_value as Limb) << bit_pos;
1721 // If zero or one-half (the hexadecimal digit 8) are followed
1722 // by non-zero, they're a little more than zero or one-half.
1723 } else if let Some(ref mut loss) = loss {
1725 if *loss == Loss::ExactlyZero {
1726 *loss = Loss::LessThanHalf;
1728 if *loss == Loss::ExactlyHalf {
1729 *loss = Loss::MoreThanHalf;
1733 loss = Some(match hex_value {
1734 0 => Loss::ExactlyZero,
1735 1..=7 => Loss::LessThanHalf,
1736 8 => Loss::ExactlyHalf,
1737 9..=15 => Loss::MoreThanHalf,
1738 _ => unreachable!(),
1741 } else if c == 'p' || c == 'P' {
1743 return Err(ParseError("Significand has no digits"));
1750 let mut chars = s[p + 1..].chars().peekable();
1752 // Adjust for the given exponent.
1753 let exp_minus = chars.peek() == Some(&'-');
1754 if exp_minus || chars.peek() == Some(&'+') {
1759 if let Some(value) = c.to_digit(10) {
1761 r.exp = r.exp.saturating_mul(10).saturating_add(value as ExpInt);
1763 return Err(ParseError("Invalid character in exponent"));
1767 return Err(ParseError("Exponent has no digits"));
1776 return Err(ParseError("Invalid character in significand"));
1780 return Err(ParseError("Significand has no digits"));
1783 // Hex floats require an exponent but not a hexadecimal point.
1785 return Err(ParseError("Hex strings require an exponent"));
1788 // Ignore the exponent if we are zero.
1789 let first_sig_digit = match first_sig_digit {
1791 None => return Ok(Status::OK.and(Self::ZERO)),
1794 // Calculate the exponent adjustment implicit in the number of
1795 // significant digits and adjust for writing the significand starting
1796 // at the most significant nibble.
1797 let exp_adjustment = if dot > first_sig_digit {
1798 ExpInt::try_from(dot - first_sig_digit).unwrap()
1800 -ExpInt::try_from(first_sig_digit - dot - 1).unwrap()
1802 let exp_adjustment = exp_adjustment
1805 .saturating_add(S::PRECISION as ExpInt)
1806 .saturating_sub(LIMB_BITS as ExpInt);
1807 r.exp = r.exp.saturating_add(exp_adjustment);
1809 Ok(r.normalize(round, loss.unwrap_or(Loss::ExactlyZero)))
1812 fn from_decimal_string(s: &str, round: Round) -> Result<StatusAnd<Self>, ParseError> {
1813 // Given a normal decimal floating point number of the form
1815 // dddd.dddd[eE][+-]ddd
1817 // where the decimal point and exponent are optional, fill out the
1818 // variables below. Exponent is appropriate if the significand is
1819 // treated as an integer, and normalized_exp if the significand
1820 // is taken to have the decimal point after a single leading
1823 // If the value is zero, first_sig_digit is None.
1825 let mut any_digits = false;
1826 let mut dec_exp = 0i32;
1828 // Without leading or trailing zeros, irrespective of the dot.
1829 let mut first_sig_digit = None;
1830 let mut last_sig_digit = 0;
1831 let mut dot = s.len();
1833 for (p, c) in s.char_indices() {
1836 return Err(ParseError("String contains multiple dots"));
1839 } else if let Some(dec_value) = c.to_digit(10) {
1843 if first_sig_digit.is_none() {
1844 first_sig_digit = Some(p);
1848 } else if c == 'e' || c == 'E' {
1850 return Err(ParseError("Significand has no digits"));
1857 let mut chars = s[p + 1..].chars().peekable();
1859 // Adjust for the given exponent.
1860 let exp_minus = chars.peek() == Some(&'-');
1861 if exp_minus || chars.peek() == Some(&'+') {
1867 if let Some(value) = c.to_digit(10) {
1869 dec_exp = dec_exp.saturating_mul(10).saturating_add(value as i32);
1871 return Err(ParseError("Invalid character in exponent"));
1875 return Err(ParseError("Exponent has no digits"));
1884 return Err(ParseError("Invalid character in significand"));
1888 return Err(ParseError("Significand has no digits"));
1891 // Test if we have a zero number allowing for non-zero exponents.
1892 let first_sig_digit = match first_sig_digit {
1894 None => return Ok(Status::OK.and(Self::ZERO)),
1897 // Adjust the exponents for any decimal point.
1898 if dot > last_sig_digit {
1899 dec_exp = dec_exp.saturating_add((dot - last_sig_digit - 1) as i32);
1901 dec_exp = dec_exp.saturating_sub((last_sig_digit - dot) as i32);
1903 let significand_digits = last_sig_digit - first_sig_digit + 1
1904 - (dot > first_sig_digit && dot < last_sig_digit) as usize;
1905 let normalized_exp = dec_exp.saturating_add(significand_digits as i32 - 1);
1907 // Handle the cases where exponents are obviously too large or too
1908 // small. Writing L for log 10 / log 2, a number d.ddddd*10^dec_exp
1909 // definitely overflows if
1911 // (dec_exp - 1) * L >= MAX_EXP
1913 // and definitely underflows to zero where
1915 // (dec_exp + 1) * L <= MIN_EXP - PRECISION
1917 // With integer arithmetic the tightest bounds for L are
1919 // 93/28 < L < 196/59 [ numerator <= 256 ]
1920 // 42039/12655 < L < 28738/8651 [ numerator <= 65536 ]
1922 // Check for MAX_EXP.
1923 if normalized_exp.saturating_sub(1).saturating_mul(42039) >= 12655 * S::MAX_EXP as i32 {
1924 // Overflow and round.
1925 return Ok(Self::overflow_result(round));
1928 // Check for MIN_EXP.
1929 if normalized_exp.saturating_add(1).saturating_mul(28738)
1930 <= 8651 * (S::MIN_EXP as i32 - S::PRECISION as i32)
1932 // Underflow to zero and round.
1934 if round == Round::TowardPositive { IeeeFloat::SMALLEST } else { IeeeFloat::ZERO };
1935 return Ok((Status::UNDERFLOW | Status::INEXACT).and(r));
1938 // A tight upper bound on number of bits required to hold an
1939 // N-digit decimal integer is N * 196 / 59. Allocate enough space
1940 // to hold the full significand, and an extra limb required by
1942 let max_limbs = limbs_for_bits(1 + 196 * significand_digits / 59);
1943 let mut dec_sig: SmallVec<[Limb; 1]> = SmallVec::with_capacity(max_limbs);
1945 // Convert to binary efficiently - we do almost all multiplication
1946 // in a Limb. When this would overflow do we do a single
1947 // bignum multiplication, and then revert again to multiplication
1949 let mut chars = s[first_sig_digit..=last_sig_digit].chars();
1952 let mut multiplier = 1;
1955 let dec_value = match chars.next() {
1956 Some('.') => continue,
1957 Some(c) => c.to_digit(10).unwrap(),
1962 val = val * 10 + dec_value as Limb;
1964 // The maximum number that can be multiplied by ten with any
1965 // digit added without overflowing a Limb.
1966 if multiplier > (!0 - 9) / 10 {
1971 // If we've consumed no digits, we're done.
1972 if multiplier == 1 {
1976 // Multiply out the current limb.
1977 let mut carry = val;
1978 for x in &mut dec_sig {
1979 let [low, mut high] = sig::widening_mul(*x, multiplier);
1982 let (low, overflow) = low.overflowing_add(carry);
1983 high += overflow as Limb;
1989 // If we had carry, we need another limb (likely but not guaranteed).
1991 dec_sig.push(carry);
1995 // Calculate pow(5, abs(dec_exp)) into `pow5_full`.
1996 // The *_calc Vec's are reused scratch space, as an optimization.
1997 let (pow5_full, mut pow5_calc, mut sig_calc, mut sig_scratch_calc) = {
1998 let mut power = dec_exp.abs() as usize;
2000 const FIRST_EIGHT_POWERS: [Limb; 8] = [1, 5, 25, 125, 625, 3125, 15625, 78125];
2002 let mut p5_scratch = smallvec![];
2003 let mut p5: SmallVec<[Limb; 1]> = smallvec![FIRST_EIGHT_POWERS[4]];
2005 let mut r_scratch = smallvec![];
2006 let mut r: SmallVec<[Limb; 1]> = smallvec![FIRST_EIGHT_POWERS[power & 7]];
2010 // Calculate pow(5,pow(2,n+3)).
2011 p5_scratch.resize(p5.len() * 2, 0);
2012 let _: Loss = sig::mul(&mut p5_scratch, &mut 0, &p5, &p5, p5.len() * 2 * LIMB_BITS);
2013 while p5_scratch.last() == Some(&0) {
2016 mem::swap(&mut p5, &mut p5_scratch);
2019 r_scratch.resize(r.len() + p5.len(), 0);
2021 sig::mul(&mut r_scratch, &mut 0, &r, &p5, (r.len() + p5.len()) * LIMB_BITS);
2022 while r_scratch.last() == Some(&0) {
2025 mem::swap(&mut r, &mut r_scratch);
2031 (r, r_scratch, p5, p5_scratch)
2034 // Attempt dec_sig * 10^dec_exp with increasing precision.
2035 let mut attempt = 0;
2037 let calc_precision = (LIMB_BITS << attempt) - 1;
2040 let calc_normal_from_limbs = |sig: &mut SmallVec<[Limb; 1]>,
2042 -> StatusAnd<ExpInt> {
2043 sig.resize(limbs_for_bits(calc_precision), 0);
2044 let (mut loss, mut exp) = sig::from_limbs(sig, limbs, calc_precision);
2046 // Before rounding normalize the exponent of Category::Normal numbers.
2047 let mut omsb = sig::omsb(sig);
2049 assert_ne!(omsb, 0);
2051 // OMSB is numbered from 1. We want to place it in the integer
2052 // bit numbered PRECISION if possible, with a compensating change in
2054 let final_exp = exp.saturating_add(omsb as ExpInt - calc_precision as ExpInt);
2056 // Shifting left is easy as we don't lose precision.
2057 if final_exp < exp {
2058 assert_eq!(loss, Loss::ExactlyZero);
2060 let exp_change = (exp - final_exp) as usize;
2061 sig::shift_left(sig, &mut exp, exp_change);
2063 return Status::OK.and(exp);
2066 // Shift right and capture any new lost fraction.
2067 if final_exp > exp {
2068 let exp_change = (final_exp - exp) as usize;
2069 loss = sig::shift_right(sig, &mut exp, exp_change).combine(loss);
2071 // Keep OMSB up-to-date.
2072 omsb = omsb.saturating_sub(exp_change);
2075 assert_eq!(omsb, calc_precision);
2077 // Now round the number according to round given the lost
2080 // As specified in IEEE 754, since we do not trap we do not report
2081 // underflow for exact results.
2082 if loss == Loss::ExactlyZero {
2083 return Status::OK.and(exp);
2086 // Increment the significand if we're rounding away from zero.
2087 if loss == Loss::MoreThanHalf || loss == Loss::ExactlyHalf && sig::get_bit(sig, 0) {
2088 // We should never overflow.
2089 assert_eq!(sig::increment(sig), 0);
2090 omsb = sig::omsb(sig);
2092 // Did the significand increment overflow?
2093 if omsb == calc_precision + 1 {
2094 let _: Loss = sig::shift_right(sig, &mut exp, 1);
2096 return Status::INEXACT.and(exp);
2100 // The normal case - we were and are not denormal, and any
2101 // significand increment above didn't overflow.
2102 Status::INEXACT.and(exp)
2106 let mut exp = unpack!(status=,
2107 calc_normal_from_limbs(&mut sig_calc, &dec_sig));
2109 let pow5_exp = unpack!(pow5_status=,
2110 calc_normal_from_limbs(&mut pow5_calc, &pow5_full));
2112 // Add dec_exp, as 10^n = 5^n * 2^n.
2113 exp += dec_exp as ExpInt;
2115 let mut used_bits = S::PRECISION;
2116 let mut truncated_bits = calc_precision - used_bits;
2118 let half_ulp_err1 = (status != Status::OK) as Limb;
2119 let (calc_loss, half_ulp_err2);
2123 sig_scratch_calc.resize(sig_calc.len() + pow5_calc.len(), 0);
2124 calc_loss = sig::mul(
2125 &mut sig_scratch_calc,
2131 mem::swap(&mut sig_calc, &mut sig_scratch_calc);
2133 half_ulp_err2 = (pow5_status != Status::OK) as Limb;
2137 sig_scratch_calc.resize(sig_calc.len(), 0);
2138 calc_loss = sig::div(
2139 &mut sig_scratch_calc,
2145 mem::swap(&mut sig_calc, &mut sig_scratch_calc);
2147 // Denormal numbers have less precision.
2148 if exp < S::MIN_EXP {
2149 truncated_bits += (S::MIN_EXP - exp) as usize;
2150 used_bits = calc_precision.saturating_sub(truncated_bits);
2152 // Extra half-ulp lost in reciprocal of exponent.
2154 2 * (pow5_status != Status::OK || calc_loss != Loss::ExactlyZero) as Limb;
2157 // Both sig::mul and sig::div return the
2158 // result with the integer bit set.
2159 assert!(sig::get_bit(&sig_calc, calc_precision - 1));
2161 // The error from the true value, in half-ulps, on multiplying two
2162 // floating point numbers, which differ from the value they
2163 // approximate by at most half_ulp_err1 and half_ulp_err2 half-ulps, is strictly less
2164 // than the returned value.
2166 // See "How to Read Floating Point Numbers Accurately" by William D Clinger.
2167 assert!(half_ulp_err1 < 2 || half_ulp_err2 < 2 || (half_ulp_err1 + half_ulp_err2 < 8));
2169 let inexact = (calc_loss != Loss::ExactlyZero) as Limb;
2170 let half_ulp_err = if half_ulp_err1 + half_ulp_err2 == 0 {
2171 inexact * 2 // <= inexact half-ulps.
2173 inexact + 2 * (half_ulp_err1 + half_ulp_err2)
2176 let ulps_from_boundary = {
2177 let bits = calc_precision - used_bits - 1;
2179 let i = bits / LIMB_BITS;
2180 let limb = sig_calc[i] & (!0 >> (LIMB_BITS - 1 - bits % LIMB_BITS));
2181 let boundary = match round {
2182 Round::NearestTiesToEven | Round::NearestTiesToAway => 1 << (bits % LIMB_BITS),
2186 let delta = limb.wrapping_sub(boundary);
2187 cmp::min(delta, delta.wrapping_neg())
2188 } else if limb == boundary {
2189 if !sig::is_all_zeros(&sig_calc[1..i]) {
2194 } else if limb == boundary.wrapping_sub(1) {
2195 if sig_calc[1..i].iter().any(|&x| x.wrapping_neg() != 1) {
2198 sig_calc[0].wrapping_neg()
2205 // Are we guaranteed to round correctly if we truncate?
2206 if ulps_from_boundary.saturating_mul(2) >= half_ulp_err {
2207 let mut r = IeeeFloat {
2210 category: Category::Normal,
2212 marker: PhantomData,
2214 sig::extract(&mut r.sig, &sig_calc, used_bits, calc_precision - used_bits);
2215 // If we extracted less bits above we must adjust our exponent
2216 // to compensate for the implicit right shift.
2217 r.exp += (S::PRECISION - used_bits) as ExpInt;
2218 let loss = Loss::through_truncation(&sig_calc, truncated_bits);
2219 return Ok(r.normalize(round, loss));
2226 /// Combine the effect of two lost fractions.
2227 fn combine(self, less_significant: Loss) -> Loss {
2228 let mut more_significant = self;
2229 if less_significant != Loss::ExactlyZero {
2230 if more_significant == Loss::ExactlyZero {
2231 more_significant = Loss::LessThanHalf;
2232 } else if more_significant == Loss::ExactlyHalf {
2233 more_significant = Loss::MoreThanHalf;
2240 /// Returns the fraction lost were a bignum truncated losing the least
2241 /// significant `bits` bits.
2242 fn through_truncation(limbs: &[Limb], bits: usize) -> Loss {
2244 return Loss::ExactlyZero;
2247 let half_bit = bits - 1;
2248 let half_limb = half_bit / LIMB_BITS;
2249 let (half_limb, rest) = if half_limb < limbs.len() {
2250 (limbs[half_limb], &limbs[..half_limb])
2254 let half = 1 << (half_bit % LIMB_BITS);
2255 let has_half = half_limb & half != 0;
2256 let has_rest = half_limb & (half - 1) != 0 || !sig::is_all_zeros(rest);
2258 match (has_half, has_rest) {
2259 (false, false) => Loss::ExactlyZero,
2260 (false, true) => Loss::LessThanHalf,
2261 (true, false) => Loss::ExactlyHalf,
2262 (true, true) => Loss::MoreThanHalf,
2267 /// Implementation details of IeeeFloat significands, such as big integer arithmetic.
2268 /// As a rule of thumb, no functions in this module should dynamically allocate.
2270 use super::{limbs_for_bits, ExpInt, Limb, Loss, LIMB_BITS};
2271 use core::cmp::Ordering;
2274 pub(super) fn is_all_zeros(limbs: &[Limb]) -> bool {
2275 limbs.iter().all(|&l| l == 0)
2278 /// One, not zero, based LSB. That is, returns 0 for a zeroed significand.
2279 pub(super) fn olsb(limbs: &[Limb]) -> usize {
2283 .find(|(_, &limb)| limb != 0)
2284 .map_or(0, |(i, limb)| i * LIMB_BITS + limb.trailing_zeros() as usize + 1)
2287 /// One, not zero, based MSB. That is, returns 0 for a zeroed significand.
2288 pub(super) fn omsb(limbs: &[Limb]) -> usize {
2292 .rfind(|(_, &limb)| limb != 0)
2293 .map_or(0, |(i, limb)| (i + 1) * LIMB_BITS - limb.leading_zeros() as usize)
2296 /// Comparison (unsigned) of two significands.
2297 pub(super) fn cmp(a: &[Limb], b: &[Limb]) -> Ordering {
2298 assert_eq!(a.len(), b.len());
2299 for (a, b) in a.iter().zip(b).rev() {
2301 Ordering::Equal => {}
2309 /// Extracts the given bit.
2310 pub(super) fn get_bit(limbs: &[Limb], bit: usize) -> bool {
2311 limbs[bit / LIMB_BITS] & (1 << (bit % LIMB_BITS)) != 0
2314 /// Sets the given bit.
2315 pub(super) fn set_bit(limbs: &mut [Limb], bit: usize) {
2316 limbs[bit / LIMB_BITS] |= 1 << (bit % LIMB_BITS);
2319 /// Clear the given bit.
2320 pub(super) fn clear_bit(limbs: &mut [Limb], bit: usize) {
2321 limbs[bit / LIMB_BITS] &= !(1 << (bit % LIMB_BITS));
2324 /// Shifts `dst` left `bits` bits, subtract `bits` from its exponent.
2325 pub(super) fn shift_left(dst: &mut [Limb], exp: &mut ExpInt, bits: usize) {
2327 // Our exponent should not underflow.
2328 *exp = exp.checked_sub(bits as ExpInt).unwrap();
2330 // Jump is the inter-limb jump; shift is the intra-limb shift.
2331 let jump = bits / LIMB_BITS;
2332 let shift = bits % LIMB_BITS;
2334 for i in (0..dst.len()).rev() {
2340 // dst[i] comes from the two limbs src[i - jump] and, if we have
2341 // an intra-limb shift, src[i - jump - 1].
2342 limb = dst[i - jump];
2346 limb |= dst[i - jump - 1] >> (LIMB_BITS - shift);
2356 /// Shifts `dst` right `bits` bits noting lost fraction.
2357 pub(super) fn shift_right(dst: &mut [Limb], exp: &mut ExpInt, bits: usize) -> Loss {
2358 let loss = Loss::through_truncation(dst, bits);
2361 // Our exponent should not overflow.
2362 *exp = exp.checked_add(bits as ExpInt).unwrap();
2364 // Jump is the inter-limb jump; shift is the intra-limb shift.
2365 let jump = bits / LIMB_BITS;
2366 let shift = bits % LIMB_BITS;
2368 // Perform the shift. This leaves the most significant `bits` bits
2369 // of the result at zero.
2370 for i in 0..dst.len() {
2373 if i + jump >= dst.len() {
2376 limb = dst[i + jump];
2379 if i + jump + 1 < dst.len() {
2380 limb |= dst[i + jump + 1] << (LIMB_BITS - shift);
2392 /// Copies the bit vector of width `src_bits` from `src`, starting at bit SRC_LSB,
2393 /// to `dst`, such that the bit SRC_LSB becomes the least significant bit of `dst`.
2394 /// All high bits above `src_bits` in `dst` are zero-filled.
2395 pub(super) fn extract(dst: &mut [Limb], src: &[Limb], src_bits: usize, src_lsb: usize) {
2400 let dst_limbs = limbs_for_bits(src_bits);
2401 assert!(dst_limbs <= dst.len());
2403 let src = &src[src_lsb / LIMB_BITS..];
2404 dst[..dst_limbs].copy_from_slice(&src[..dst_limbs]);
2406 let shift = src_lsb % LIMB_BITS;
2407 let _: Loss = shift_right(&mut dst[..dst_limbs], &mut 0, shift);
2409 // We now have (dst_limbs * LIMB_BITS - shift) bits from `src`
2410 // in `dst`. If this is less that src_bits, append the rest, else
2411 // clear the high bits.
2412 let n = dst_limbs * LIMB_BITS - shift;
2414 let mask = (1 << (src_bits - n)) - 1;
2415 dst[dst_limbs - 1] |= (src[dst_limbs] & mask) << (n % LIMB_BITS);
2416 } else if n > src_bits && src_bits % LIMB_BITS > 0 {
2417 dst[dst_limbs - 1] &= (1 << (src_bits % LIMB_BITS)) - 1;
2420 // Clear high limbs.
2421 for x in &mut dst[dst_limbs..] {
2426 /// We want the most significant PRECISION bits of `src`. There may not
2427 /// be that many; extract what we can.
2428 pub(super) fn from_limbs(dst: &mut [Limb], src: &[Limb], precision: usize) -> (Loss, ExpInt) {
2429 let omsb = omsb(src);
2431 if precision <= omsb {
2432 extract(dst, src, precision, omsb - precision);
2433 (Loss::through_truncation(src, omsb - precision), omsb as ExpInt - 1)
2435 extract(dst, src, omsb, 0);
2436 (Loss::ExactlyZero, precision as ExpInt - 1)
2440 /// For every consecutive chunk of `bits` bits from `limbs`,
2441 /// going from most significant to the least significant bits,
2442 /// call `f` to transform those bits and store the result back.
2443 pub(super) fn each_chunk<F: FnMut(Limb) -> Limb>(limbs: &mut [Limb], bits: usize, mut f: F) {
2444 assert_eq!(LIMB_BITS % bits, 0);
2445 for limb in limbs.iter_mut().rev() {
2447 for i in (0..LIMB_BITS / bits).rev() {
2448 r |= f((*limb >> (i * bits)) & ((1 << bits) - 1)) << (i * bits);
2454 /// Increment in-place, return the carry flag.
2455 pub(super) fn increment(dst: &mut [Limb]) -> Limb {
2457 *x = x.wrapping_add(1);
2466 /// Decrement in-place, return the borrow flag.
2467 pub(super) fn decrement(dst: &mut [Limb]) -> Limb {
2469 *x = x.wrapping_sub(1);
2478 /// `a += b + c` where `c` is zero or one. Returns the carry flag.
2479 pub(super) fn add(a: &mut [Limb], b: &[Limb], mut c: Limb) -> Limb {
2482 for (a, &b) in a.iter_mut().zip(b) {
2483 let (r, overflow) = a.overflowing_add(b);
2484 let (r, overflow2) = r.overflowing_add(c);
2486 c = (overflow | overflow2) as Limb;
2492 /// `a -= b + c` where `c` is zero or one. Returns the borrow flag.
2493 pub(super) fn sub(a: &mut [Limb], b: &[Limb], mut c: Limb) -> Limb {
2496 for (a, &b) in a.iter_mut().zip(b) {
2497 let (r, overflow) = a.overflowing_sub(b);
2498 let (r, overflow2) = r.overflowing_sub(c);
2500 c = (overflow | overflow2) as Limb;
2506 /// `a += b` or `a -= b`. Does not preserve `b`.
2507 pub(super) fn add_or_sub(
2515 // Are we bigger exponent-wise than the RHS?
2516 let bits = *a_exp - b_exp;
2518 // Determine if the operation on the absolute values is effectively
2519 // an addition or subtraction.
2520 // Subtraction is more subtle than one might naively expect.
2521 if *a_sign ^ b_sign {
2522 let (reverse, loss);
2525 reverse = cmp(a_sig, b_sig) == Ordering::Less;
2526 loss = Loss::ExactlyZero;
2527 } else if bits > 0 {
2528 loss = shift_right(b_sig, &mut 0, (bits - 1) as usize);
2529 shift_left(a_sig, a_exp, 1);
2532 loss = shift_right(a_sig, a_exp, (-bits - 1) as usize);
2533 shift_left(b_sig, &mut 0, 1);
2537 let borrow = (loss != Loss::ExactlyZero) as Limb;
2539 // The code above is intended to ensure that no borrow is necessary.
2540 assert_eq!(sub(b_sig, a_sig, borrow), 0);
2541 a_sig.copy_from_slice(b_sig);
2544 // The code above is intended to ensure that no borrow is necessary.
2545 assert_eq!(sub(a_sig, b_sig, borrow), 0);
2548 // Invert the lost fraction - it was on the RHS and subtracted.
2550 Loss::LessThanHalf => Loss::MoreThanHalf,
2551 Loss::MoreThanHalf => Loss::LessThanHalf,
2555 let loss = if bits > 0 {
2556 shift_right(b_sig, &mut 0, bits as usize)
2558 shift_right(a_sig, a_exp, -bits as usize)
2560 // We have a guard bit; generating a carry cannot happen.
2561 assert_eq!(add(a_sig, b_sig, 0), 0);
2566 /// `[low, high] = a * b`.
2568 /// This cannot overflow, because
2570 /// `(n - 1) * (n - 1) + 2 * (n - 1) == (n - 1) * (n + 1)`
2572 /// which is less than n<sup>2</sup>.
2573 pub(super) fn widening_mul(a: Limb, b: Limb) -> [Limb; 2] {
2574 let mut wide = [0, 0];
2576 if a == 0 || b == 0 {
2580 const HALF_BITS: usize = LIMB_BITS / 2;
2582 let select = |limb, i| (limb >> (i * HALF_BITS)) & ((1 << HALF_BITS) - 1);
2585 let mut x = [select(a, i) * select(b, j), 0];
2586 shift_left(&mut x, &mut 0, (i + j) * HALF_BITS);
2587 assert_eq!(add(&mut wide, &x, 0), 0);
2594 /// `dst = a * b` (for normal `a` and `b`). Returns the lost fraction.
2595 pub(super) fn mul<'a>(
2602 // Put the narrower number on the `a` for less loops below.
2603 if a.len() > b.len() {
2604 mem::swap(&mut a, &mut b);
2607 for x in &mut dst[..b.len()] {
2611 for i in 0..a.len() {
2613 for j in 0..b.len() {
2614 let [low, mut high] = widening_mul(a[i], b[j]);
2617 let (low, overflow) = low.overflowing_add(carry);
2618 high += overflow as Limb;
2620 // And now `dst[i + j]`, and store the new low part there.
2621 let (low, overflow) = low.overflowing_add(dst[i + j]);
2622 high += overflow as Limb;
2627 dst[i + b.len()] = carry;
2630 // Assume the operands involved in the multiplication are single-precision
2631 // FP, and the two multiplicants are:
2632 // a = a23 . a22 ... a0 * 2^e1
2633 // b = b23 . b22 ... b0 * 2^e2
2634 // the result of multiplication is:
2635 // dst = c48 c47 c46 . c45 ... c0 * 2^(e1+e2)
2636 // Note that there are three significant bits at the left-hand side of the
2637 // radix point: two for the multiplication, and an overflow bit for the
2638 // addition (that will always be zero at this point). Move the radix point
2639 // toward left by two bits, and adjust exponent accordingly.
2642 // Convert the result having "2 * precision" significant-bits back to the one
2643 // having "precision" significant-bits. First, move the radix point from
2644 // poision "2*precision - 1" to "precision - 1". The exponent need to be
2645 // adjusted by "2*precision - 1" - "precision - 1" = "precision".
2646 *exp -= precision as ExpInt + 1;
2648 // In case MSB resides at the left-hand side of radix point, shift the
2649 // mantissa right by some amount to make sure the MSB reside right before
2650 // the radix point (i.e., "MSB . rest-significant-bits").
2652 // Note that the result is not normalized when "omsb < precision". So, the
2653 // caller needs to call IeeeFloat::normalize() if normalized value is
2655 let omsb = omsb(dst);
2656 if omsb <= precision { Loss::ExactlyZero } else { shift_right(dst, exp, omsb - precision) }
2659 /// `quotient = dividend / divisor`. Returns the lost fraction.
2660 /// Does not preserve `dividend` or `divisor`.
2662 quotient: &mut [Limb],
2664 dividend: &mut [Limb],
2665 divisor: &mut [Limb],
2668 // Normalize the divisor.
2669 let bits = precision - omsb(divisor);
2670 shift_left(divisor, &mut 0, bits);
2671 *exp += bits as ExpInt;
2673 // Normalize the dividend.
2674 let bits = precision - omsb(dividend);
2675 shift_left(dividend, exp, bits);
2678 let olsb_divisor = olsb(divisor);
2679 if olsb_divisor == precision {
2680 quotient.copy_from_slice(dividend);
2681 return Loss::ExactlyZero;
2684 // Ensure the dividend >= divisor initially for the loop below.
2685 // Incidentally, this means that the division loop below is
2686 // guaranteed to set the integer bit to one.
2687 if cmp(dividend, divisor) == Ordering::Less {
2688 shift_left(dividend, exp, 1);
2689 assert_ne!(cmp(dividend, divisor), Ordering::Less)
2692 // Helper for figuring out the lost fraction.
2693 let lost_fraction = |dividend: &[Limb], divisor: &[Limb]| match cmp(dividend, divisor) {
2694 Ordering::Greater => Loss::MoreThanHalf,
2695 Ordering::Equal => Loss::ExactlyHalf,
2697 if is_all_zeros(dividend) {
2705 // Try to perform a (much faster) short division for small divisors.
2706 let divisor_bits = precision - (olsb_divisor - 1);
2707 macro_rules! try_short_div {
2708 ($W:ty, $H:ty, $half:expr) => {
2709 if divisor_bits * 2 <= $half {
2710 // Extract the small divisor.
2711 let _: Loss = shift_right(divisor, &mut 0, olsb_divisor - 1);
2712 let divisor = divisor[0] as $H as $W;
2714 // Shift the dividend to produce a quotient with the unit bit set.
2715 let top_limb = *dividend.last().unwrap();
2716 let mut rem = (top_limb >> (LIMB_BITS - (divisor_bits - 1))) as $H;
2717 shift_left(dividend, &mut 0, divisor_bits - 1);
2719 // Apply short division in place on $H (of $half bits) chunks.
2720 each_chunk(dividend, $half, |chunk| {
2721 let chunk = chunk as $H;
2722 let combined = ((rem as $W) << $half) | (chunk as $W);
2723 rem = (combined % divisor) as $H;
2724 (combined / divisor) as $H as Limb
2726 quotient.copy_from_slice(dividend);
2728 return lost_fraction(&[(rem as Limb) << 1], &[divisor as Limb]);
2733 try_short_div!(u32, u16, 16);
2734 try_short_div!(u64, u32, 32);
2735 try_short_div!(u128, u64, 64);
2737 // Zero the quotient before setting bits in it.
2738 for x in &mut quotient[..limbs_for_bits(precision)] {
2743 for bit in (0..precision).rev() {
2744 if cmp(dividend, divisor) != Ordering::Less {
2745 sub(dividend, divisor, 0);
2746 set_bit(quotient, bit);
2748 shift_left(dividend, &mut 0, 1);
2751 lost_fraction(dividend, divisor)