Double precision (64-bit) representation of numeri

2019-03-31 19:09发布

问题:

R FAQ states that:

The only numbers that can be represented exactly in R’s numeric type are integers and fractions whose denominator is a power of 2. All other numbers are internally rounded to (typically) 53 binary digits accuracy.

R uses IEEE 754 double-precision floating-point numbers which is

  • 1 bit for sign
  • 11 bits for exponent
  • 52 bits for mantissa (or significand)

which sums up to 64-bits.

For the numeric number 0.1, R represents

sprintf("%.60f", 0.1)
[1] "0.100000000000000005551115123125782702118158340454101562500000"

Double (IEEE754 Double precision 64-bit) gives us this binary representation for 0.1 :

00111111 10111001 10011001 10011001 10011001 10011001 10011001 10011010

How we can get this representation in R and how does it relate to the output given by sprintf in our example?

回答1:

The answer to the question raised by @chux in the comments is "yes"; R supports the %a format:

sprintf("%a", 0.1)
#> [1] "0x1.999999999999ap-4"

If you want to access the underlying bit pattern, you will have to reinterpret the double as a 64bit integer. For this task one can use C++ via Rcpp:

Rcpp::cppFunction('void print_hex(double x) {
    uint64_t y;
    static_assert(sizeof x == sizeof y, "Size does not match!");
    std::memcpy(&y, &x, sizeof y);
    Rcpp::Rcout << std::hex << y << std::endl;
}', plugins = "cpp11", includes = "#include <cstdint>")
print_hex(0.1)
#> 3fb999999999999a

This hexadecimal representation is identical to your binary representation. How does one get to the decimal representation?

  • The first bit is zero, hence the sign is positive
  • The exponent is 0x3fb, i.e. 1019 in decimal. Given the exponent bias this corresponds to an actual exponent of -4.
  • The mantissa is 0x1999999999999a × 2^-52 including the implicit 1, i.e. 2^−52 × 7,205,759,403,792,794.
  • In total this gives 2^−56 × 7,205,759,403,792,794:

    sprintf("%.60f", 2^-56 * 7205759403792794)
    #> [1] "0.100000000000000005551115123125782702118158340454101562500000"
    


回答2:

Take for example 0.3 into account. Run in R console

> sprintf("%a", 0.3)
[1] "0x1.3333333333333p-2"

Mantissa or Significand

The hex representation 3333333333333 to binary would give us the mantissa (or significand) part. That is

0011001100110011001100110011001100110011001100110011

Exponent

The exponent part (11 bits) should be the offset from 2^(11-1) - 1 = 1023 so as the trailing 3 is p-2 (in the output given by sprintf) we have

-2 + 1023 = 1021

and its binary representation fixed in 11 bits is

01111111101

Sign

As for the sign bit, its 0 for positive and 1 otherwise

Double Precision Representation

So the complete representation is

0 | 01111111101 | 0011001100110011001100110011001100110011001100110011

Another example:

> sprintf("%a", -2.94)
[1] "-0x1.7851eb851eb85p+1"

# Mantissa or Significand
(7851eb851eb85) # base 16 
(0111100001010001111010111000010100011110101110000101) # base 2

# Exponent
1 + 1023 = 1024 # base 10
10000000000 # base 2

# So the complete representation is
1 | 10000000000 | 0111100001010001111010111000010100011110101110000101