1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
// Copyright 2015 The Rust Project Developers. See the COPYRIGHT
// file at the top-level directory of this distribution and at
// http://rust-lang.org/COPYRIGHT.
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.

//! Bit fiddling on positive IEEE 754 floats. Negative numbers aren't and needn't be handled.
//! Normal floating point numbers have a canonical representation as (frac, exp) such that the
//! value is 2^exp * (1 + sum(frac[N-i] / 2^i)) where N is the number of bits. Subnormals are
//! slightly different and weird, but the same principle applies.
//!
//! Here, however, we represent them as (sig, k) with f positive, such that the value is f * 2^e.
//! Besides making the "hidden bit" explicit, this changes the exponent by the so-called
//! mantissa shift.
//!
//! Put another way, normally floats are written as (1) but here they are written as (2):
//!
//! 1. `1.101100...11 * 2^m`
//! 2. `1101100...11 * 2^n`
//!
//! We call (1) the **fractional representation** and (2) the **integral representation**.
//!
//! Many functions in this module only handle normal numbers. The dec2flt routines conservatively
//! take the universally-correct slow path (Algorithm M) for very small and very large numbers.
//! That algorithm needs only next_float() which does handle subnormals and zeros.
use prelude::v1::*;
use u32;
use cmp::Ordering::{Less, Equal, Greater};
use ops::{Mul, Div, Neg};
use fmt::{Debug, LowerExp};
use mem::transmute;
use num::diy_float::Fp;
use num::FpCategory::{Infinite, Zero, Subnormal, Normal, Nan};
use num::Float;
use num::dec2flt::num::{self, Big};
use num::dec2flt::table;

#[derive(Copy, Clone, Debug)]
pub struct Unpacked {
    pub sig: u64,
    pub k: i16,
}

impl Unpacked {
    pub fn new(sig: u64, k: i16) -> Self {
        Unpacked { sig: sig, k: k }
    }
}

/// A helper trait to avoid duplicating basically all the conversion code for `f32` and `f64`.
///
/// See the parent module's doc comment for why this is necessary.
///
/// Should **never ever** be implemented for other types or be used outside the dec2flt module.
/// Inherits from `Float` because there is some overlap, but all the reused methods are trivial.
/// The "methods" (pseudo-constants) with default implementation should not be overriden.
pub trait RawFloat : Float + Copy + Debug + LowerExp
                    + Mul<Output=Self> + Div<Output=Self> + Neg<Output=Self>
{
    // suffix of "2" because Float::infinity is deprecated
    #[allow(deprecated)]
    fn infinity2() -> Self {
        Float::infinity()
    }

    // suffix of "2" because Float::nan is deprecated
    #[allow(deprecated)]
    fn nan2() -> Self {
        Float::nan()
    }

    // suffix of "2" because Float::zero is deprecated
    fn zero2() -> Self;

    // suffix of "2" because Float::integer_decode is deprecated
    #[allow(deprecated)]
    fn integer_decode2(self) -> (u64, i16, i8) {
        Float::integer_decode(self)
    }

    /// Get the raw binary representation of the float.
    fn transmute(self) -> u64;

    /// Transmute the raw binary representation into a float.
    fn from_bits(bits: u64) -> Self;

    /// Decode the float.
    fn unpack(self) -> Unpacked;

    /// Cast from a small integer that can be represented exactly.  Panic if the integer can't be
    /// represented, the other code in this module makes sure to never let that happen.
    fn from_int(x: u64) -> Self;

    /// Get the value 10^e from a pre-computed table. Panics for e >= ceil_log5_of_max_sig().
    fn short_fast_pow10(e: usize) -> Self;

    // FIXME Everything that follows should be associated constants, but taking the value of an
    // associated constant from a type parameter does not work (yet?)
    // A possible workaround is having a `FloatInfo` struct for all the constants, but so far
    // the methods aren't painful enough to rewrite.

    /// What the name says. It's easier to hard code than juggling intrinsics and
    /// hoping LLVM constant folds it.
    fn ceil_log5_of_max_sig() -> i16;

    // A conservative bound on the decimal digits of inputs that can't produce overflow or zero or
    /// subnormals. Probably the decimal exponent of the maximum normal value, hence the name.
    fn max_normal_digits() -> usize;

    /// When the most significant decimal digit has a place value greater than this, the number
    /// is certainly rounded to infinity.
    fn inf_cutoff() -> i64;

    /// When the most significant decimal digit has a place value less than this, the number
    /// is certainly rounded to zero.
    fn zero_cutoff() -> i64;

    /// The number of bits in the exponent.
    fn exp_bits() -> u8;

    /// The number of bits in the singificand, *including* the hidden bit.
    fn sig_bits() -> u8;

    /// The number of bits in the singificand, *excluding* the hidden bit.
    fn explicit_sig_bits() -> u8 {
        Self::sig_bits() - 1
    }

    /// The maximum legal exponent in fractional representation.
    fn max_exp() -> i16 {
        (1 << (Self::exp_bits() - 1)) - 1
    }

    /// The minimum legal exponent in fractional representation, excluding subnormals.
    fn min_exp() -> i16 {
        -Self::max_exp() + 1
    }

    /// `MAX_EXP` for integral representation, i.e., with the shift applied.
    fn max_exp_int() -> i16 {
        Self::max_exp() - (Self::sig_bits() as i16 - 1)
    }

    /// `MAX_EXP` encoded (i.e., with offset bias)
    fn max_encoded_exp() -> i16 {
        (1 << Self::exp_bits()) - 1
    }

    /// `MIN_EXP` for integral representation, i.e., with the shift applied.
    fn min_exp_int() -> i16 {
        Self::min_exp() - (Self::sig_bits() as i16 - 1)
    }

    /// The maximum normalized singificand in integral representation.
    fn max_sig() -> u64 {
        (1 << Self::sig_bits()) - 1
    }

    /// The minimal normalized significand in integral representation.
    fn min_sig() -> u64 {
        1 << (Self::sig_bits() - 1)
    }
}

