[Rcpp-devel] Rcpp RMpfr mpreal to mpfr wrapping
Dirk Eddelbuettel
edd at debian.org
Sun Nov 18 20:59:44 CET 2012
Hi Thell,
On 18 November 2012 at 13:03, Thell Fowler wrote:
| Hi all.
|
| I've been doing some testing on how to get arbitrary precision numbers out of
| c++ and into R. The following example inline Rcpp code generates an RcppMpfr
| cxxfunction plugin for wrapping <mpreal> to the Rmpfr S4 "mpfr1" object.
Nice.
| There are local include paths that need adjusting (for now) to your
| environment, and the mpreal header only c++ wrapper for mpfr ( http://
| www.holoborodko.com/pavel/mpfr/ ). I'll be making this into a package but
| wanted to ask about how to do something first...
Just so that I understand: that is another C++ layer on top of GNU mpfr ?
Yet you do link with -lmpfr below...
| Rmpfr defines 'mpfr1' and 'mpfr' where 'mpfr' is usable within R and is a list
| type containing 'mpfr1' objects. So the Rcpp::wrap is setup to get and
| populate an S4 "mpfr1" object, but to be usable on the R side the S4 "mpfr"
| needs to be used. I can get and populate the ".Data" field in the
| mpfr::mpreal wrapper using Rcpp::List but was wondering if the Rcpp::wrap could
| be setup to capture an S4 class type.
"I guess so".
I haven't tried. We do get the whole SEXP, and we have some code to access
and set S4 slots.
| BTW:: The timing comparison shows a huge performance gain (the R_mpfr_fac call
| returns an mpfr1 object only and speed is comparable but slightly slower then
| using Rcpp to do the same) ::
|
| test
| elapsed relative
| 5 factorial(93) 0.05
| 1.0
| 6 .Call(Rmpfr:::R_mpfr_fac, 93, 255) 0.08 1.6
| 2 mpreal_factorial(93, 255) 0.20
| 4.0
| 4 new("mpfr", .Call(Rmpfr:::R_mpfr_fac, 93, 255)) 2.14 42.8
| 3 Rmpfr::factorialMpfr(93, 255) 2.37 47.4
| 1 factorial(mpfr(93, 255)) 13.26
| 265.2
| --
|
| I'd like to have wrap<mpfr::mpreal> wrap to an 'mpfr1' and wrap<S4<mpfr1>>
| (imagined def) wrap to an mpfr list.
|
| suppressMessages( require( inline) )
| suppressMessages( require( Rcpp ) )
| suppressMessages( require( Rmpfr ))
|
| #### Plugin Setup ####
| # For Rcpp wrapping.
|
| include.before <- '
| #include <cmath>
| #include <Iterator>
| #include <RcppCommon.h>
| #include <mpreal.h>
|
| namespace Rcpp {
| template <> SEXP wrap( const mpfr::mpreal& ) ;
| }
| '
|
| include.after <- '
| // Definitions from Rmpfr conversion utilities.
| #if GMP_NUMB_BITS == 32
|
| # define R_mpfr_nr_ints nr_limbs
| # define R_mpfr_exp_size 1
|
| #elif GMP_NUMB_BITS == 64
|
| # define R_mpfr_nr_ints (2*nr_limbs)
| # define R_mpfr_exp_size 2
|
| #endif
|
| # define RIGHT_HALF(_LONG_) ((long long)(_LONG_) & 0x00000000FFFFFFFF)
| # define LEFT_SHIFT(_LONG_) (((long long)(_LONG_)) << 32)
|
| namespace Rcpp {
|
| template<> SEXP wrap( const mpfr::mpreal& v )
| {
| mpfr_srcptr pv( v.mpfr_srcptr() );
|
| S4 vS4("mpfr1");
| vS4.slot("prec") = wrap( (int)pv->_mpfr_prec );
| vS4.slot("sign") = wrap( (int)pv->_mpfr_sign);
|
| IntegerVector d( std::ceil( (double)pv->_mpfr_prec / (double)
| mp_bits_per_limb ) );
| mp_size_t i = 0;
|
| if( GMP_NUMB_BITS == 32 ) {
|
| vS4.slot("exp") = wrap( (int)pv->_mpfr_exp );
|
| for( auto & e : d ) {
| e = (int) pv->_mpfr_d[i];
| ++i;
| }
|
| } else {
|
| IntegerVector exp(2);
| exp[0] = (int) RIGHT_HALF(pv->_mpfr_exp);
| exp[1] = (int) (pv->_mpfr_exp >> 32);
| vS4.slot("exp") = wrap( exp );
|
| for(i=0; i < d.size(); i++) {
| d[2*i] = (int) RIGHT_HALF(pv->_mpfr_d[i]);
| d[2*i+1]= (int) (pv->_mpfr_d[i] >> 32);
| }
|
| }
|
| vS4.slot("d") = wrap( d );
|
| S4 ans("mpfr");
| ans.slot(".Data") = List::create( wrap(vS4) );
|
| return wrap( ans );
| //return wrap( vS4 );
| }
|
| }
| '
|
| #### Plugin Definition ####
| RcppMpfr <- function()
| {
| plugin <- Rcpp:::Rcpp.plugin.maker(
| include.before= include.before,
| include.after= include.after,
| LinkingTo=c('Rmpfr',"Rcpp"),
I don't see Rmpfr exporting header files, so not sure if you need Rmpfr here?
| libs="-Lc:/RtoolsLocal/R-2.15.0/lib/i386 -lmpfr -lgmp" )
| settings <- plugin()
| settings$env$PKG_CPPFLAGS=paste(
| "-std=c++0x",
| "-Ic:/users/almostautomated/src/RcppMpfr++/mpreal/inst/include",
You can probably make that local to your package if you keep it in RcppMpfr++/inst/include/...
Or is this another package? I'm lost here....
| "-Ic:/RtoolsLocal/R-2.15.0/include" )
R will already provide this.
| settings$env$PKG_CXXFLAGS="-IC:/Users/almostautomated/R/win-library/2.15/
| Rcpp/include"
|
| settings
| }
|
| registerPlugin(name="RcppMpfr", plugin=RcppMpfr )
|
|
| #### Sample Function Setup ####
| headers <- '
| #include <iostream>
| '
|
| sources <- '
| using namespace mpfr;
| using namespace Rcpp;
|
| mpreal::set_default_prec( as< unsigned long >( prec ) );
| return wrap( fac_ui( as< double >( x ) ) );
| '
|
| #### Sample Function Defintion ####
| mpreal_factorial <-
| cxxfunction( sig=signature(x='numeric', prec='numeric'),
| includes=headers,
| body=sources,
| plugin="RcppMpfr",
| verbose=FALSE)
|
|
| #### Sample Function Benchmark ####
| timings <-
| benchmark( factorial( mpfr(93,255) ),
| mpreal_factorial(93,255),
| Rmpfr::factorialMpfr(93,255),
| new("mpfr", .Call(Rmpfr:::R_mpfr_fac, 93, 255)),
| factorial(93),
| .Call(Rmpfr:::R_mpfr_fac, 93, 255),
| replications=1000,
| order="elapsed",
| columns=c("test", "elapsed", "relative"))
|
Looks good. I am not quite sure what your question is, though. Your proof of
concept seems to working, so why not build a package ?
Dirk
|
| Sincerely,
| Thell
|
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
--
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com
More information about the Rcpp-devel
mailing list