1 // Copyright 2015 The Rust Project Developers. See the COPYRIGHT
2 // file at the top-level directory of this distribution and at
3 // http://rust-lang.org/COPYRIGHT.
5 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
6 // http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
7 // <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
8 // option. This file may not be copied, modified, or distributed
9 // except according to those terms.
11 //! Bit fiddling on positive IEEE 754 floats. Negative numbers aren't and needn't be handled.
12 //! Normal floating point numbers have a canonical representation as (frac, exp) such that the
13 //! value is 2<sup>exp</sup> * (1 + sum(frac[N-i] / 2<sup>i</sup>)) where N is the number of bits.
14 //! Subnormals are slightly different and weird, but the same principle applies.
16 //! Here, however, we represent them as (sig, k) with f positive, such that the value is f *
17 //! 2<sup>e</sup>. Besides making the "hidden bit" explicit, this changes the exponent by the
18 //! so-called mantissa shift.
20 //! Put another way, normally floats are written as (1) but here they are written as (2):
22 //! 1. `1.101100...11 * 2^m`
23 //! 2. `1101100...11 * 2^n`
25 //! We call (1) the **fractional representation** and (2) the **integral representation**.
27 //! Many functions in this module only handle normal numbers. The dec2flt routines conservatively
28 //! take the universally-correct slow path (Algorithm M) for very small and very large numbers.
29 //! That algorithm needs only next_float() which does handle subnormals and zeros.
30 use cmp::Ordering::{Less, Equal, Greater};
31 use convert::{TryFrom, TryInto};
32 use ops::{Add, Mul, Div, Neg};
33 use fmt::{Debug, LowerExp};
34 use num::diy_float::Fp;
35 use num::FpCategory::{Infinite, Zero, Subnormal, Normal, Nan};
37 use num::dec2flt::num::{self, Big};
38 use num::dec2flt::table;
40 #[derive(Copy, Clone, Debug)]
47 pub fn new(sig: u64, k: i16) -> Self {
52 /// A helper trait to avoid duplicating basically all the conversion code for `f32` and `f64`.
54 /// See the parent module's doc comment for why this is necessary.
56 /// Should **never ever** be implemented for other types or be used outside the dec2flt module.
69 /// Type used by `to_bits` and `from_bits`.
70 type Bits: Add<Output = Self::Bits> + From<u8> + TryFrom<u64>;
72 /// Raw transmutation to integer.
73 fn to_bits(self) -> Self::Bits;
75 /// Raw transmutation from integer.
76 fn from_bits(v: Self::Bits) -> Self;
78 /// Returns the category that this number falls into.
79 fn classify(self) -> FpCategory;
81 /// Returns the mantissa, exponent and sign as integers.
82 fn integer_decode(self) -> (u64, i16, i8);
85 fn unpack(self) -> Unpacked;
87 /// Cast from a small integer that can be represented exactly. Panic if the integer can't be
88 /// represented, the other code in this module makes sure to never let that happen.
89 fn from_int(x: u64) -> Self;
91 /// Get the value 10<sup>e</sup> from a pre-computed table.
92 /// Panics for `e >= CEIL_LOG5_OF_MAX_SIG`.
93 fn short_fast_pow10(e: usize) -> Self;
95 /// What the name says. It's easier to hard code than juggling intrinsics and
96 /// hoping LLVM constant folds it.
97 const CEIL_LOG5_OF_MAX_SIG: i16;
99 // A conservative bound on the decimal digits of inputs that can't produce overflow or zero or
100 /// subnormals. Probably the decimal exponent of the maximum normal value, hence the name.
101 const MAX_NORMAL_DIGITS: usize;
103 /// When the most significant decimal digit has a place value greater than this, the number
104 /// is certainly rounded to infinity.
105 const INF_CUTOFF: i64;
107 /// When the most significant decimal digit has a place value less than this, the number
108 /// is certainly rounded to zero.
109 const ZERO_CUTOFF: i64;
111 /// The number of bits in the exponent.
114 /// The number of bits in the significand, *including* the hidden bit.
117 /// The number of bits in the significand, *excluding* the hidden bit.
118 const EXPLICIT_SIG_BITS: u8;
120 /// The maximum legal exponent in fractional representation.
123 /// The minimum legal exponent in fractional representation, excluding subnormals.
126 /// `MAX_EXP` for integral representation, i.e., with the shift applied.
127 const MAX_EXP_INT: i16;
129 /// `MAX_EXP` encoded (i.e., with offset bias)
130 const MAX_ENCODED_EXP: i16;
132 /// `MIN_EXP` for integral representation, i.e., with the shift applied.
133 const MIN_EXP_INT: i16;
135 /// The maximum normalized significand in integral representation.
138 /// The minimal normalized significand in integral representation.
142 // Mostly a workaround for #34344.
143 macro_rules! other_constants {
145 const EXPLICIT_SIG_BITS: u8 = Self::SIG_BITS - 1;
146 const MAX_EXP: i16 = (1 << (Self::EXP_BITS - 1)) - 1;
147 const MIN_EXP: i16 = -Self::MAX_EXP + 1;
148 const MAX_EXP_INT: i16 = Self::MAX_EXP - (Self::SIG_BITS as i16 - 1);
149 const MAX_ENCODED_EXP: i16 = (1 << Self::EXP_BITS) - 1;
150 const MIN_EXP_INT: i16 = Self::MIN_EXP - (Self::SIG_BITS as i16 - 1);
151 const MAX_SIG: u64 = (1 << Self::SIG_BITS) - 1;
152 const MIN_SIG: u64 = 1 << (Self::SIG_BITS - 1);
154 const INFINITY: Self = $crate::$type::INFINITY;
155 const NAN: Self = $crate::$type::NAN;
156 const ZERO: Self = 0.0;
160 impl RawFloat for f32 {
163 const SIG_BITS: u8 = 24;
164 const EXP_BITS: u8 = 8;
165 const CEIL_LOG5_OF_MAX_SIG: i16 = 11;
166 const MAX_NORMAL_DIGITS: usize = 35;
167 const INF_CUTOFF: i64 = 40;
168 const ZERO_CUTOFF: i64 = -48;
169 other_constants!(f32);
171 /// Returns the mantissa, exponent and sign as integers.
172 fn integer_decode(self) -> (u64, i16, i8) {
173 let bits = self.to_bits();
174 let sign: i8 = if bits >> 31 == 0 { 1 } else { -1 };
175 let mut exponent: i16 = ((bits >> 23) & 0xff) as i16;
176 let mantissa = if exponent == 0 {
177 (bits & 0x7fffff) << 1
179 (bits & 0x7fffff) | 0x800000
181 // Exponent bias + mantissa shift
182 exponent -= 127 + 23;
183 (mantissa as u64, exponent, sign)
186 fn unpack(self) -> Unpacked {
187 let (sig, exp, _sig) = self.integer_decode();
188 Unpacked::new(sig, exp)
191 fn from_int(x: u64) -> f32 {
192 // rkruppe is uncertain whether `as` rounds correctly on all platforms.
193 debug_assert!(x as f32 == fp_to_float(Fp { f: x, e: 0 }));
197 fn short_fast_pow10(e: usize) -> Self {
198 table::F32_SHORT_POWERS[e]
201 fn classify(self) -> FpCategory { self.classify() }
202 fn to_bits(self) -> Self::Bits { self.to_bits() }
203 fn from_bits(v: Self::Bits) -> Self { Self::from_bits(v) }
207 impl RawFloat for f64 {
210 const SIG_BITS: u8 = 53;
211 const EXP_BITS: u8 = 11;
212 const CEIL_LOG5_OF_MAX_SIG: i16 = 23;
213 const MAX_NORMAL_DIGITS: usize = 305;
214 const INF_CUTOFF: i64 = 310;
215 const ZERO_CUTOFF: i64 = -326;
216 other_constants!(f64);
218 /// Returns the mantissa, exponent and sign as integers.
219 fn integer_decode(self) -> (u64, i16, i8) {
220 let bits = self.to_bits();
221 let sign: i8 = if bits >> 63 == 0 { 1 } else { -1 };
222 let mut exponent: i16 = ((bits >> 52) & 0x7ff) as i16;
223 let mantissa = if exponent == 0 {
224 (bits & 0xfffffffffffff) << 1
226 (bits & 0xfffffffffffff) | 0x10000000000000
228 // Exponent bias + mantissa shift
229 exponent -= 1023 + 52;
230 (mantissa, exponent, sign)
233 fn unpack(self) -> Unpacked {
234 let (sig, exp, _sig) = self.integer_decode();
235 Unpacked::new(sig, exp)
238 fn from_int(x: u64) -> f64 {
239 // rkruppe is uncertain whether `as` rounds correctly on all platforms.
240 debug_assert!(x as f64 == fp_to_float(Fp { f: x, e: 0 }));
244 fn short_fast_pow10(e: usize) -> Self {
245 table::F64_SHORT_POWERS[e]
248 fn classify(self) -> FpCategory { self.classify() }
249 fn to_bits(self) -> Self::Bits { self.to_bits() }
250 fn from_bits(v: Self::Bits) -> Self { Self::from_bits(v) }
253 /// Convert an Fp to the closest machine float type.
254 /// Does not handle subnormal results.
255 pub fn fp_to_float<T: RawFloat>(x: Fp) -> T {
256 let x = x.normalize();
257 // x.f is 64 bit, so x.e has a mantissa shift of 63
260 panic!("fp_to_float: exponent {} too large", e)
261 } else if e > T::MIN_EXP {
262 encode_normal(round_normal::<T>(x))
264 panic!("fp_to_float: exponent {} too small", e)
268 /// Round the 64-bit significand to T::SIG_BITS bits with half-to-even.
269 /// Does not handle exponent overflow.
270 pub fn round_normal<T: RawFloat>(x: Fp) -> Unpacked {
271 let excess = 64 - T::SIG_BITS as i16;
272 let half: u64 = 1 << (excess - 1);
273 let (q, rem) = (x.f >> excess, x.f & ((1 << excess) - 1));
274 assert_eq!(q << excess | rem, x.f);
275 // Adjust mantissa shift
276 let k = x.e + excess;
279 } else if rem == half && (q % 2) == 0 {
281 } else if q == T::MAX_SIG {
282 Unpacked::new(T::MIN_SIG, k + 1)
284 Unpacked::new(q + 1, k)
288 /// Inverse of `RawFloat::unpack()` for normalized numbers.
289 /// Panics if the significand or exponent are not valid for normalized numbers.
290 pub fn encode_normal<T: RawFloat>(x: Unpacked) -> T {
291 debug_assert!(T::MIN_SIG <= x.sig && x.sig <= T::MAX_SIG,
292 "encode_normal: significand not normalized");
293 // Remove the hidden bit
294 let sig_enc = x.sig & !(1 << T::EXPLICIT_SIG_BITS);
295 // Adjust the exponent for exponent bias and mantissa shift
296 let k_enc = x.k + T::MAX_EXP + T::EXPLICIT_SIG_BITS as i16;
297 debug_assert!(k_enc != 0 && k_enc < T::MAX_ENCODED_EXP,
298 "encode_normal: exponent out of range");
299 // Leave sign bit at 0 ("+"), our numbers are all positive
300 let bits = (k_enc as u64) << T::EXPLICIT_SIG_BITS | sig_enc;
301 T::from_bits(bits.try_into().unwrap_or_else(|_| unreachable!()))
304 /// Construct a subnormal. A mantissa of 0 is allowed and constructs zero.
305 pub fn encode_subnormal<T: RawFloat>(significand: u64) -> T {
306 assert!(significand < T::MIN_SIG, "encode_subnormal: not actually subnormal");
307 // Encoded exponent is 0, the sign bit is 0, so we just have to reinterpret the bits.
308 T::from_bits(significand.try_into().unwrap_or_else(|_| unreachable!()))
311 /// Approximate a bignum with an Fp. Rounds within 0.5 ULP with half-to-even.
312 pub fn big_to_fp(f: &Big) -> Fp {
313 let end = f.bit_length();
314 assert!(end != 0, "big_to_fp: unexpectedly, input is zero");
315 let start = end.saturating_sub(64);
316 let leading = num::get_bits(f, start, end);
317 // We cut off all bits prior to the index `start`, i.e., we effectively right-shift by
318 // an amount of `start`, so this is also the exponent we need.
319 let e = start as i16;
320 let rounded_down = Fp { f: leading, e }.normalize();
321 // Round (half-to-even) depending on the truncated bits.
322 match num::compare_with_half_ulp(f, start) {
323 Less => rounded_down,
324 Equal if leading % 2 == 0 => rounded_down,
325 Equal | Greater => match leading.checked_add(1) {
326 Some(f) => Fp { f, e }.normalize(),
327 None => Fp { f: 1 << 63, e: e + 1 },
332 /// Find the largest floating point number strictly smaller than the argument.
333 /// Does not handle subnormals, zero, or exponent underflow.
334 pub fn prev_float<T: RawFloat>(x: T) -> T {
336 Infinite => panic!("prev_float: argument is infinite"),
337 Nan => panic!("prev_float: argument is NaN"),
338 Subnormal => panic!("prev_float: argument is subnormal"),
339 Zero => panic!("prev_float: argument is zero"),
341 let Unpacked { sig, k } = x.unpack();
342 if sig == T::MIN_SIG {
343 encode_normal(Unpacked::new(T::MAX_SIG, k - 1))
345 encode_normal(Unpacked::new(sig - 1, k))
351 // Find the smallest floating point number strictly larger than the argument.
352 // This operation is saturating, i.e., next_float(inf) == inf.
353 // Unlike most code in this module, this function does handle zero, subnormals, and infinities.
354 // However, like all other code here, it does not deal with NaN and negative numbers.
355 pub fn next_float<T: RawFloat>(x: T) -> T {
357 Nan => panic!("next_float: argument is NaN"),
358 Infinite => T::INFINITY,
359 // This seems too good to be true, but it works.
360 // 0.0 is encoded as the all-zero word. Subnormals are 0x000m...m where m is the mantissa.
361 // In particular, the smallest subnormal is 0x0...01 and the largest is 0x000F...F.
362 // The smallest normal number is 0x0010...0, so this corner case works as well.
363 // If the increment overflows the mantissa, the carry bit increments the exponent as we
364 // want, and the mantissa bits become zero. Because of the hidden bit convention, this
365 // too is exactly what we want!
366 // Finally, f64::MAX + 1 = 7eff...f + 1 = 7ff0...0 = f64::INFINITY.
367 Zero | Subnormal | Normal => {
368 T::from_bits(x.to_bits() + T::Bits::from(1u8))