impl RawFloat for f32 {
    fn zero2() -> Self {
        0.0
    }

    fn sig_bits() -> u8 {
        24
    }

    fn exp_bits() -> u8 {
        8
    }

    fn ceil_log5_of_max_sig() -> i16 {
        11
    }

    fn transmute(self) -> u64 {
        let bits: u32 = unsafe { transmute(self) };
        bits as u64
    }

    fn from_bits(bits: u64) -> f32 {
        assert!(bits < u32::MAX as u64, "f32::from_bits: too many bits");
        unsafe { transmute(bits as u32) }
    }

    fn unpack(self) -> Unpacked {
        let (sig, exp, _sig) = self.integer_decode2();
        Unpacked::new(sig, exp)
    }

    fn from_int(x: u64) -> f32 {
        // rkruppe is uncertain whether `as` rounds correctly on all platforms.
        debug_assert!(x as f32 == fp_to_float(Fp { f: x, e: 0 }));
        x as f32
    }

    fn short_fast_pow10(e: usize) -> Self {
        table::F32_SHORT_POWERS[e]
    }

    fn max_normal_digits() -> usize {
        35
    }

    fn inf_cutoff() -> i64 {
        40
    }

    fn zero_cutoff() -> i64 {
        -48
    }
}


impl RawFloat for f64 {
    fn zero2() -> Self {
        0.0
    }

    fn sig_bits() -> u8 {
        53
    }

    fn exp_bits() -> u8 {
        11
    }

    fn ceil_log5_of_max_sig() -> i16 {
        23
    }

    fn transmute(self) -> u64 {
        let bits: u64 = unsafe { transmute(self) };
        bits
    }

    fn from_bits(bits: u64) -> f64 {
        unsafe { transmute(bits) }
    }

    fn unpack(self) -> Unpacked {
        let (sig, exp, _sig) = self.integer_decode2();
        Unpacked::new(sig, exp)
    }

    fn from_int(x: u64) -> f64 {
        // rkruppe is uncertain whether `as` rounds correctly on all platforms.
        debug_assert!(x as f64 == fp_to_float(Fp { f: x, e: 0 }));
        x as f64
    }

    fn short_fast_pow10(e: usize) -> Self {
        table::F64_SHORT_POWERS[e]
    }

    fn max_normal_digits() -> usize {
        305
    }

    fn inf_cutoff() -> i64 {
        310
    }

    fn zero_cutoff() -> i64 {
        -326
    }

}

/// Convert an Fp to the closest f64. Only handles number that fit into a normalized f64.
pub fn fp_to_float<T: RawFloat>(x: Fp) -> T {
    let x = x.normalize();
    // x.f is 64 bit, so x.e has a mantissa shift of 63
    let e = x.e + 63;
    if e > T::max_exp() {
        panic!("fp_to_float: exponent {} too large", e)
    }  else if e > T::min_exp() {
        encode_normal(round_normal::<T>(x))
    } else {
        panic!("fp_to_float: exponent {} too small", e)
    }
}

