N Выберите K. Функция вызывает сбой Rcpp

Question

Я написал функцию «n select k» на C++, которая взаимодействует с R через Rcpp.По какой-то причине я получаю ошибку времени выполнения «деление на ноль».Это происходит, когда я пытаюсь оценить 30, выбираю 2.

Я пробовал оценивать каждую строку вручную (с помощью evalCpp), но до сих пор не понимаю, где происходит деление на ноль.Может быть, кто-нибудь мог бы указать мне на это или предложить лучший способ написания и выбрать 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 о tgamma() в C++ 11 cmath заголовок, вот исправленная версия, которая отлично работает для меня:

#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.Эта функция довольно быстро стремится к бесконечности:

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

Уже при x == 31 у вас очень большое число.

При преобразовании этого очень большого двойного значения обратно в int результаты имеют неопределенное поведение (4.9 Преобразования с плавающим интегралом [conv.fpint]):

Prvalue типа с плавающей запятой может быть преобразована в Prvalue целочисленного типа.Преобразование усекается;то есть дробная часть отбрасывается.Поведение не определено, если усеченное значение не может быть представлено в типе назначения.

В моей системе это преобразование (с вводом {30, 2}) приводит к int со значением -2147483648.Это легко увидеть, вставив несколько операторов печати:

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 использует, но обобщает на любой делитель, чтобы получить результат, который всегда является точным, если он может вписаться в предоставленный целочисленный тип.