N Escolha K Função trava Rcpp

Question

Eu escrevi uma função 'n escolha k' em C++, que faz interface com R via Rcpp.Por algum motivo, estou recebendo um erro de tempo de execução 'dividir por zero'.Acontece quando tento avaliar 30 e escolher 2.

Tentei avaliar cada linha manualmente (com evalCpp) e ainda estou confuso sobre onde está acontecendo a divisão por zero.Talvez alguém possa me apontar isso ou sugerir uma maneira melhor de escrever e escolher K?

Aqui está o 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:

• Não há tgamma em std, pelo que posso ver

• O próprio R como um choose função, então eu faria apenas o que está abaixo

• R também tem a distribuição gama, etc., então você também pode fazer isso manualmente

• Por que você simplesmente não imprimiu os valores factN, factK, factDiff ?

Solução Rcpp simples:

#include <Rcpp.h>

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

Exemplo:

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

Editar: Seguindo o comentário de @Blastfurnace sobre tgamma() no C++11 cmath cabeçalho, aqui está uma versão reparada que funciona bem para mim:

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

Exemplo de uso:

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

Então std::tgamma(x) calcula a função gama de x.Esta função vai para o infinito rapidamente:

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

Já em x == 31, você tem um número muito grande.

Ao converter este duplo muito grande de volta para int, os resultados são um comportamento indefinido (4.9 Conversões de integrais flutuantes [conv.fpint]):

Um prvalue de um tipo de ponto flutuante pode ser convertido em um prvalor de um tipo inteiro.A conversão é truncada;Ou seja, a parte fracionária é descartada.O comportamento é indefinido se o valor truncado não puder ser representado no tipo de destino.

No meu sistema esta conversão (com uma entrada de {30, 2}) resulta em um int com o valor -2147483648.Isso é facilmente observado inserindo algumas instruções de impressão:

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 mim produz:

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

Como pode ser visto, o UB resulta em uma divisão por zero, que também é UB.E parece muito semelhante ao comportamento que você está vendo.

A solução para este problema é calcular coisas usando apenas aritmética integral, e de tal forma que os cálculos intermediários não transbordem se o resultado final for representável no tipo integral.Isso implica o uso de uma função do Máximo Divisor Comum.

O código-fonte aberto que faz isso está disponível aqui:

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

Procure por "count_each_combination".Seu chooseC pode ser codificado em termos de count_each_combination igual a:

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

Agora chooseC(30, 2) retornará 435.Se count_each_combination não consegue armazenar o resultado em um int, a std::overflow_error será lançado.

Se você quiser restringir seu chooseC para k == 2, ou talvez faça isso apenas temporariamente, apenas para entender melhor o algoritmo, observe que a fórmula para contar combinações é:

Quando k == 2, isso simplifica para:

n*(n-1)/2

Agora também n é par, ou n-1 é par.Você pode descobrir qual e depois dividir esse número por 2, sem erro de truncamento, e então multiplicar o resultado pelo número que não foi dividido por 2.Assim você obtém o resultado exato sem possibilidade de erro de truncamento, nem overflow intermediário, usando apenas aritmética integral.Esta é a técnica que count_each_combination usa, mas generalizado para qualquer divisor, para fornecer um resultado que é sempre exato se puder caber no tipo integral fornecido.