nextafter功能在R中的实现
nextafter功能在R中的实现
提问于 2014-02-27 00:24:11
在R中有没有功能的实现，这样就可以从给定的浮点数中获得下一个可表示的浮点数。这将类似于C标准库中的nextafter函数。像number + .Machine$double.eps这样的方案通常不起作用。
回答 3
Stack Overflow用户
发布于 2014-02-27 00:46:39
不能，但有两种方法可以做到：
使用C的
- 该函数不返回值。所有的工作都是作为“副作用”来完成的(改变arguments).
- All的值，参数就是指针。偶数标量是R.
中的向量(长度为1)
然后将该函数编译为共享库：
R CMD SHLIB foo.c用于类UNIX的OSs。可以使用dyn.load("foo.so")调用共享库。然后，可以使用.C()函数从R内部调用该函数
.C("foo", ...)从R调用C的一个更深入的处理是here。
使用R的
number + .Machine$double.eps是可行的，但您必须考虑边缘情况，例如if x - y < .Machine$double.eps或if x == y。我会像这样写这个函数：
nextafter <- function(x, y){
# Appropriate type checking and bounds checking goes here
delta = y - x
if(x > 0){
factor = 2^floor(log2(x)) + ifelse(x >= 4, 1, 0)
} else if (x < 0) {
factor = 65
}
if (delta > .Machine$double.eps){
return(x + factor * .Machine$double.eps)
} else if (delta < .Machine$double.eps){
return(x - factor * .Machine$double.eps)
} else {
return(x)
}
}UPDATE对于大于2的数字，前面的代码未按预期执行。有一个因子需要与.Machine$double.eps相乘，才能使其足够大，从而导致数字不同。它与2加1的最接近的幂有关。你可以通过下面的代码来了解它是如何工作的：
n <- -100
factor <- vector('numeric', 100)
for(i in 1:n){
j = 0
while(TRUE){
j = j + 1
if(i - j * .Machine$double.eps != i) break()
}
factor[i] = j
} Stack Overflow用户
发布于 2014-02-27 01:01:08
如果您更喜欢Rcpp：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double nextAfter(double x, double y) {
return nextafter(x, y);
}然后在R中：
sprintf("%.20f", 1)
#[1] "1.00000000000000000000"
sprintf("%.20f", nextAfter(1, 2))
#[1] "1.00000000000000022204"Stack Overflow用户
发布于 2021-11-03 02:15:12
我不确定Christopher Louden的答案是否适用于所有的值，但这里有一个经典方法的纯R版本(递增/递减整数位)。R不能方便地在doubles和integer之间进行转换，也没有64位的整数类型，因此需要编写大量代码。
doubleToRaw <- function(d) writeBin(d, raw());
rawToDouble <- function(r) readBin(r, numeric());
int64inc <- function(lo, hi) {
if (lo == 0xffffffff) { hi <- hi + 1; lo <- 0; } else { lo <- lo + 1; }
return(c(lo, hi));
}
int64dec <- function(lo, hi) {
if (lo == 0) { hi <- hi - 1; lo <- 0xffffffff; } else { lo <- lo - 1; }
return(c(lo, hi));
}
nextafter <- function(x, y) {
if (is.nan(x + y))
return(NaN);
if (x == y)
return(x);
if (x == 0)
return(sign(y) * rawToDouble(as.raw(c(0, 0, 0, 0, 0, 0, 0, 1))));
ints <- packBits(rawToBits(doubleToRaw(x)), "integer")
if ((y > x) == (x > 0))
ints <- int64inc(ints[1], ints[2])
else
ints <- int64dec(ints[1], ints[2]);
return(rawToDouble(packBits(intToBits(ints), "raw")))
}页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/22047238
复制相关文章
相似问题
腾讯云开发者
