N Scegli K Funzione si arresta anomentarsi RCPP

Question

Ho scritto una funzione 'n Scegliere K' in C ++, che interfaccia con r via RCPP.Per qualche motivo sto ottenendo un errore di runtime "divide per zero".Succede quando provo a valutare 30 Scegli 2.

Ho provato a valutare manualmente ogni riga (con EvalCPP), e sono ancora perplesso di dove sta accadendo il divario per zero.Forse qualcuno potrebbe indicarmi a me o suggerire un modo migliore per scrivere n scegliere k?

Ecco il codice:

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

.

• Non c'è TGAMMA in STD per quanto posso vedere

• è come una funzione choose, quindi farei solo ciò che è sotto

• r ha anche la distribuzione gamma ecc. Quindi puoi farlo a mano pure

• Perché non hai appena stampato i valori factN, factK, factDiff?

Semplice soluzione RCPP:

#include <Rcpp.h>

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

.

Esempio:

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

.

Modifica: Seguendo il commento di @Blastfurnace Informazioni su tgamma() nell'intestazione GeneracoDetaGcode C ++ 11, ecco una versione riparata che funziona bene per me:

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

.

Esempio di esempio:

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

Quindi std::tgamma(x) calcola la funzione gamma di x. Questa funzione va a infinito abbastanza rapidamente:

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

Già a X== 31, hai un numero molto grande.

Quando si convertono questo doppio grande schienale in int, i risultati sono un comportamento indefinito (4.9 conversioni flottanti-integrali [conv.fppint]):

.

Un pruralue di un tipo di punto flottante può essere convertito in una prvalue di un Tipo intero. La conversione Trunjates; cioè, la parte frazionata è scartato. Il comportamento non è definito se il valore troncato non può essere rappresentato nel tipo di destinazione.

Sul mio sistema questa conversione (con un ingresso di {30, 2}) si traduce in int con il valore -2147483648. Questo è facilmente osservato inserendo alcune dichiarazioni di stampa:

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

.

che per me uscirà:

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

.

Come si può vedere, l'UB infine si traduce in un divario per zero, che è anche UB. E sembra molto simile al comportamento che stai vedendo.

La soluzione a questo problema è quella di calcolare le cose utilizzando solo l'aritmetica integrale e in modo tale che i calcoli intermedi non eccessivano se il risultato finale è rappresentato nel tipo integrale. Ciò comporta l'uso di una maggiore funzione di divisor comune.

Codice open source che è disponibile qui:

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

Ricerca di "conte_each_combination". Il tuo chooseC può essere codificato in termini di count_each_combination come:

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

.

Ora chooseC(30, 2) restituirà 435. Se count_each_combination non è in grado di memorizzare il risultato in un int, verrà lanciato un std::overflow_error.

Se si desidera vincolare il chooseC a k == 2, o forse farlo solo temporaneamente solo per capire meglio l'algoritmo, nota che la formula per le combinazioni di conteggio è:

Quando k == 2, questo semplifica a:

n*(n-1)/2

.

Ora n è pari o n-1 è pari. È possibile scoprire quale, quindi dividere quel numero di 2, senza errori di troncamento, e quindi moltiplicare il risultato del numero che non è stato diviso per 2. Quindi si ottiene il risultato esatto senza possibilità di errore di troncamento, né overflow intermedio , utilizzando solo aritmetico integrale. Questa è la tecnica che genera count_each_combination, ma generalizzato a qualsiasi divisore, per fornire un risultato sempre esatto se può adattarsi al tipo integrale in dotazione.