N Choisir La Fonction K Se Bloque Rcpp

Question

J'ai écrit un 'n choisir k' en fonction en C++, qui est à l'interface avec la R via Rcpp.Pour une raison que je suis une "division par zéro" erreur d'exécution.Il se passe lorsque j'essaie d'évaluer les 30 choisissez 2.

J'ai essayé de l'évaluation de chaque ligne manuellement (avec evalCpp), et je suis toujours perplexe sur l'endroit où la division par zéro qui se passe.Peut-être que quelqu'un pourrait me le ou suggérer une meilleure façon d'écrire n choisir K?

Voici le code:

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

Brièvement:

• Il n'y a pas de tgamma en std aussi loin que je peux voir

• R lui-même comme un choose fonction donc, je voudrais juste faire ce qui est ci-dessous

• R a également la distribution gamma, etc afin que vous pouvez le faire à la main ainsi

• Pourquoi ne pas vous suffit d'imprimer les valeurs factN, factK, factDiff ?

Simple Rcpp solution:

#include <Rcpp.h>

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

Exemple:

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

Edit: Suivant le commentaire de @Blastfurnace sur tgamma() dans le C++11 cmath l'en-tête, voici une version réparée qui fonctionne très bien pour moi:

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

Exemple d'utilisation:

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

Donc std::tgamma(x) calcule la fonction gamma de x.Cette fonction tend vers l'infini assez rapidement:

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

Déjà à x == 31, vous avez un très grand nombre.

Lors de la conversion de ce très grand double retour de type int, les résultats sont un comportement indéfini (4.9 Flottant intégrale des conversions [conv.fpint]):

Un prvalue d'un type à virgule flottante peuvent être convertis à un prvalue d'un type entier.La conversion de trun - cates;qui est, la partie fractionnaire est mis au rebut.Le comportement est indéfini si le tronc de la valeur ne peut pas être représenté dans le type de destination.

Sur mon système, cette conversion (avec une entrée de {30, 2}) résultats dans un int avec la valeur de -2147483648.C'est facile à observer par l'insertion de certains des instructions d'impression:

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

qui, pour moi, les résultats:

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

Comme on le voit, l'UB à l'issue d'une division par zéro, ce qui est également UB.Et des sons très semblable au comportement que vous voyez.

La solution à ce problème est de calculer les choses en utilisant seulement partie intégrante de l'arithmétique, et de telle manière que les calculs intermédiaires ne débordent pas si le résultat final est représentable dans le type intégral.Cela implique l'utilisation d'un Plus grand Commun Diviseur de la fonction.

Le code Open source qui est fait, est disponible ici:

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

Recherche pour "count_each_combination".Votre chooseC peut être codée en termes de count_each_combination comme suit:

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

Maintenant chooseC(30, 2) sera de retour 435.Si count_each_combination est incapable de stocker le résultat dans un int, un std::overflow_error sera levée.

Si vous souhaitez limiter votre chooseC pour k == 2, ou peut-être le faire juste temporairement juste pour mieux comprendre l'algorithme, notez que la formule pour le comptage des combinaisons:

Lorsque k == 2, cela se simplifie en:

n*(n-1)/2

Maintenant n est même, ou n-1 est même.Vous pouvez découvrir les, et ensuite diviser ce nombre par 2, sans erreur de troncation, puis multiplier le résultat par le nombre, ce qui n'était pas divisé par 2.Ainsi, vous obtenez le résultat exact, sans possibilité d'erreur de troncation, ni intermédiaire de dépassement, en utilisant seulement partie intégrante de l'arithmétique.C'est la technique qui count_each_combination utilise, mais généralisé à tout diviseur, afin de livrer un résultat qui est toujours exacte si elle peut s'adapter à l'intégrale fournie type.