0

I am a C++ beginner hoping to find some help here.

To motivate my question, I wish to write a function that performs convolutions in C++ via either infinite precision or rational arithmetic. Since I wish to call this function from R (via Rcpp), I started from the following approach (which is based on Dirk's example).

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector convolve5cpp(Rcpp::NumericVector & xa,
                                 Rcpp::NumericVector & xb) {
  int n_xa = xa.size();
  int n_xb = xb.size();
  Rcpp::NumericVector xab(n_xa + n_xb - 1);
  
  Rcpp::Range r(0, n_xb-1);
  for (int i=0; i<n_xa; i++, r++) {
    xab[r] += xa[i] * xb;
  }
  return xab;
}

Realising I have to cast an infinite-precision rational type for my C++ variables, I next try to perform type re-casting, relying on boost libraries. Starting with a simple example involving scalars (shown below)...

// [[Rcpp::depends(BH)]]
namespace mp = boost::multiprecision;

mp::mpq_rational divideQ(int a, int b) {
  mp::mpq_rational c = a / b;
  return c;
}

I hit an error "cannot convert type to SEXP".

I am lost at this point. I am aware of the gmp R package providing bigq class objects in R, but don't know if this can be leveraged to convert mpq_rationals into SEXP.

My questions are thus:

  • How to convert a boost library rational type to a valid R object? (Ideally with bigq class)
  • Does my reasoning make sense, or are there other approaches to this? (For example, does there already exist an implementation of convolve5cpp above for vectors that performs internal operations with infinite-precision?)

Many thanks!


(30th Dec 2021 Update)

Thanks to all who left helpful comments below. After some trial and error, I went ahead and tried locally editing some source files in the gmp library (available here).

I wrote the following function, which closely follows the example by Dirk above. I mirrored the structure of operator+, another function defined for addition of two bigrationals.

/**
 * \brief Return  conv( a, b )
 */
bigrational operatorCONV(const bigrational& lhs, const bigrational& rhs)
{
  /* get sizes of each vector */
  int n_xa = lhs.size();
  int n_xb = rhs.size();

  /* initialize big rational of size */
  bigrational xab(n_xa + n_xb - 1);
  
  /* fill in the entries */
  Rcpp::Range r(0, n_xb-1);
  for (int i=0; i<n_xa; i++, r++) {
    xab[r] = bigrationalR::create_bigrational(xab[r], bigrationalR::create_bigrational(lhs[i], rhs, mpq_mul), mpq_add);
  }

  /* return */
  return xab;
}

This didn't work. When building the package, I ran into the following errors, which look like they fall under two broad categories.

  1. No member named 'size' in 'bigrational': int n_xa = lhs.size().
  2. Type 'const bigrational' does not provide a subscript operator: bigrationalR::create_bigrational(lhs[i], rhs, mpq_mul)
  3. Use of undeclared identifier 'Rcpp': Rcpp::Range r(0, n_xb-1)
  4. Use of undeclared identifier 'r': bigrationalR::create_bigrational(xab[r],...

Namely, 1 and 2 arise because I treat lhs and rhs as vectors of bigrationals; 3 and 4 arise because I use an object recognized only by Rcpp that wasn't included in the file.

How would one go about addressing the two broad problems? For 1 and 2, I do not understand how add.bigq (which calls bigrational_add, a function that takes in two bigrationals) can take in two vectors but somehow a convolution function like the above is not "flexible" like that. For 3 and 4, is there a stdlib alternative to Rcpp::Range that can be called relatively simply, or do I have to create my own object type?

Thanks again!

  • Briefly, it doesn't (autoMAGICally) work that way. There is no (automatic) conversion (programmed) for `mp::mpq_rational` to `SEXP`, and the compiler tells you so. You _can_ however add those (as we say how in the Rcpp documentation) but that is not a beginner's task. – Dirk Eddelbuettel Dec 12 '21 at 01:36
  • I see. Thanks for writing back, @DirkEddelbuettel. Could you point me to a simple example or where this is mentioned specifically in the documentation, so that I can try to build on it? – youngtableaux Dec 12 '21 at 01:43
  • 2
    Most specifically the ['extending Rcpp' vignette](https://cloud.r-project.org/web/packages/Rcpp/vignettes/Rcpp-extending.pdf); the problem is also (more gingerly) discussed with some examples at the [Rcpp Gallery](https://gallery.rcpp.org). – Dirk Eddelbuettel Dec 12 '21 at 01:53
  • 2
    Specific examples I've found usefule for this are [here](https://gallery.rcpp.org/articles/custom-templated-wrap-and-as-for-seamingless-interfaces/) and [here](https://gallery.rcpp.org/articles/custom-as-and-wrap-example/), and I've had a go at implementing custom wrappers [here](https://github.com/SymbolixAU/mapboxGeometry/blob/master/inst/include/mapboxgeometry.hpp) (among other examples) – SymbolixAU Dec 17 '21 at 01:18
  • Thanks for both your suggestions. – youngtableaux Dec 30 '21 at 22:48

0 Answers0