(Rcpp)armadillos abs() function outputs false values when used with c++ double while std::abs works

Consider this R function:

r_abs <- function(x,y,z){
    2 * abs((x >= y) - z)
}

And those 2 RcppArmadillo equivalents:

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppProgress)]]
#include <RcppArmadillo.h>
#include <string>
#include <progress.hpp>
using namespace arma;

// [[Rcpp::export]]
double arma_abs(const double &x, const double &y, const double &z)
{
    return 2 * abs((x >= y) - z);
}

// [[Rcpp::export]]
double std_abs(const double &x, const double &y, const double &z)
{
    return 2 * std::abs((x >= y) - z);
}

Now here is the problem:

sourceCpp(".test/test.cpp")

x <- 1
y <- 2
z <- 0.5

r_abs(x,y,z)
[1] 1
arma_abs(x, y, z)
[1] 0
std_abs(x, y, z)
[1] 1

Why does arma_abs() output 0 here? The problem only occurs for -1 < z < 1.

I really appreciate any help you can provide.

Solution

This is a tricky and somewhat tedious issue with C and C++ partially fixed by the corred std::abs() added more recently.

If you're on a Unix system try man 3 abs. Your use of the (unqualified) abs() is equivalent to ::abs() and it gets int abs(int) from stdlib.h.

And the cast to int makes the argument -0.5 truncate to zero which times two is still zero. In short, arma has nothing to do with this. (And Armadillo generally only has functions for its vector, matrix, cube, ... types and not atomic C/C++ types like double).

For completeness, my version of your code with the R part embedded below. I removed the unused references RcppProgress.

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

//using namespace arma;

// [[Rcpp::export]]
double arma_abs(const double &x, const double &y, const double &z) {
  auto val = (x >= y) - z;
  Rcpp::Rcout << "val " << val << std::endl;
  Rcpp::Rcout << "int(val) " << int(val) << std::endl;
  return 2 * ::abs((x >= y) - z);
}

// [[Rcpp::export]]
double std_abs(const double &x, const double &y, const double &z) {
    return 2 * std::abs((x >= y) - z);
}

/*** R

r_abs <- function(x,y,z){
    2 * abs((x >= y) - z)
}

x <- 1
y <- 2
z <- 0.5

r_abs(x,y,z)           # 1
arma_abs(x, y, z)      # 0
std_abs(x, y, z)       # 1
*/

Edit: Typo corrected. Above I wrote std::fabs() when I meant std::abs(). fabs() is what we always had to use to avoid the very issue of a cast to int seen here.