/// Round the 64-bit significand to 53 bit with half-to-even. Does not handle exponent overflow.
pub fn round_normal<T: RawFloat>(x: Fp) -> Unpacked {
    let excess = 64 - T::sig_bits() as i16;
    let half: u64 = 1 << (excess - 1);
    let (q, rem) = (x.f >> excess, x.f & ((1 << excess) - 1));
    assert_eq!(q << excess | rem, x.f);
    // Adjust mantissa shift
    let k = x.e + excess;
    if rem < half {
        Unpacked::new(q, k)
    } else if rem == half && (q % 2) == 0 {
        Unpacked::new(q, k)
    } else if q == T::max_sig() {
        Unpacked::new(T::min_sig(), k + 1)
    } else {
        Unpacked::new(q + 1, k)
    }
}

/// Inverse of `RawFloat::unpack()` for normalized numbers.
/// Panics if the significand or exponent are not valid for normalized numbers.
pub fn encode_normal<T: RawFloat>(x: Unpacked) -> T {
    debug_assert!(T::min_sig() <= x.sig && x.sig <= T::max_sig(),
        "encode_normal: significand not normalized");
    // Remove the hidden bit
    let sig_enc = x.sig & !(1 << T::explicit_sig_bits());
    // Adjust the exponent for exponent bias and mantissa shift
    let k_enc = x.k + T::max_exp() + T::explicit_sig_bits() as i16;
    debug_assert!(k_enc != 0 && k_enc < T::max_encoded_exp(),
        "encode_normal: exponent out of range");
    // Leave sign bit at 0 ("+"), our numbers are all positive
    let bits = (k_enc as u64) << T::explicit_sig_bits() | sig_enc;
    T::from_bits(bits)
}

/// Construct the subnormal. A mantissa of 0 is allowed and constructs zero.
pub fn encode_subnormal<T: RawFloat>(significand: u64) -> T {
    assert!(significand < T::min_sig(), "encode_subnormal: not actually subnormal");
    // Encoded exponent is 0, the sign bit is 0, so we just have to reinterpret the bits.
    T::from_bits(significand)
}

/// Approximate a bignum with an Fp. Rounds within 0.5 ULP with half-to-even.
pub fn big_to_fp(f: &Big) -> Fp {
    let end = f.bit_length();
    assert!(end != 0, "big_to_fp: unexpectedly, input is zero");
    let start = end.saturating_sub(64);
    let leading = num::get_bits(f, start, end);
    // We cut off all bits prior to the index `start`, i.e., we effectively right-shift by
    // an amount of `start`, so this is also the exponent we need.
    let e = start as i16;
    let rounded_down = Fp { f: leading, e: e }.normalize();
    // Round (half-to-even) depending on the truncated bits.
    match num::compare_with_half_ulp(f, start) {
        Less => rounded_down,
        Equal if leading % 2 == 0 => rounded_down,
        Equal | Greater => match leading.checked_add(1) {
            Some(f) => Fp { f: f, e: e }.normalize(),
            None => Fp { f: 1 << 63, e: e + 1 },
        }
    }
}

/// Find the largest floating point number strictly smaller than the argument.
/// Does not handle subnormals, zero, or exponent underflow.
pub fn prev_float<T: RawFloat>(x: T) -> T {
    match x.classify() {
        Infinite => panic!("prev_float: argument is infinite"),
        Nan => panic!("prev_float: argument is NaN"),
        Subnormal => panic!("prev_float: argument is subnormal"),
        Zero => panic!("prev_float: argument is zero"),
        Normal => {
            let Unpacked { sig, k } = x.unpack();
            if sig == T::min_sig() {
                encode_normal(Unpacked::new(T::max_sig(), k - 1))
            } else {
                encode_normal(Unpacked::new(sig - 1, k))
            }
        }
    }
}

// Find the smallest floating point number strictly larger than the argument.
// This operation is saturating, i.e. next_float(inf) == inf.
// Unlike most code in this module, this function does handle zero, subnormals, and infinities.
// However, like all other code here, it does not deal with NaN and negative numbers.
pub fn next_float<T: RawFloat>(x: T) -> T {
    match x.classify() {
        Nan => panic!("next_float: argument is NaN"),
        Infinite => T::infinity2(),
        // This seems too good to be true, but it works.
        // 0.0 is encoded as the all-zero word. Subnormals are 0x000m...m where m is the mantissa.
        // In particular, the smallest subnormal is 0x0...01 and the largest is 0x000F...F.
        // The smallest normal number is 0x0010...0, so this corner case works as well.
        // If the increment overflows the mantissa, the carry bit increments the exponent as we
        // want, and the mantissa bits become zero. Because of the hidden bit convention, this
        // too is exactly what we want!
        // Finally, f64::MAX + 1 = 7eff...f + 1 = 7ff0...0 = f64::INFINITY.
        Zero | Subnormal | Normal => {
            let bits: u64 = x.transmute();
            T::from_bits(bits + 1)
        }
    }
}