Rcpp: Bug: gamma function presenting strange values

[nettoyoussef]

First of all, I would like to thank you for this amazing package.
I am pretty sure I would have a hard time using the R's standard API to develop C/C++ code.

Bug Report

Not sure why, but the gamma function seems to be presenting values that contrast with R.
Check the MRE below:

// [[Rcpp::export]]
void check_gamma(double z){
  Rprintf("gamma: %f\n", gamma(z));
  Rprintf("Rf_gammafn: %f\n", Rf_gammafn(z));
}

Then, in a R session:

check_gamma(5.5)
# gamma: 3.957814
# Rf_gammafn: 52.342778

Alternatively, using the cppFunction:

cppFunction("
double check_gamma2(double z){
  return gamma(z);
}
")

check_gamma2(5.5)
#> [1] 3.957814
gamma(5.5) # base R version
#> [1] 52.34278

My session info:

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.utf8       LC_NUMERIC=C             
 [3] LC_TIME=en_US.utf8        LC_COLLATE=en_US.utf8    
 [5] LC_MONETARY=en_US.utf8    LC_MESSAGES=C            
 [7] LC_PAPER=en_US.utf8       LC_NAME=C                
 [9] LC_ADDRESS=C              LC_TELEPHONE=C           
[11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C      

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] Rcpp_1.0.3

loaded via a namespace (and not attached):
[1] compiler_3.6.1 tools_3.6.1   

[eddelbuettel]

Be careful -- I just got in from a run and don't _this minute_ have time to look into this but that is IIRC a known and documented issue. _Our interface_ follows the _C API interface of R_ which may have different parameterisations than the _R interface_. Look at the header; I believe we have (old) prior tickets here and at SO. I can look later if you find nothing.

[nettoyoussef]

Hi Dirk, thanks for the super fast response!

Yes, I am aware of those different parameterizations, but I am using above the gamma function, not the distribution.

A strange thing happens when I change the type of input.
Maybe I am doing something wrong?

# using a numeric vector works
cppFunction("
Rcpp::NumericVector check_gamma2(Rcpp::NumericVector z){
  return gamma(z);
}
")

check_gamma2(5.5)
#> [1] 52.34278

#using a double it produces the unexpect result
cppFunction("
double check_gamma2(double z){
  return gamma(z);
}
")
check_gamma2(5.5)
#> [1] 3.957814

If this is a problem of type coercion, and I just messed things up, feel free to close this issue.
I thought that Rcpp::wrap would work automagically in such context :) .

[eddelbuettel]

Sorry I am still between things. Does this help?

> Rcpp::cppFunction("NumericVector g(NumericVector x) { return Rcpp::gamma(x); }")
> gamma(seq(-3,3) + 0.001)   # simpler version of what help(gamma) has in examples
[1]  -166.876402   500.462330 -1000.424197   999.423772     0.999424     1.000423     2.001847
> g(seq(-3,3) + 0.001)
[1]  -166.876402   500.462330 -1000.424197   999.423772     0.999424     1.000423     2.001847
> 

[eddelbuettel]

The key here is that double gamma(double) is not us but the system libraries. Hence the Rf_ prefix for R, and/or the R:: namespace for the ones we reexport via the Rcpp package. Be careful with flattening namespaces, I usually recommend against. On my (Ubuntu 20.04) system, man 3 gamma leads to some warnings and a recommendation to use tgamma or lgamma instead.

[nettoyoussef]

Thanks a lot, Dirk, you are of course, correct!
Sorry to disturb you with this, I am closing this issue now.

[eddelbuettel]

No worries we've all been there. I cannot tell you how often (under an older C++ standard) I mistakenly called abs() which wasn't even defined for double (but fabs()) was. In short: mind the signature, and pay attention to the namespace.