La función N Elija K bloquea Rcpp

Question

Escribí una función 'n elegir k' en C++, que interactúa con R a través de Rcpp.Por alguna razón, aparece un error de tiempo de ejecución de "división por cero".Sucede cuando intento evaluar 30, elige 2.

Intenté evaluar cada línea manualmente (con evalCpp) y todavía estoy desconcertado acerca de dónde ocurre la división por cero.¿Quizás alguien podría señalarme esto o sugerirme una mejor forma de escribir y elegir K?

Aquí está el código:

// [[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

Brevemente:

• Por lo que puedo ver, no hay tgamma en std.

• R en sí mismo como choose funciona así que simplemente haría lo que está debajo

• R también tiene distribución gamma, etc., por lo que también puedes hacerlo a mano.

• ¿Por qué no imprimiste los valores? factN, factK, factDiff ?

Solución Rcpp simple:

#include <Rcpp.h>

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

Ejemplo:

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

Editar: Siguiendo el comentario de @Blastfurnace sobre tgamma() en el C ++ 11 cmath encabezado, aquí hay una versión reparada que funciona bien para mí:

#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); 
}

Uso de ejemplo:

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

Entonces std::tgamma(x) calcula la función gamma de x.Esta función llega al infinito con bastante rapidez:

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

Ya en x == 31, tienes un número muy grande.

Al convertir este doble muy grande a int, los resultados son un comportamiento indefinido (4.9 Conversiones integrales flotantes [conv.fpint]):

Se puede convertir un pralentador de un tipo de punto flotante a un pralente de un tipo entero.La conversión se trunca;es decir, la parte fraccional se descarta.El comportamiento no está definido si el valor truncado no puede representarse en el tipo de destino.

En mi sistema, esta conversión (con una entrada de {30, 2}) da como resultado un int con el valor -2147483648.Esto se observa fácilmente insertando algunas declaraciones impresas:

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); 
}

que para mí produce:

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

Como puede verse, la UB resulta finalmente en una división por cero, que también es UB.Y suena muy similar al comportamiento que estás viendo.

La solución a este problema es calcular cosas usando sólo aritmética integral, y de tal manera que los cálculos intermedios no se desborden si el resultado final es representable en el tipo integral.Esto implica el uso de una función del máximo común divisor.

El código fuente abierto que hace esto está disponible aquí:

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

Busque "count_each_combination".Su chooseC se puede codificar en términos de count_each_combination al igual que:

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);
}

Ahora chooseC(30, 2) devolverá 435.Si count_each_combination no puede almacenar el resultado en un int, a std::overflow_error será arrojado.

Si desea limitar su chooseC a k == 2, o quizás hacerlo solo temporalmente para comprender mejor el algoritmo, tenga en cuenta que la fórmula para contar combinaciones es:

Cuando k == 2, esto se simplifica a:

n*(n-1)/2

Ahora tampoco n es par, o n-1 incluso.Puedes descubrir cuál y luego dividir ese número por 2, sin error de truncamiento, y luego multiplicar el resultado por el número que no se dividió por 2.De esta forma se obtiene el resultado exacto sin posibilidad de error de truncamiento, ni desbordamiento intermedio, utilizando únicamente aritmética integral.Esta es la técnica que count_each_combination utiliza, pero generalizado a cualquier divisor, para entregar un resultado que siempre es exacto si puede encajar en el tipo integral proporcionado.