R 中数值的双精度（64 位）表示（符号、指数、有效数）答案

【问题标题】：Double precision (64-bit) representation of numeric value in R (sign, exponent, significand)R 中数值的双精度（64 位）表示（符号、指数、有效数）
【发布时间】：2018-05-07 15:44:10
【问题描述】：

R FAQ 声明：

唯一可以在 R 的数字类型中精确表示的数字是整数和分母是 2 的幂的分数。所有其他数字在内部四舍五入到（通常）53 位二进制精度。

R 使用 IEEE 754 双精度浮点数

1 位符号
11 位的指数
52 位用于尾数（或有效位）

总计为 64 位。

对于数字0.1，R代表

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

Double (IEEE754 Double precision 64-bit) 为我们提供了0.1 的二进制表示：

00111111 10111001 10011001 10011001 10011001 10011001 10011001 10011010

我们如何在 R 中获得这种表示形式，它与我们示例中 sprintf 给出的输出有何关系？

【问题讨论】：

R 是否支持sprintf("%a", 0.1)？

标签： r floating-point numeric

【解决方案1】：

@chux 在 cmets 中提出的问题的答案是“是”； R 支持%a 格式：

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

如果要访问底层位模式，则必须将双精度重新解释为 64 位整数。对于这项任务，可以通过 Rcpp 使用 C++：

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

此十六进制表示与您的二进制表示相同。如何获得十进制表示？

第一位为零，因此符号为正
指数为 0x3fb，即十进制的 1019。鉴于 exponent bias，这对应于 -4 的实际指数。
尾数为 0x1999999999999a × 2^-52，包括 implicit 1，即 2^−52 × 7,205,759,403,792,794。

总共得到 2^−56 × 7,205,759,403,792,794：

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

【讨论】：

【解决方案2】：

以0.3 为例。在 R 控制台中运行

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

尾数或有效位

十六进制表示 3333333333333 到二进制会给我们尾数（或有效数）部分。那就是

0011001100110011001100110011001100110011001100110011

指数

指数部分（11 位）应该是从 2^(11-1) - 1 = 1023 的偏移量，因此尾随 3 是 p-2（在 sprintf 给出的输出中）我们有

-2 + 1023 = 1021

其二进制表示固定为11位

01111111101

签名

符号位为0，否则为1

双精度表示

所以完整的表示是

0 | 01111111101 | 0011001100110011001100110011001100110011001100110011

另一个例子：

> 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

【讨论】：

如何从base16传递到base2？
base16 中的每个数字都相当于该数字的 4 位二进制表示。例如，3 ==> 0011

【解决方案3】：

从十进制到标准化双精度：

library(BMS)

from10toNdp <- function(my10baseNumber) {
out <- list()

# Handle special cases (0, Inf, -Inf)
if (my10baseNumber %in% c(0,Inf,-Inf)) {
if (my10baseNumber==0)    { out <- "0000000000000000000000000000000000000000000000000000000000000000" }
if (my10baseNumber==Inf)  { out <- "0111111111110000000000000000000000000000000000000000000000000000" }
if (my10baseNumber==-Inf) { out <- "1111111111110000000000000000000000000000000000000000000000000000" }
} else {

signBit <- 0 # assign initial value

from10to2 <- function(deciNumber) {
  binaryVector <- rep(0, 1 + floor(log(deciNumber, 2)))
  while (deciNumber >= 2) {
    theExpo <- floor(log(deciNumber, 2))
    binaryVector[1 + theExpo] <- 1
    deciNumber <- deciNumber - 2^theExpo  }
  binaryVector[1] <- deciNumber %% 2
  paste(rev(binaryVector), collapse = "")}

#Sign bit
if (my10baseNumber<0) { signBit <- 1 
} else { signBit <- 0 }

# Biased Exponent
BiasedExponent <- strsplit(from10to2(as.numeric(substr(sprintf("%a", my10baseNumber), which(strsplit( sprintf("%a", my10baseNumber), "")[[1]]=="p")+1, length( strsplit( sprintf("%a", my10baseNumber), "")[[1]]))) + 1023), "")[[1]] 
BiasedExponent <- paste(BiasedExponent, collapse='')
if (nchar(BiasedExponent)<11) {BiasedExponent <-  paste(c(  rep(0,11-nchar(BiasedExponent)), BiasedExponent),collapse='')    }

# Significand
significand <- BMS::hex2bin(substr( sprintf("%a", my10baseNumber) , which(strsplit( sprintf("%a", my10baseNumber), "")[[1]]=="x")+3, which(strsplit( sprintf("%a", my10baseNumber), "")[[1]]=="p")-1))

significand <- paste(significand, collapse='')
if (nchar(significand)<52) {significand <-  paste(c( significand,rep(0,52-nchar(significand))),collapse='')    }

out <- paste(c(signBit, BiasedExponent, significand), collapse='')
}

out
}

因此，

from10toNdp(0.1)
# "0011111110111001100110011001100110011001100110011001100110011010"

【讨论】：