N Выберите K. Функция вызывает сбой Rcpp
-
20-12-2019 - |
Вопрос
Я написал функцию «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);
}
Решение
Кратко:
Насколько я могу судить, в 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>
Другие советы
Так 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
использует, но обобщает на любой делитель, чтобы получить результат, который всегда является точным, если он может вписаться в предоставленный целочисленный тип.