使用gmp和机器极限处理大整数

Question

我想知道在R中是否可以使用比.Machine$double.xmax(~1.79e308)更大的整数。 我认为通过在R中使用例如Rmpfr或gmp库，您可以分配任何大小的值，直到您系统上RAM的限制。 我认为这比.Machine$double.xmax大，但显然不是这样。

> require( gmp )
> as.bigz( .Machine$double.xmax )
Big Integer ('bigz') :
[1] 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368
> as.bigz( 1e309 )
Big Integer ('bigz') :
[1] NA
> 

有没有人能解释一下为什么使用64位内存寻址的计算机无法存储大于1.79e308的值？抱歉，我没有计算机科学背景，但我正在努力学习。

谢谢。

- Simon O'Hanlon

---

你可能会发现阅读文档很有帮助，特别是 ?as.bigz 的注释部分。 - joran

谢谢Joran。我错过了最后一行。令人恼火的是，现在我不能使用科学计数法！ - Simon O'Hanlon

真的，但你可以只需执行 as.bigz(10)^309。事实上，你可以这样做："%e%" <- function(x,y) as.bigz(x) * 10^as.bigz(y); 1%e%309 - Ben Bolker

PS：我的巧妙黑科技只适用于整数 x，所以如果你想要 1.5e309，你需要类似 15%e%308 的东西。 - Ben Bolker

请注意，http://r-bc.googlecode.com 上的 bc R 包也不直接支持科学计数法，但可以处理：e <- bc(10); 1.5 * e ^ 309。 - G. Grothendieck

@BenBolker @g-grothendieck 感谢您的有益评论和不错的技巧！我之前没有听说过 bc，所以我会去了解一下。谢谢。 - Simon O'Hanlon

Answer 1

Rmpfr可以使用mpfr_set_str进行字符串转换...

Rmpfr是一种用于高精度计算的R语言包。它能够使用mpfr_set_str函数来进行字符串转换。

val <- mpfr("1e309")

## 1 'mpfr' number of precision  17   bits 
## [1] 9.999997e308

# set a precision (assume base 10)...
est_prec <- function(e) floor( e/log10(2) ) + 1

val <- mpfr("1e309", est_prec(309) )

## 1 'mpfr' number of precision  1027   bits 
## [1]1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

.mpfr2bigz(val)

## Big Integer ('bigz') :
## [1] 1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

# extract exponent from a scientific notation string
get_exp <- function( sci ) as.numeric( gsub("^.*e",'', sci) )

# Put it together
sci2bigz <- function( str ) {
  .mpfr2bigz( mpfr( str, est_prec( get_exp( str ) ) ) )
}

val <- sci2bigz( paste0( format( Const("pi", 1027) ), "e309") )

identical( val, .mpfr2bigz( Const("pi",1027)*mpfr(10,1027)^309 ) )

## [1] TRUE

## Big Integer ('bigz') :
## [1] 3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587004

关于为什么要存储比 .Machine$double.xmax 更大的数字，IEEE规范中有关浮点数编码的文档、R常见问题解答和维基百科都涉及了所有术语，但我发现仅定义这些术语（使用'?'.Machine'）会很有帮助...

double.xmax（最大的规格化浮点数）=
(1-double.neg.eps) * double.base ^ double.max.exp，其中

1. double.neg.eps（一个较小的正浮点数x，使得1-x！= 1）= double.base ^ double.neg.ulp.digits，其中
   • double.neg.ulp.digits = 最大的负整数，使得1-double.base ^ i！= 1
2. double.max.exp = 溢出的最小正 double.base 幂
3. double.base（浮点表示法的基数）= 2（二进制）。

考虑可以从另一个区分有限浮点数的区间；IEEE规范告诉我们，在二进制64位数中有11位用于指数，因此我们有最大指数2^(11-1)-1=1023但我们想要最大溢出的指数，所以double.max.exp是1024。

# Maximum number of representations
# double.base ^ double.max.exp
base <- mpfr(2, 2048)
max.exp <- mpfr( 1024, 2048 )

# This is where the big part of the 1.79... comes from
base^max.exp

## 1 'mpfr' number of precision  2048   bits 
## [1] 179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137216

# Smallest definitive unit.
# Find the largest negative integer...
neg.ulp.digits <- -64; while( ( 1 - 2^neg.ulp.digits ) == 1 ) 
  neg.ulp.digits <<- neg.ulp.digits + 1

neg.ulp.digits

## [1] -53

# It makes a real small number...
neg.eps <- base^neg.ulp.digits

neg.eps

## 1 'mpfr' number of precision  2048   bits 
## [1] 1.11022302462515654042363166809082031250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-16

# Largest difinitive floating point number less than 1
# times the number of representations
xmax <- (1-neg.eps) * base^max.exp

xmax

## 1 'mpfr' number of precision  2048   bits 
## [1] 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368

identical( asNumeric(xmax), .Machine$double.xmax )

## [1] TRUE

- Thell