N 选择 K 函数崩溃 Rcpp

Question

我在 C++ 中编写了一个“n 选择 k”函数，它通过 Rcpp 与 R 交互。出于某种原因，我收到“除以零”运行时错误。当我尝试评估 30 选择 2 时会发生这种情况。

我尝试过手动评估每一行（使用 evalCpp），但我仍然对除以零的位置感到困惑。也许有人可以向我指出这一点，或者建议一种更好的写 n 选择 K 的方法？

代码如下：

// [[Rcpp::export]]                                                                                                                                  
int chooseC(int n, int k) {                                                                                                                         
  if (k > n) {                                                                                                                                      
    std::cout << "Error. k cannot be greater than n." << std::endl;                                                                                 
    return 0;                                                                                                                                       
  }                                                                                                                                                 
  int factN = std::tgamma(n + 1);                                                                                                                   
  int factK = std::tgamma(k + 1);                                                                                                                   
  int factDiff = std::tgamma(n - k + 1);                                                                                                            
  return factN/(factK*factDiff);                                                                                                                    
} 

Answer 1

简单地说：

• 据我所知，std 中没有 tgamma

• R 本身是一个choose 函数，所以我会做下面的事情

• R 还具有伽马分布等，因此您也可以手动执行此操作

• 为什么不直接打印值 factN、factK、factDiff？

简单的 Rcpp 解决方案：

#include <Rcpp.h>

// [[Rcpp::export]]  
double chooseC(double n, double k) {
  return Rf_choose(n, k);
}

例子：

R> chooseC(5,2)     
[1] 10
R> 

编辑：根据@Blastfurnace 在C++11 cmath 标头中对tgamma() 的评论，这是一个修复后的版本，对我来说很好：

#include <Rcpp.h>
#include <cmath>

// [[Rcpp::plugins(cpp11)]]

// [[Rcpp::export]] 
int chooseCtake2(int n, int k) {
  if (k > n) {
    Rcpp::stop("Error. k cannot be greater than n.");
  }
  int factN = std::tgamma(n + 1);
  int factK = std::tgamma(k + 1);
  int factDiff = std::tgamma(n - k + 1);
  return factN/(factK*factDiff); 
}

使用示例：

R> sourceCpp("/tmp/chooseC.cpp")
R> chooseCtake2(2,3)
Error: Error. k cannot be greater than n.
R> chooseCtake2(5,2)
[1] 10
R> 

Answer 2

所以std::tgamma(x) 计算 x 的 gamma 函数。这个函数很快地趋于无穷：

http://www.wolframalpha.com/share/clip?f=d41d8cd98f00b204e9800998ecf8427et5pmak8jtn

已经在 x == 31，你有一个非常大的数字。

将这个非常大的 double 转换回 int 时，结果是未定义的行为（4.9 浮点整数转换 [conv.fpint]）：

浮点类型的纯右值可以转换为 整数类型。转换截断；也就是小数部分 被丢弃。如果截断的值不能，则行为未定义 以目标类型表示。

在我的系统上，这种转换（输入为 {30, 2}）会生成一个值为 -2147483648 的 int。通过插入一些打印语句很容易观察到这一点：

int
chooseC(int n, int k)
{
    if (k > n)
    {                                                                                                                                      
        std::cout << "Error. k cannot be greater than n.\n";
        return 0;                                                                                                                                       
    }                                                                                                                                                 
    int factN = std::tgamma(n + 1);
    std::cout << "factN = " << factN << '\n';
    int factK = std::tgamma(k + 1);
    std::cout << "factK = " << factK << '\n';
    int factDiff = std::tgamma(n - k + 1);
    std::cout << "factDiff = " << factDiff << '\n';
    std::cout << "factK*factDiff = " << factK*factDiff << '\n';
    return factN/(factK*factDiff); 
}

对我来说输出：

factN = -2147483648
factK = 2
factDiff = -2147483648
factK*factDiff = 0

可以看出，UB 最终导致除以零，这也是 UB。并且听起来与您看到的行为非常相似。

这个问题的解决方案是只使用整数运算来计算事物，并且如果最终结果可以用整数类型表示，则中间计算不会溢出。这需要使用最大公约数函数。

执行此操作的开源代码可在此处获得：

http://howardhinnant.github.io/combinations.html

搜索“count_each_combination”。您的chooseC 可以用count_each_combination 编码，如下所示：

int
chooseC(int n, int k)
{
    if (k > n)
    {                                                                                                                                      
        std::cout << "Error. k cannot be greater than n.\n";
        return 0;                                                                                                                                       
    }                                                                                                                                                 
    return count_each_combination(n-k, k);
}

现在chooseC(30, 2) 将返回 435。如果count_each_combination 无法将结果存储在int 中，则会抛出std::overflow_error。

如果您想将chooseC 限制为k == 2，或者只是为了更好地理解算法而暂时这样做，请注意计算组合的公式为：

当k == 2 时，这简化为：

n*(n-1)/2

现在n 是偶数，或者n-1 是偶数。您可以发现哪个，然后将该数字除以 2，没有截断错误，然后将结果乘以未除以 2 的数字。因此，您得到准确的结果，没有截断错误的可能性，也没有中间溢出，仅使用积分算术。这是count_each_combination 使用的技术，但可以推广到任何除数，以提供始终准确的结果（如果它可以适合提供的整数类型）。