[Rcpp-devel] Rcpp RMpfr mpreal to mpfr wrapping

[Thell Fowler]

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.  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...

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.

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"),
    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",
    "-Ic:/RtoolsLocal/R-2.15.0/include" )

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"))


Sincerely,
Thell
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20121118/3a75e530/attachment-0001.html>