src/libcore/num/dec2flt/rawfp.rs

   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.
   4 //
   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.
  10
  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.
  15 //!
  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.
  19 //!
  20 //! Put another way, normally floats are written as (1) but here they are written as (2):
  21 //!
  22 //! 1. `1.101100...11 * 2^m`
  23 //! 2. `1101100...11 * 2^n`
  24 //!
  25 //! We call (1) the **fractional representation** and (2) the **integral representation**.
  26 //!
  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 u32;
  31 use cmp::Ordering::{Less, Equal, Greater};
  32 use ops::{Mul, Div, Neg};
  33 use fmt::{Debug, LowerExp};
  34 use mem::transmute;
  35 use num::diy_float::Fp;
  36 use num::FpCategory::{Infinite, Zero, Subnormal, Normal, Nan};
  37 use num::Float;
  38 use num::dec2flt::num::{self, Big};
  39 use num::dec2flt::table;
  40
  41 #[derive(Copy, Clone, Debug)]
  42 pub struct Unpacked {
  43     pub sig: u64,
  44     pub k: i16,
  45 }
  46
  47 impl Unpacked {
  48     pub fn new(sig: u64, k: i16) -> Self {
  49         Unpacked { sig: sig, k: k }
  50     }
  51 }
  52
  53 /// A helper trait to avoid duplicating basically all the conversion code for `f32` and `f64`.
  54 ///
  55 /// See the parent module's doc comment for why this is necessary.
  56 ///
  57 /// Should **never ever** be implemented for other types or be used outside the dec2flt module.
  58 /// Inherits from `Float` because there is some overlap, but all the reused methods are trivial.
  59 pub trait RawFloat : Float + Copy + Debug + LowerExp
  60                     + Mul<Output=Self> + Div<Output=Self> + Neg<Output=Self>
  61 {
  62     const INFINITY: Self;
  63     const NAN: Self;
  64     const ZERO: Self;
  65
  66     /// Returns the mantissa, exponent and sign as integers.
  67     fn integer_decode(self) -> (u64, i16, i8);
  68
  69     /// Get the raw binary representation of the float.
  70     fn transmute(self) -> u64;
  71
  72     /// Transmute the raw binary representation into a float.
  73     fn from_bits(bits: u64) -> Self;
  74
  75     /// Decode the float.
  76     fn unpack(self) -> Unpacked;
  77
  78     /// Cast from a small integer that can be represented exactly.  Panic if the integer can't be
  79     /// represented, the other code in this module makes sure to never let that happen.
  80     fn from_int(x: u64) -> Self;
  81
  82     /// Get the value 10<sup>e</sup> from a pre-computed table.
  83     /// Panics for `e >= CEIL_LOG5_OF_MAX_SIG`.
  84     fn short_fast_pow10(e: usize) -> Self;
  85
  86     /// What the name says. It's easier to hard code than juggling intrinsics and
  87     /// hoping LLVM constant folds it.
  88     const CEIL_LOG5_OF_MAX_SIG: i16;
  89
  90     // A conservative bound on the decimal digits of inputs that can't produce overflow or zero or
  91     /// subnormals. Probably the decimal exponent of the maximum normal value, hence the name.
  92     const MAX_NORMAL_DIGITS: usize;
  93
  94     /// When the most significant decimal digit has a place value greater than this, the number
  95     /// is certainly rounded to infinity.
  96     const INF_CUTOFF: i64;
  97
  98     /// When the most significant decimal digit has a place value less than this, the number
  99     /// is certainly rounded to zero.
 100     const ZERO_CUTOFF: i64;
 101
 102     /// The number of bits in the exponent.
 103     const EXP_BITS: u8;
 104
 105     /// The number of bits in the singificand, *including* the hidden bit.
 106     const SIG_BITS: u8;
 107
 108     /// The number of bits in the singificand, *excluding* the hidden bit.
 109     const EXPLICIT_SIG_BITS: u8;
 110
 111     /// The maximum legal exponent in fractional representation.
 112     const MAX_EXP: i16;
 113
 114     /// The minimum legal exponent in fractional representation, excluding subnormals.
 115     const MIN_EXP: i16;
 116
 117     /// `MAX_EXP` for integral representation, i.e., with the shift applied.
 118     const MAX_EXP_INT: i16;
 119
 120     /// `MAX_EXP` encoded (i.e., with offset bias)
 121     const MAX_ENCODED_EXP: i16;
 122
 123     /// `MIN_EXP` for integral representation, i.e., with the shift applied.
 124     const MIN_EXP_INT: i16;
 125
 126     /// The maximum normalized singificand in integral representation.
 127     const MAX_SIG: u64;
 128
 129     /// The minimal normalized significand in integral representation.
 130     const MIN_SIG: u64;
 131 }
 132
 133 // Mostly a workaround for #34344.
 134 macro_rules! other_constants {
 135     ($type: ident) => {
 136         const EXPLICIT_SIG_BITS: u8 = Self::SIG_BITS - 1;
 137         const MAX_EXP: i16 = (1 << (Self::EXP_BITS - 1)) - 1;
 138         const MIN_EXP: i16 = -Self::MAX_EXP + 1;
 139         const MAX_EXP_INT: i16 = Self::MAX_EXP - (Self::SIG_BITS as i16 - 1);
 140         const MAX_ENCODED_EXP: i16 = (1 << Self::EXP_BITS) - 1;
 141         const MIN_EXP_INT: i16 = Self::MIN_EXP - (Self::SIG_BITS as i16 - 1);
 142         const MAX_SIG: u64 = (1 << Self::SIG_BITS) - 1;
 143         const MIN_SIG: u64 = 1 << (Self::SIG_BITS - 1);
 144
 145         const INFINITY: Self = $crate::$type::INFINITY;
 146         const NAN: Self = $crate::$type::NAN;
 147         const ZERO: Self = 0.0;
 148     }
 149 }
 150
 151 impl RawFloat for f32 {
 152     const SIG_BITS: u8 = 24;
 153     const EXP_BITS: u8 = 8;
 154     const CEIL_LOG5_OF_MAX_SIG: i16 = 11;
 155     const MAX_NORMAL_DIGITS: usize = 35;
 156     const INF_CUTOFF: i64 = 40;
 157     const ZERO_CUTOFF: i64 = -48;
 158     other_constants!(f32);
 159
 160     /// Returns the mantissa, exponent and sign as integers.
 161     fn integer_decode(self) -> (u64, i16, i8) {
 162         let bits: u32 = unsafe { transmute(self) };
 163         let sign: i8 = if bits >> 31 == 0 { 1 } else { -1 };
 164         let mut exponent: i16 = ((bits >> 23) & 0xff) as i16;
 165         let mantissa = if exponent == 0 {
 166             (bits & 0x7fffff) << 1
 167         } else {
 168             (bits & 0x7fffff) | 0x800000
 169         };
 170         // Exponent bias + mantissa shift
 171         exponent -= 127 + 23;
 172         (mantissa as u64, exponent, sign)
 173     }
 174
 175     fn transmute(self) -> u64 {
 176         let bits: u32 = unsafe { transmute(self) };
 177         bits as u64
 178     }
 179
 180     fn from_bits(bits: u64) -> f32 {
 181         assert!(bits < u32::MAX as u64, "f32::from_bits: too many bits");
 182         unsafe { transmute(bits as u32) }
 183     }
 184
 185     fn unpack(self) -> Unpacked {
 186         let (sig, exp, _sig) = self.integer_decode();
 187         Unpacked::new(sig, exp)
 188     }
 189
 190     fn from_int(x: u64) -> f32 {
 191         // rkruppe is uncertain whether `as` rounds correctly on all platforms.
 192         debug_assert!(x as f32 == fp_to_float(Fp { f: x, e: 0 }));
 193         x as f32
 194     }
 195
 196     fn short_fast_pow10(e: usize) -> Self {
 197         table::F32_SHORT_POWERS[e]
 198     }
 199 }
 200
 201
 202 impl RawFloat for f64 {
 203     const SIG_BITS: u8 = 53;
 204     const EXP_BITS: u8 = 11;
 205     const CEIL_LOG5_OF_MAX_SIG: i16 = 23;
 206     const MAX_NORMAL_DIGITS: usize = 305;
 207     const INF_CUTOFF: i64 = 310;
 208     const ZERO_CUTOFF: i64 = -326;
 209     other_constants!(f64);
 210
 211     /// Returns the mantissa, exponent and sign as integers.
 212     fn integer_decode(self) -> (u64, i16, i8) {
 213         let bits: u64 = unsafe { transmute(self) };
 214         let sign: i8 = if bits >> 63 == 0 { 1 } else { -1 };
 215         let mut exponent: i16 = ((bits >> 52) & 0x7ff) as i16;
 216         let mantissa = if exponent == 0 {
 217             (bits & 0xfffffffffffff) << 1
 218         } else {
 219             (bits & 0xfffffffffffff) | 0x10000000000000
 220         };
 221         // Exponent bias + mantissa shift
 222         exponent -= 1023 + 52;
 223         (mantissa, exponent, sign)
 224     }
 225
 226     fn transmute(self) -> u64 {
 227         let bits: u64 = unsafe { transmute(self) };
 228         bits
 229     }
 230
 231     fn from_bits(bits: u64) -> f64 {
 232         unsafe { transmute(bits) }
 233     }
 234
 235     fn unpack(self) -> Unpacked {
 236         let (sig, exp, _sig) = self.integer_decode();
 237         Unpacked::new(sig, exp)
 238     }
 239
 240     fn from_int(x: u64) -> f64 {
 241         // rkruppe is uncertain whether `as` rounds correctly on all platforms.
 242         debug_assert!(x as f64 == fp_to_float(Fp { f: x, e: 0 }));
 243         x as f64
 244     }
 245
 246     fn short_fast_pow10(e: usize) -> Self {
 247         table::F64_SHORT_POWERS[e]
 248     }
 249 }
 250
 251 /// Convert an Fp to the closest machine float type.
 252 /// Does not handle subnormal results.
 253 pub fn fp_to_float<T: RawFloat>(x: Fp) -> T {
 254     let x = x.normalize();
 255     // x.f is 64 bit, so x.e has a mantissa shift of 63
 256     let e = x.e + 63;
 257     if e > T::MAX_EXP {
 258         panic!("fp_to_float: exponent {} too large", e)
 259     }  else if e > T::MIN_EXP {
 260         encode_normal(round_normal::<T>(x))
 261     } else {
 262         panic!("fp_to_float: exponent {} too small", e)
 263     }
 264 }
 265
 266 /// Round the 64-bit significand to T::SIG_BITS bits with half-to-even.
 267 /// Does not handle exponent overflow.
 268 pub fn round_normal<T: RawFloat>(x: Fp) -> Unpacked {
 269     let excess = 64 - T::SIG_BITS as i16;
 270     let half: u64 = 1 << (excess - 1);
 271     let (q, rem) = (x.f >> excess, x.f & ((1 << excess) - 1));
 272     assert_eq!(q << excess | rem, x.f);
 273     // Adjust mantissa shift
 274     let k = x.e + excess;
 275     if rem < half {
 276         Unpacked::new(q, k)
 277     } else if rem == half && (q % 2) == 0 {
 278         Unpacked::new(q, k)
 279     } else if q == T::MAX_SIG {
 280         Unpacked::new(T::MIN_SIG, k + 1)
 281     } else {
 282         Unpacked::new(q + 1, k)
 283     }
 284 }
 285
 286 /// Inverse of `RawFloat::unpack()` for normalized numbers.
 287 /// Panics if the significand or exponent are not valid for normalized numbers.
 288 pub fn encode_normal<T: RawFloat>(x: Unpacked) -> T {
 289     debug_assert!(T::MIN_SIG <= x.sig && x.sig <= T::MAX_SIG,
 290         "encode_normal: significand not normalized");
 291     // Remove the hidden bit
 292     let sig_enc = x.sig & !(1 << T::EXPLICIT_SIG_BITS);
 293     // Adjust the exponent for exponent bias and mantissa shift
 294     let k_enc = x.k + T::MAX_EXP + T::EXPLICIT_SIG_BITS as i16;
 295     debug_assert!(k_enc != 0 && k_enc < T::MAX_ENCODED_EXP,
 296         "encode_normal: exponent out of range");
 297     // Leave sign bit at 0 ("+"), our numbers are all positive
 298     let bits = (k_enc as u64) << T::EXPLICIT_SIG_BITS | sig_enc;
 299     T::from_bits(bits)
 300 }
 301
 302 /// Construct a subnormal. A mantissa of 0 is allowed and constructs zero.
 303 pub fn encode_subnormal<T: RawFloat>(significand: u64) -> T {
 304     assert!(significand < T::MIN_SIG, "encode_subnormal: not actually subnormal");
 305     // Encoded exponent is 0, the sign bit is 0, so we just have to reinterpret the bits.
 306     T::from_bits(significand)
 307 }
 308
 309 /// Approximate a bignum with an Fp. Rounds within 0.5 ULP with half-to-even.
 310 pub fn big_to_fp(f: &Big) -> Fp {
 311     let end = f.bit_length();
 312     assert!(end != 0, "big_to_fp: unexpectedly, input is zero");
 313     let start = end.saturating_sub(64);
 314     let leading = num::get_bits(f, start, end);
 315     // We cut off all bits prior to the index `start`, i.e., we effectively right-shift by
 316     // an amount of `start`, so this is also the exponent we need.
 317     let e = start as i16;
 318     let rounded_down = Fp { f: leading, e: e }.normalize();
 319     // Round (half-to-even) depending on the truncated bits.
 320     match num::compare_with_half_ulp(f, start) {
 321         Less => rounded_down,
 322         Equal if leading % 2 == 0 => rounded_down,
 323         Equal | Greater => match leading.checked_add(1) {
 324             Some(f) => Fp { f: f, e: e }.normalize(),
 325             None => Fp { f: 1 << 63, e: e + 1 },
 326         }
 327     }
 328 }
 329
 330 /// Find the largest floating point number strictly smaller than the argument.
 331 /// Does not handle subnormals, zero, or exponent underflow.
 332 pub fn prev_float<T: RawFloat>(x: T) -> T {
 333     match x.classify() {
 334         Infinite => panic!("prev_float: argument is infinite"),
 335         Nan => panic!("prev_float: argument is NaN"),
 336         Subnormal => panic!("prev_float: argument is subnormal"),
 337         Zero => panic!("prev_float: argument is zero"),
 338         Normal => {
 339             let Unpacked { sig, k } = x.unpack();
 340             if sig == T::MIN_SIG {
 341                 encode_normal(Unpacked::new(T::MAX_SIG, k - 1))
 342             } else {
 343                 encode_normal(Unpacked::new(sig - 1, k))
 344             }
 345         }
 346     }
 347 }
 348
 349 // Find the smallest floating point number strictly larger than the argument.
 350 // This operation is saturating, i.e. next_float(inf) == inf.
 351 // Unlike most code in this module, this function does handle zero, subnormals, and infinities.
 352 // However, like all other code here, it does not deal with NaN and negative numbers.
 353 pub fn next_float<T: RawFloat>(x: T) -> T {
 354     match x.classify() {
 355         Nan => panic!("next_float: argument is NaN"),
 356         Infinite => T::INFINITY,
 357         // This seems too good to be true, but it works.
 358         // 0.0 is encoded as the all-zero word. Subnormals are 0x000m...m where m is the mantissa.
 359         // In particular, the smallest subnormal is 0x0...01 and the largest is 0x000F...F.
 360         // The smallest normal number is 0x0010...0, so this corner case works as well.
 361         // If the increment overflows the mantissa, the carry bit increments the exponent as we
 362         // want, and the mantissa bits become zero. Because of the hidden bit convention, this
 363         // too is exactly what we want!
 364         // Finally, f64::MAX + 1 = 7eff...f + 1 = 7ff0...0 = f64::INFINITY.
 365         Zero | Subnormal | Normal => {
 366             let bits: u64 = x.transmute();
 367             T::from_bits(bits + 1)
 368         }
 369     }
 370 }