N Wählen Sie K Funktion stürzt Rcpp ab

Question

Ich habe in C++ eine Funktion „n Choose k“ geschrieben, die über Rcpp mit R verbunden ist.Aus irgendeinem Grund erhalte ich den Laufzeitfehler „Teile durch Null“.Es passiert, wenn ich versuche, 30 auszuwerten und 2 auszuwählen.

Ich habe versucht, jede Zeile manuell auszuwerten (mit evalCpp), und bin immer noch verwirrt, wo die Division durch Null stattfindet.Vielleicht könnte mich jemand darauf hinweisen oder eine bessere Schreibweise vorschlagen, um K zu wählen?

Hier ist der 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

Knapp:

• Soweit ich sehen kann, gibt es in std kein Tgamma

• R selbst als choose Funktion, also würde ich einfach das tun, was unten steht

• R verfügt auch über die Gammaverteilung usw., sodass Sie dies auch manuell tun können

• Warum haben Sie die Werte nicht einfach ausgedruckt? factN, factK, factDiff ?

Einfache Rcpp-Lösung:

#include <Rcpp.h>

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

Beispiel:

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

Bearbeiten: Im Anschluss an den Kommentar von @Blastfurnace über tgamma() im C++11 cmath Header, hier ist eine reparierte Version, die für mich gut funktioniert:

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

Beispielanwendung:

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

Also std::tgamma(x) berechnet die Gammafunktion von x.Diese Funktion geht ziemlich schnell gegen Unendlich:

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

Bereits bei x == 31 haben Sie eine sehr große Zahl.

Bei der Konvertierung dieses sehr großen Doubles zurück in int kommt es zu undefiniertem Verhalten (4.9 Floating-Integral-Konvertierungen [conv.fpint]):

Ein Prvalue eines schwebenden Punkttyps kann in einen Prorue eines Ganzzahltyps umgewandelt werden.Die Konvertierung wird abgeschnitten;Das heißt, der Bruchteil wird verworfen.Das Verhalten ist undefiniert, wenn der abgeschnittene Wert im Zieltyp nicht dargestellt werden kann.

Auf meinem System führt diese Konvertierung (mit einer Eingabe von {30, 2}) zu einem int mit dem Wert -2147483648.Dies lässt sich leicht beobachten, indem man einige print-Anweisungen einfügt:

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

was für mich Folgendes ausgibt:

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

Wie man sieht, ergibt sich aus der UB letztlich eine Division durch Null, die ebenfalls UB ist.Und klingt dem Verhalten, das Sie sehen, sehr ähnlich.

Die Lösung dieses Problems besteht darin, Dinge ausschließlich mit Integralarithmetik zu berechnen, und zwar so, dass die Zwischenberechnungen nicht überlaufen, wenn das Endergebnis im Integraltyp darstellbar ist.Dies erfordert die Verwendung einer Funktion des größten gemeinsamen Teilers.

Open-Source-Code, der dies tut, ist hier verfügbar:

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

Suchen Sie nach „count_each_combination“.Dein chooseC kann in Bezug auf codiert werden count_each_combination etwa so:

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

Jetzt chooseC(30, 2) gibt 435 zurück.Wenn count_each_combination ist nicht in der Lage, das Ergebnis in einem zu speichern int, A std::overflow_error wird geworfen.

Wenn Sie Ihre einschränken möchten chooseC Zu k == 2, oder vielleicht nur vorübergehend tun, um den Algorithmus besser zu verstehen, beachten Sie, dass die Formel zum Zählen von Kombinationen lautet:

Wann k == 2, dies vereinfacht sich zu:

n*(n-1)/2

Jetzt auch nicht n ist gerade, oder n-1 ist gerade.Sie können herausfinden, was das ist, und diese Zahl dann ohne Kürzungsfehler durch 2 dividieren und dann das Ergebnis mit der Zahl multiplizieren, die nicht durch 2 geteilt wurde.Somit erhalten Sie das exakte Ergebnis ohne die Möglichkeit eines Kürzungsfehlers oder eines zwischenzeitlichen Überlaufs, da Sie nur die Integralarithmetik verwenden.Dies ist die Technik, die count_each_combination verwendet, jedoch auf einen beliebigen Teiler verallgemeinert, um ein Ergebnis zu liefern, das immer genau ist, wenn es in den angegebenen Integraltyp passen kann.