0

I am currently implementing Clenshaw's algorithm with Rcpp to speed up my previous implementation in R. My current implementation is as follows (note that I am using RcppParallel for other functions defined in the same source file; RcppParalell is not used in this specific function, but I've left the headers in case this is somehow relevant):

#include <Rcpp.h>
#include <RcppParallel.h>
using namespace Rcpp;
using namespace RcppParallel;

// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
NumericVector clenshawAllDerivatives(double t, int N, double Ta, double Tb, NumericVector Coeffs, int derivativesOrder) {
    double tau = (2*t-Ta-Tb)/(Tb-Ta);
    double helperValues[derivativesOrder + 1][3];
    double scale;
    for(double i = N; i > 1; i--) {
        helperValues[0][2] = helperValues[0][1];
        helperValues[0][1] = helperValues[0][0];
        helperValues[0][0] = 2*tau*helperValues[0][1]-helperValues[0][2] + Coeffs[i - 1];
        scale=2.0;
        for(int j = 1; j <= derivativesOrder; j++) {
            helperValues[j][2] = helperValues[j][1];
            helperValues[j][1] = helperValues[j][0];
            helperValues[j][0] = scale*helperValues[j-1][1] + 2*tau*helperValues[j][1] - helperValues[j][2];
            scale += 2.0;
        }
    }
    NumericVector output(derivativesOrder + 1);
    output[0] = tau*helperValues[0][0] - helperValues[0][1] + Coeffs[0];
    scale = 1.0;
    double scale2initial = ((Tb-Ta)/2 * 86400.0), scale2 = scale2initial;
    for(int j = 1; j <= derivativesOrder; j++) {
        output[j] = (scale*helperValues[j-1][0] + tau*helperValues[j][0] - helperValues[j][1]) / scale2;
        scale += 1.0;
        scale2 = scale2 * scale2initial;
    }
    return output;
}

An example of application of the function, with example input values:

clenshawAllDerivatives(59568.5, 11, 59568, 59584, c(-1.281626e+06, -4.069960e+03, 2.725817e+01, -9.715712e-02, -1.115373e-03, -5.121949e-04, -9.068147e-05, -6.829206e-06, 1.657523e-07 , 1.406006e-07, 2.273966e-08), 1)

When run multiple times, this returns most often the expected correct output of c(-1.277790e+06, -6.037188e-03). However, sometimes it returns instead wrong values, typically very high numbers.

Any help to identify the cause of this unexpected behavior would be greatly appreciated!

Rafa
  • 157
  • 6
  • 1
    `helperValues[0][2] = helperValues[0][1];` is using `helperValues[0][1]` uninitialized – 463035818_is_not_an_ai Sep 27 '22 at 08:42
  • Thanks a lot for pointing this out! I thought declaring the `helperValues` array with dimensions `[derivativesOrder + 1][3]` would initialize it to all 0s by default? Is this not the case? – Rafa Sep 27 '22 at 08:43
  • 1
    No. And `double helperValues[derivativesOrder + 1][3];` is not standard C++. [Why aren't variable-length arrays part of the C++ standard?](https://stackoverflow.com/questions/1887097/why-arent-variable-length-arrays-part-of-the-c-standard) – 463035818_is_not_an_ai Sep 27 '22 at 08:44
  • 1
    Not the issue, but use `R_xlen_t` as type for indices (of Rcpp vectors). Also, it is good practice to define constants as `const`. – Roland Sep 27 '22 at 08:47
  • Thanks a lot, I understand now why I was getting such strange behavior. I am requiring C++11, it seems it compiles fine as long as I don't initialize the variable-length array and just declare it. I have gotten around this by having `memset(helperValues, 0, sizeof helperValues);` just after array declaration. Is there a cleaner way to do this, or is it good practice in that way? – Rafa Sep 27 '22 at 09:32

0 Answers0