c++ - base-R seq 的 Rcpp 版本丢弃值

Question

我写了一个 Rcpp 版本的 base-R seq 函数。

library(Rcpp)

cppFunction('NumericVector seqC(double x, double y, double by) {

  // length of result vector
  int nRatio = (y - x) / by;
  NumericVector anOut(nRatio + 1);

  // compute sequence
  int n = 0;
  for (double i = x; i <= y; i = i + by) {
    anOut[n] = i;
    n += 1;
  }

  return anOut;
}')

对于以下测试，它工作得很好。

seqC(1, 11, 2)
[1]  1  3  5  7  9 11

seqC(1, 10, 2)
[1]  1  3  5  7  9 11

此外，它(有时)在传递带有十进制数字的值而不是整数。

seqC(0.43, 0.45, 0.001)
[1] 0.430 0.431 0.432 0.433 0.434 0.435 0.436 0.437 0.438 0.439 0.440 0.441 0.442 0.443 0.444 0.445 0.446 0.447 0.448 0.449 0.450

但是，自上次以来，有时该功能似乎没有按预期工作序列的条目被删除(或者更确切地说，输出 vector anOut没有合适的大小)，这 - 根据我相当缺乏的 C++ 技能，可能归因于某种舍入错误。

seqC(0.53, 0.59, 0.001)
 [1] 0.530 0.531 0.532 0.533 0.534 0.535 0.536 0.537 0.538 0.539 0.540 0.541 0.542 0.543 0.544 0.545 0.546 0.547 0.548 0.549 0.550 0.551
[23] 0.552 0.553 0.554 0.555 0.556 0.557 0.558 0.559 0.560 0.561 0.562 0.563 0.564 0.565 0.566 0.567 0.568 0.569 0.570 0.571 0.572 0.573
[45] 0.574 0.575 0.576 0.577 0.578 0.579 0.580 0.581 0.582 0.583 0.584 0.585 0.586 0.587 0.588 0.589

例如，在最后一个示例中，缺少最后一个值 (0.590)。做有人知道如何解决这个问题吗？

Answer 1

正如其他人所指出的，您遇到的问题基本上是浮点运算错误。一个常见的解决方法是将 double 缩放到足够大的整数，执行任务，然后将结果调整为输入的原始比例。我采用了与@RHertel 略有不同的方法，让缩放量 (adjust) 由增量的精度决定，而不是使用固定量，但想法本质上是相同的。

---

#include <Rcpp.h>

struct add_multiple {
  int incr;
  int count;
  add_multiple(int incr)
    : incr(incr), count(0)
    {}
  inline int operator()(int d) {
    return d + incr * count++;
  }
};

// [[Rcpp::export]]
Rcpp::NumericVector rcpp_seq(double from_, double to_, double by_ = 1.0) {
  int adjust = std::pow(10, std::ceil(std::log10(10 / by_)) - 1);
  int from = adjust * from_;
  int to = adjust * to_;
  int by = adjust * by_;

  std::size_t n = ((to - from) / by) + 1;
  Rcpp::IntegerVector res = Rcpp::rep(from, n);
  add_multiple ftor(by);

  std::transform(res.begin(), res.end(), res.begin(), ftor);
  return Rcpp::NumericVector(res) / adjust;
}

/*** R

all.equal(seq(.53, .59, .001), seqC(.53, .59, .001)) &&
  all.equal(seq(.53, .59, .001), rcpp_seq(.53, .59, .001))
# [1] TRUE

all.equal(seq(.53, .54, .000001), seqC(.53, .54, .000001)) &&
  all.equal(seq(.53, .54, .000001), rcpp_seq(.53, .54, .000001))
# [1] TRUE 

microbenchmark::microbenchmark(
  "seq" = seq(.53, .54, .000001),
  "seqC" = seqC(0.53, 0.54, 0.000001),
  "rcpp_seq" = rcpp_seq(0.53, 0.54, 0.000001),
  times = 100L)
Unit: microseconds
     expr        min          lq        mean     median         uq        max neval
      seq    896.190   1015.7940   1167.4708   1132.466   1221.624   1651.571   100
     seqC 212293.307 219527.6590 226933.4329 223384.592 227860.410 398462.561   100
 rcpp_seq    182.848    194.1665    225.4338    227.396    244.942    320.436   100

*/ 

---

seqC 是@RHertel 修改后的实现，它产生了正确的结果。 FWIW 我认为此功能的缓慢性能主要是由于在 NumericVector 类型上使用 push_back 造成的，Rcpp 开发人员强烈建议不要这样做。