Regression with Armadillo 15.0.0 -- and 14.6.* under Intel

[eddelbuettel]

The eBsc package package is experiencing a very strange regression when going from Armadillo 14.6.* (either .0 or .3) to 15.0.0. So this looks like a change to worse. Oddly, the same error also appears under the Intel compiler where (copying here in case the link disappears) we currently see

* checking examples ... [9s/11s] ERROR
Running examples in ‘eBsc-Ex.R’ failed
The error most likely occurred in:

> ### Name: plot.eBsc
> ### Title: Plot fitted components
> ### Aliases: plot.eBsc
> ### Keywords: plot
> 
> ### ** Examples
> 
> 
> library(eBsc)
> n <- 250
> sigma <- 0.05
> Basis <- list()
> for(i in 1:6) Basis[[i]] <- drbasis(nn = n, qq = i)
Error: rank(): failed
Execution halted

Here the example for plot.eBsc fails with insuffcient rank(). We can trace this back to (R function) drbasis() which calls C++ function drbasis(). For arguments nn=250, qq=5 this fails deterministically.

We isolated the core of the qq=5 code (which is rather neatly inside a case statement) having it return the matrix for which the rank() fails, and have it print the final ten column of the first ten rows. The resulting C++ file has zero dependencies on R or R packages. No changes were made apart from reindentation, and removing of extra code not needed for a minimally viable replicable example as provided here.

For convenience we just point at the include directories provided by what is currently PR #473 which has both Armadillo 15.0.0 (for compilations under C++14 or newer) and Armadillo 14.6.3 (for C++11 compilations as more than two hundred CRAN packages still hardcode).

Note that the following example can be reproduced anywhere with both Armadillo 15.0.0 or equivalent, and the older Armadillo 14.6.* by adjusting the -I... location. current is the released upstream version 15.0.0, legacy the preceding 14.6.3. We do not need to link with -larmadillo or -lblas or -llapack etc as we execute no such functionality as only vector / matrix code is used.

edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ g++ -Wall -pedantic \   # indented for display
    -I/usr/local/lib/R/site-library/RcppArmadillo/include/legacy  \
    -o ebsc_rank_issue_main ebsc_rank_issue_main.cpp -lm 
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ ./ebsc_rank_issue_main 
M[0:9, 240:249]
   0.2000   0.2000   0.2000   0.2000   0.2000   0.2000   0.2000   0.2000   0.2000   0.2000
   0.0698   0.0704   0.0710   0.0717   0.0722   0.0728   0.0734   0.0739   0.0744   0.0749
  -0.0796  -0.0804  -0.0811  -0.0818  -0.0823  -0.0829  -0.0834  -0.0838  -0.0841  -0.0844
   0.0803   0.0816   0.0829   0.0840   0.0849   0.0858   0.0865   0.0871   0.0876   0.0879
  -0.0759  -0.0781  -0.0802  -0.0821  -0.0838  -0.0852  -0.0865  -0.0875  -0.0883  -0.0888
   0.0688   0.0722   0.0754   0.0783   0.0809   0.0831   0.0850   0.0866   0.0878   0.0887
  -0.0602  -0.0650  -0.0695  -0.0735  -0.0772  -0.0804  -0.0831  -0.0854  -0.0871  -0.0884
   0.0505   0.0568   0.0626   0.0680   0.0729   0.0772   0.0809   0.0839   0.0863   0.0881
  -0.0397  -0.0476  -0.0550  -0.0618  -0.0680  -0.0735  -0.0783  -0.0823  -0.0854  -0.0876
   0.0283   0.0377   0.0467   0.0550   0.0626   0.0695   0.0754   0.0804   0.0843   0.0871
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ g++ -Wall -pedantic \   # indented for display
    -I/usr/local/lib/R/site-library/RcppArmadillo/include/current  \
    -o ebsc_rank_issue_main ebsc_rank_issue_main.cpp -lm 
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ ./ebsc_rank_issue_main 
M[0:9, 240:249]
      nan      nan      nan      nan      nan      nan      nan      nan      nan      nan
   0.0698      nan      nan      nan      nan      nan      nan      nan      nan      nan
  -0.0796  -0.0804      nan      nan      nan      nan      nan      nan      nan      nan
   0.0803   0.0816   0.0829      nan      nan      nan      nan      nan      nan      nan
  -0.0759  -0.0781  -0.0802  -0.0821      nan      nan      nan      nan      nan      nan
   0.0688   0.0722   0.0754   0.0783   0.0809      nan      nan      nan      nan      nan
  -0.0602  -0.0650  -0.0695  -0.0735  -0.0772  -0.0804      nan      nan      nan      nan
   0.0505   0.0568   0.0626   0.0680   0.0729   0.0772   0.0809      nan      nan      nan
  -0.0397  -0.0476  -0.0550  -0.0618  -0.0680  -0.0735  -0.0783  -0.0823      nan      nan
   0.0283   0.0377   0.0467   0.0550   0.0626   0.0695   0.0754   0.0804   0.0843      nan
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ 

