[Rcpp-devel] Numerical precision in rotations with Eigen

[Rafael Ayala Hernandez]

Hi,

I have implemented a function to rotate a 3D vector a given angle around a given axis (basically wrapping the functionality provided by Eigen::AngleAxis) as an Rcpp function.
Below is an extract from the source file:

#include <Rcpp.h>
#include <RcppEigen.h>
#include <Eigen/Eigen>
#include <Eigen/Geometry>
using namespace Rcpp;
using Rcpp::as;
using Eigen::Map;
using Eigen::VectorXd;
using Eigen::Vector3d;
using Eigen::Matrix;
using Eigen::MatrixXd;
using Eigen::Matrix3d;
using Eigen::Dynamic;
using Eigen::AngleAxisd;

// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
NumericVector rotate3DVectorAngleAxis(NumericVector x, NumericVector axis, double angle) {
    // Note: x should be 1 single vector to rotate
    Map<VectorXd> xEigen(as<Map<VectorXd> >(x));
    Map<VectorXd> axisEigen(as<Map<VectorXd> >(axis));
    Matrix3d rotation;
    rotation = AngleAxisd(angle, axisEigen.normalized());
    Vector3d rotatedVector;
    rotatedVector = rotation.matrix() * xEigen;
    return wrap(rotatedVector);
}


The function executes with no problems and gives the expected output.

However, I have noticed that there is times where a value of 0 would be expected for some vector components, but instead a very small value is returned. For example, executing the following on my machine:

rotate3DVectorAngleAxis(c(0,0,1), c(1,0,0), pi/2)

Which represents a rotation of 90 degrees around the X axis, and therefore the expected output value would be c(0,1,0)

However, it instead results in an output vector of c(0.000000e+00, -1.000000e+00,  6.123234e-17)

I.e., the value for the Z component is 6.123234e-17. 

I guess this is somehow related to machine precision, but is there an exact cause that could be possibly fixed?
Additionally, why does the deviation from 0 only seem to happen for the Z component, and not for the X component?
Varying the amount of rotation also seems to lead to a cumulative error on the same component. E.g., rotating by 4001 times pi/2 (i.e., 1000 complete 360 degrees rotations plus a 90 degree rotation):

rotate3DVectorAngleAxis(c(0,0,1), c(1,0,0), 14401*pi/2)

results in c(0.000000e+00, -1.000000e+00, 4.776933e-13)

Which seems again to be accumulating in the same Z component.

Is there anything I could improve in my implementation to avoid these numerical errors?

For reference, .Machine$double.eps on my machine is 2.220446e-16

Thanks a lot in advance

Best wishes,

Rafa