The self-contained C++ code is below.

// The eBsc package on CRAN (at version 4.17) reports an error under the Intel compiler.
// This error now also appears with RcppArmadillo 15.0.0. In the example for plot.eBsc,
// the function drbasis() is called with n=250 and a second parameter p that varies from
// 1 to 6. In the case of 5, NA values return leading an rank failure. The drbasis()
// function calls an underlying C function drbasis from file drbasisC.cpp which has case
// statements for the differenct values of p.  We extrace the one for 5 here.

#include <armadillo>

arma::mat compute(int nn) {
    using namespace arma;
    using namespace std;
    const double E = exp(1);
    arma::vec t(nn);
    arma::vec k(nn);
    typedef complex<double> dcomp;
    double Pi = M_PI;
    complex<double> comp0 (0,1);
    complex<double> comp1 (-1,0);
    complex<double> comp2 (0,-0.5);
    complex<double> comp3 (1,-1);
    complex<double> comp4 (-1,1);
    complex<double> comp5 (0,0.5);
    complex<double> comp6 (0,0.25);
    int const0 = 1;
    double const1 = -sqrt(3);
    double const2 = 2*sqrt(3);
    double const3 = -sqrt(5);
    double const4 = 6*sqrt(5);
    double const5 = -6*sqrt(5);
    double const6 = -sqrt(7);
    double const7 = 12*sqrt(7);
    double const8 = -30*sqrt(7);
    double const9 = 20*sqrt(7);
    double const10 = 3;
    double const11 = -60;
    double const12 = 270;
    double const13 = -420;
    double const14 = 210;
    /*
    double const15 = -sqrt(11);
    double const16 = 30*sqrt(11);
    double const17 = -210*sqrt(11);
    double const18 = 560*sqrt(11);
    double const19 = 630*sqrt(11);
    double const20 = 252*sqrt(11);
    */

    for (int j=0;j<nn;j++) {
        t(j)= j/double((nn-1));
    }
    for (int j=0;j<nn;j++) {
        k(j)= j+1;
    }

    arma::mat M5(nn,nn);
    for (int i=0;i<nn; i++) {
        for (int j=0;j<nn;j++) {
            if (j==0) {
                M5(i,j)=const0;
            } else if (j==1) {
                M5(i,j)= const1 + const2*t(i);
            } else if (j==2) {
                M5(i,j)= const3 + const4*t(i) + const5* pow(t(i),2);
            } else if (j==3) {
                M5(i,j)= const6 + const7*t(i) + const8* pow(t(i),2) + const9* pow(t(i),3);
            } else if (j==4) {
                M5(i,j)= const10 + const11*t(i) + const12* pow(t(i),2) + const13* pow(t(i),3) + const14* pow(t(i),4);
            } else {
                dcomp ans = (sqrt(2)*(1 + sqrt(5)/2) -
                             (comp5*(sqrt(10 - 2*sqrt(5)) +
                                     sqrt(2*(5 + sqrt(5)))))/sqrt(2))*(pow(comp1,1 + k(j))/pow(E,pow(comp1,0.1)* (dcomp(-3 + k(j)))*Pi*(1 - t(i))) + pow(E,-(pow(comp1,0.1)*(dcomp(-3 + k(j)))* Pi*t(i)))) +
                    (-(1/sqrt(2)) -
                     (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/
                     sqrt(2))*(pow(comp1,1 + k(j))/pow(E,pow(comp1,0.3)*(dcomp (-3 + k(j)))*Pi*(1 - t(i))) + pow(E,-(pow(comp1,0.3)*(dcomp (-3 + k(j)))*Pi*t(i)))) +
                    (sqrt(2)*(1 + sqrt(5)/2) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
                    *(pow(comp1,1 + k(j))/ pow(E,(sqrt(0.625 + sqrt(5)/8) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k(j)))*Pi*
                                               (1 - t(i))) +pow(E,-((sqrt(0.625 + sqrt(5)/8.) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k(j)))*
                                                                    Pi*t(i)))) +(-(1/sqrt(2)) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
                    *(pow(comp1,1 + k(j))/pow(E,(dcomp(sqrt(0.625 - sqrt(5)/8)) - comp6*(1 + sqrt(5)))*
                                              (dcomp (-3 + k(j)))*Pi*(1 - t(i))) +pow(E,-((sqrt(0.625 - sqrt(5)/8) - comp6*(1 + sqrt(5)))* (dcomp (-3 + k(j)))*Pi*t(i))))- sqrt(2)*cos((-3 + k(j))*Pi*t(i));
                M5(i,j) =ans.real();
            }
            M5(i,j) = M5(i,j)/sqrt(nn);
        }
    }
    return M5;
}

int main(int argc, char *argv[]) {
    constexpr int nn = 250;
    arma::mat M = compute(nn);
    M.submat(0, 240, 9, 249).print("M[0:9, 240:249]");
    exit(0);
}

I have not found any compiler switches, or compilation standard changes, affecting this. As the result also appears under the Intel compiler for the older Armadillo, I suspect the package has code with a latent frailty revealed by either the Intel compiler, or the newer Armadillo. The eBsc maintainer has been contacted by email a few days ago, I have yet to hear back.

[eddelbuettel]

However, when we further reduce the example to not include any Armadillo headers, the issue goes away. That would point to Armadillo being involved.

edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ g++ -o issue474_ebsc_rank_main_no_arma issue474_ebsc_rank_main_no_arma.cpp -lm
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ ./issue474_ebsc_rank_main_no_arma 
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ 

For simplicity, this computes just the first row of the matrix (which failed rank()) and prints just the final ten elements of that first row.

Code below.

// The eBsc package on CRAN (at version 4.17) reports an error under the Intel compiler.
// This error now also appears with RcppArmadillo 15.0.0. In the example for plot.eBsc,
// the function drbasis() is called with n=250 and a second parameter p that varies from
// 1 to 6. In the case of 5, NA values return leading an rank failure. The drbasis()
// function calls an underlying C function drbasis from file drbasisC.cpp which has case
// statements for the differenct values of p.  We extrace the one for 5 here.

#include <cmath>
#include <complex>
#include <iostream>
#include <vector>

std::vector<double> compute(int nn) {
    using namespace std;
    const double E = exp(1);
    vector<double> t(nn);
    vector<double> k(nn);
    typedef complex<double> dcomp;
    double Pi = M_PI;
    complex<double> comp0 (0,1);
    complex<double> comp1 (-1,0);
    complex<double> comp2 (0,-0.5);
    complex<double> comp3 (1,-1);
    complex<double> comp4 (-1,1);
    complex<double> comp5 (0,0.5);
    complex<double> comp6 (0,0.25);
    int const0 = 1;
    double const1 = -sqrt(3);
    double const2 = 2*sqrt(3);
    double const3 = -sqrt(5);
    double const4 = 6*sqrt(5);
    double const5 = -6*sqrt(5);
    double const6 = -sqrt(7);
    double const7 = 12*sqrt(7);
    double const8 = -30*sqrt(7);
    double const9 = 20*sqrt(7);
    double const10 = 3;
    double const11 = -60;
    double const12 = 270;
    double const13 = -420;
    double const14 = 210;
    /*
    double const15 = -sqrt(11);
    double const16 = 30*sqrt(11);
    double const17 = -210*sqrt(11);
    double const18 = 560*sqrt(11);
    double const19 = 630*sqrt(11);
    double const20 = 252*sqrt(11);
    */

    for (int j=0;j<nn;j++) {
        t[j]= j/double((nn-1));
    }
    for (int j=0;j<nn;j++) {
        k[j]= j+1;
    }

    vector<double> M5(nn);
    constexpr int i = 0;
    for (int j=0;j<nn;j++) {
        if (j==0) {
            M5[j] = const0;
        } else if (j==1) {
            M5[j] = const1 + const2*t[i];
        } else if (j==2) {
            M5[j] = const3 + const4*t[i] + const5* pow(t[i],2);
        } else if (j==3) {
            M5[j] = const6 + const7*t[i] + const8* pow(t[i],2) + const9* pow(t[i],3);
        } else if (j==4) {
            M5[j] = const10 + const11*t[i] + const12* pow(t[i],2) + const13* pow(t[i],3) + const14* pow(t[i],4);
        } else {
            dcomp ans = (sqrt(2)*(1 + sqrt(5)/2) -
                         (comp5*(sqrt(10 - 2*sqrt(5)) +
                                 sqrt(2*(5 + sqrt(5)))))/sqrt(2))*(pow(comp1,1 + k[j])/pow(E,pow(comp1,0.1)* (dcomp(-3 + k[j]))*Pi*(1 - t[i])) + pow(E,-(pow(comp1,0.1)*(dcomp(-3 + k[j]))* Pi*t[i]))) +
                (-(1/sqrt(2)) -
                 (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/
                 sqrt(2))*(pow(comp1,1 + k[j])/pow(E,pow(comp1,0.3)*(dcomp (-3 + k[j]))*Pi*(1 - t[i])) + pow(E,-(pow(comp1,0.3)*(dcomp (-3 + k[j]))*Pi*t[i]))) +
                (sqrt(2)*(1 + sqrt(5)/2) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
                *(pow(comp1,1 + k[j])/ pow(E,(sqrt(0.625 + sqrt(5)/8) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k[j]))*Pi*
                                           (1 - t[i])) +pow(E,-((sqrt(0.625 + sqrt(5)/8.) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k[j]))*
                                                                Pi*t[i]))) +(-(1/sqrt(2)) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
                *(pow(comp1,1 + k[j])/pow(E,(dcomp(sqrt(0.625 - sqrt(5)/8)) - comp6*(1 + sqrt(5)))*
                                          (dcomp (-3 + k[j]))*Pi*(1 - t[i])) +pow(E,-((sqrt(0.625 - sqrt(5)/8) - comp6*(1 + sqrt(5)))* (dcomp (-3 + k[j]))*Pi*t[i])))- sqrt(2)*cos((-3 + k[j])*Pi*t[i]);
            M5[j] =ans.real();
        }
        M5[j] = M5[j]/sqrt(nn);
    }
    return M5;
}

int main(int argc, char *argv[]) {
    constexpr int nn = 250;
    std::vector<double> M = compute(nn);
    for (int i=240; i<250; i++) std::cout << M[i] << " ";
    std::cout << std::endl;
    exit(0);
}

[coatless]

@eddelbuettel briefly glancing at the code:

arma::mat would be a double-double; but, dcomp is typedef complex<double>. The inside use of ambiguous sqrt()/pow() without casting might be creating issues.

Something like:

static_cast<double>(int);

std::pow(2.0, - 1.0);

I'm thinking this way as I think Armadillo 15.0.0 added FP16 support and the jump over to C++14 as min standard.

[eddelbuettel]

Yes -- switching to std::pow() etc was the first thing I tried when I ran into this two or so weeks ago. It is not it.

arma::mat is used as a matrix is filled, but we only need the first row to demonstrate that its last ten elements become nan so the example is valid. See the plain C++ example in the second comment. The complex<double> is also not the cause---if I glance some more at the code example you will see "everything" happens at a tuple <i,j> level. That should work. As it did with Armadillo up to and including 14.6.3 (but alas not under the Intel compiler).

FP16 has nothing to do with this, nor has the C++14 standard AFAICT. Maybe it is single element access?

[eddelbuettel]

Ok, truly minimal example is taking the last C++ 'no arma' example and just including #include <armadillo> (with an appropriate -I to point to the 15.0.0 release) and bam it blows up even when all the data structures are standard C++ and not Armadillo.

edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ g++ -I/usr/local/lib/R/site-library/RcppArmadillo/include/current  -o issue474_ebsc_rank_main_min_arma issue474_ebsc_rank_main_min_arma.cpp -lm
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ ./issue474_ebsc_rank_main_min_arma 
-nan -nan -nan -nan -nan -nan -nan -nan -nan -nan 
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ 

Code below

// The eBsc package on CRAN (at version 4.17) reports an error under the Intel compiler.
// This error now also appears with RcppArmadillo 15.0.0. In the example for plot.eBsc,
// the function drbasis() is called with n=250 and a second parameter p that varies from
// 1 to 6. In the case of 5, NA values return leading an rank failure. The drbasis()
// function calls an underlying C function drbasis from file drbasisC.cpp which has case
// statements for the differenct values of p.  We extrace the one for 5 here.

#include <armadillo>

std::vector<double> compute(int nn) {
    using namespace std;
    const double E = exp(1);
    vector<double> t(nn);
    vector<double> k(nn);
    typedef complex<double> dcomp;
    double Pi = M_PI;
    complex<double> comp0 (0,1);
    complex<double> comp1 (-1,0);
    complex<double> comp2 (0,-0.5);
    complex<double> comp3 (1,-1);
    complex<double> comp4 (-1,1);
    complex<double> comp5 (0,0.5);
    complex<double> comp6 (0,0.25);
    int const0 = 1;
    double const1 = -sqrt(3);
    double const2 = 2*sqrt(3);
    double const3 = -sqrt(5);
    double const4 = 6*sqrt(5);
    double const5 = -6*sqrt(5);
    double const6 = -sqrt(7);
    double const7 = 12*sqrt(7);
    double const8 = -30*sqrt(7);
    double const9 = 20*sqrt(7);
    double const10 = 3;
    double const11 = -60;
    double const12 = 270;
    double const13 = -420;
    double const14 = 210;
    /*
    double const15 = -sqrt(11);
    double const16 = 30*sqrt(11);
    double const17 = -210*sqrt(11);
    double const18 = 560*sqrt(11);
    double const19 = 630*sqrt(11);
    double const20 = 252*sqrt(11);
    */

    for (int j=0;j<nn;j++) {
        t[j]= j/double((nn-1));
    }
    for (int j=0;j<nn;j++) {
        k[j]= j+1;
    }

    vector<double> M5(nn);
    constexpr int i = 0;
    for (int j=0;j<nn;j++) {
        if (j==0) {
            M5[j] = const0;
        } else if (j==1) {
            M5[j] = const1 + const2*t[i];
        } else if (j==2) {
            M5[j] = const3 + const4*t[i] + const5* pow(t[i],2);
        } else if (j==3) {
            M5[j] = const6 + const7*t[i] + const8* pow(t[i],2) + const9* pow(t[i],3);
        } else if (j==4) {
            M5[j] = const10 + const11*t[i] + const12* pow(t[i],2) + const13* pow(t[i],3) + const14* pow(t[i],4);
        } else {
            dcomp ans = (sqrt(2)*(1 + sqrt(5)/2) -
                         (comp5*(sqrt(10 - 2*sqrt(5)) +
                                 sqrt(2*(5 + sqrt(5)))))/sqrt(2))*(pow(comp1,1 + k[j])/pow(E,pow(comp1,0.1)* (dcomp(-3 + k[j]))*Pi*(1 - t[i])) + pow(E,-(pow(comp1,0.1)*(dcomp(-3 + k[j]))* Pi*t[i]))) +
                (-(1/sqrt(2)) -
                 (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/
                 sqrt(2))*(pow(comp1,1 + k[j])/pow(E,pow(comp1,0.3)*(dcomp (-3 + k[j]))*Pi*(1 - t[i])) + pow(E,-(pow(comp1,0.3)*(dcomp (-3 + k[j]))*Pi*t[i]))) +
                (sqrt(2)*(1 + sqrt(5)/2) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
                *(pow(comp1,1 + k[j])/ pow(E,(sqrt(0.625 + sqrt(5)/8) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k[j]))*Pi*
                                           (1 - t[i])) +pow(E,-((sqrt(0.625 + sqrt(5)/8.) - comp6*(-1 + sqrt(5)))*(dcomp (-3 + k[j]))*
                                                                Pi*t[i]))) +(-(1/sqrt(2)) + (comp5*(sqrt(10 - 2*sqrt(5)) + sqrt(2*(5 + sqrt(5)))))/sqrt(2))
                *(pow(comp1,1 + k[j])/pow(E,(dcomp(sqrt(0.625 - sqrt(5)/8)) - comp6*(1 + sqrt(5)))*
                                          (dcomp (-3 + k[j]))*Pi*(1 - t[i])) +pow(E,-((sqrt(0.625 - sqrt(5)/8) - comp6*(1 + sqrt(5)))* (dcomp (-3 + k[j]))*Pi*t[i])))- sqrt(2)*cos((-3 + k[j])*Pi*t[i]);
            M5[j] =ans.real();
        }
        M5[j] = M5[j]/sqrt(nn);
    }
    return M5;
}

int main(int argc, char *argv[]) {
    constexpr int nn = 250;
    std::vector<double> M = compute(nn);
    for (int i=240; i<250; i++) std::cout << M[i] << " ";
    std::cout << std::endl;
    exit(0);
}

[eddelbuettel]

And the truly minimal example of using only the header persists even when we clean up the namespace use and prefix everything from the standard library with std::.

edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ g++ -I/usr/local/lib/R/site-library/RcppArmadillo/include/current  -o issue474_ebsc_rank_main_min_arma_std issue474_ebsc_rank_main_min_arma_std.cpp -lm
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ ./issue474_ebsc_rank_main_min_arma_std 
-nan -nan -nan -nan -nan -nan -nan -nan -nan -nan 
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ 

Code below.

// The eBsc package on CRAN (at version 4.17) reports an error under the Intel compiler.
// This error now also appears with RcppArmadillo 15.0.0. In the example for plot.eBsc,
// the function drbasis() is called with n=250 and a second parameter p that varies from
// 1 to 6. In the case of 5, NA values return leading an rank failure. The drbasis()
// function calls an underlying C function drbasis from file drbasisC.cpp which has case
// statements for the differenct values of p.  We extrace the one for 5 here.

#if defined(USE_ARMA)
#include <armadillo>
#else
#include <cmath>
#include <complex>
#include <iostream>
#include <vector>
#endif

std::vector<double> compute(int nn) {
    const double E = std::exp(1);
    std::vector<double> t(nn);
    std::vector<double> k(nn);
    typedef std::complex<double> dcomp;
    double Pi = M_PI;
    std::complex<double> comp0 (0,1);
    std::complex<double> comp1 (-1,0);
    std::complex<double> comp2 (0,-0.5);
    std::complex<double> comp3 (1,-1);
    std::complex<double> comp4 (-1,1);
    std::complex<double> comp5 (0,0.5);
    std::complex<double> comp6 (0,0.25);
    int const0 = 1;
    double const1 = -std::sqrt(3);
    double const2 = 2*std::sqrt(3);
    double const3 = -std::sqrt(5);
    double const4 = 6*std::sqrt(5);
    double const5 = -6*std::sqrt(5);
    double const6 = -std::sqrt(7);
    double const7 = 12*std::sqrt(7);
    double const8 = -30*std::sqrt(7);
    double const9 = 20*std::sqrt(7);
    double const10 = 3;
    double const11 = -60;
    double const12 = 270;
    double const13 = -420;
    double const14 = 210;
    /*
    double const15 = -std::sqrt(11);
    double const16 = 30*std::sqrt(11);
    double const17 = -210*std::sqrt(11);
    double const18 = 560*std::sqrt(11);
    double const19 = 630*std::sqrt(11);
    double const20 = 252*std::sqrt(11);
    */

    for (int j=0;j<nn;j++) {
        t[j]= j/double((nn-1));
    }
    for (int j=0;j<nn;j++) {
        k[j]= j+1;
    }

    std::vector<double> M5(nn);
    constexpr int i = 0;
    for (int j=0;j<nn;j++) {
        if (j==0) {
            M5[j] = const0;
        } else if (j==1) {
            M5[j] = const1 + const2*t[i];
        } else if (j==2) {
            M5[j] = const3 + const4*t[i] + const5* std::pow(t[i],2);
        } else if (j==3) {
            M5[j] = const6 + const7*t[i] + const8* std::pow(t[i],2) + const9* std::pow(t[i],3);
        } else if (j==4) {
            M5[j] = const10 + const11*t[i] + const12* std::pow(t[i],2) + const13* std::pow(t[i],3) + const14* std::pow(t[i],4);
        } else {
            dcomp ans = (std::sqrt(2)*(1 + std::sqrt(5)/2) -
                         (comp5*(std::sqrt(10 - 2*std::sqrt(5)) +
                                 std::sqrt(2*(5 + std::sqrt(5)))))/std::sqrt(2))*(std::pow(comp1,1 + k[j])/std::pow(E,std::pow(comp1,0.1)* (dcomp(-3 + k[j]))*Pi*(1 - t[i])) + std::pow(E,-(std::pow(comp1,0.1)*(dcomp(-3 + k[j]))* Pi*t[i]))) +
                (-(1/std::sqrt(2)) -
                 (comp5*(std::sqrt(10 - 2*std::sqrt(5)) + std::sqrt(2*(5 + std::sqrt(5)))))/
                 std::sqrt(2))*(std::pow(comp1,1 + k[j])/std::pow(E,std::pow(comp1,0.3)*(dcomp (-3 + k[j]))*Pi*(1 - t[i])) + std::pow(E,-(std::pow(comp1,0.3)*(dcomp (-3 + k[j]))*Pi*t[i]))) +
                (std::sqrt(2)*(1 + std::sqrt(5)/2) + (comp5*(std::sqrt(10 - 2*std::sqrt(5)) + std::sqrt(2*(5 + std::sqrt(5)))))/std::sqrt(2))
                *(std::pow(comp1,1 + k[j])/ std::pow(E,(std::sqrt(0.625 + std::sqrt(5)/8) - comp6*(-1 + std::sqrt(5)))*(dcomp (-3 + k[j]))*Pi*
                                           (1 - t[i])) +std::pow(E,-((std::sqrt(0.625 + std::sqrt(5)/8.) - comp6*(-1 + std::sqrt(5)))*(dcomp (-3 + k[j]))*
                                                                Pi*t[i]))) +(-(1/std::sqrt(2)) + (comp5*(std::sqrt(10 - 2*std::sqrt(5)) + std::sqrt(2*(5 + std::sqrt(5)))))/std::sqrt(2))
                *(std::pow(comp1,1 + k[j])/std::pow(E,(dcomp(std::sqrt(0.625 - std::sqrt(5)/8)) - comp6*(1 + std::sqrt(5)))*
                                          (dcomp (-3 + k[j]))*Pi*(1 - t[i])) +std::pow(E,-((std::sqrt(0.625 - std::sqrt(5)/8) - comp6*(1 + std::sqrt(5)))* (dcomp (-3 + k[j]))*Pi*t[i])))- std::sqrt(2)*cos((-3 + k[j])*Pi*t[i]);
            M5[j] =ans.real();
        }
        M5[j] = M5[j]/std::sqrt(nn);
    }
    return M5;
}

int main(int argc, char *argv[]) {
    constexpr int nn = 250;
    std::vector<double> M = compute(nn);
    for (int i=240; i<250; i++) std::cout << M[i] << " ";
    std::cout << std::endl;
    exit(0);
}

PS Made a second minimal change and allow to fall back to not using Armadillo via a #define. If we do, it passes:

edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ g++ -o issue474_ebsc_rank_main_min_arma_std issue474_ebsc_rank_main_min_arma_std.cpp -lm
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ ./issue474_ebsc_rank_main_min_arma_std 
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ 

but as above (now plus -D....)

edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ g++ -DUSE_ARMA -I/usr/local/lib/R/site-library/RcppArmadillo/include/current  -o issue474_ebsc_rank_main_min_arma_std issue474_ebsc_rank_main_min_arma_std.cpp -lm
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ ./issue474_ebsc_rank_main_min_arma_std 
-nan -nan -nan -nan -nan -nan -nan -nan -nan -nan 
edd@paul:~/git/rcpparmadillo/local(feature/arma_15.0)$ 

[eddelbuettel]

Now also reported over at the Armadillo repo in issue #283 there.

[eddelbuettel]

The files used here have been added to the local/ directory of the main branch,

[eddelbuettel]

This turns out to have been a genuine bug in Armadillo fixed by 15.0.1.