Same Rcpp function returns different output if print statement added

Question

A C++ function I wrote using Rcpp gives different output depending on whether or not I have a Rcout or Rprintf statement in the code. The code below returns 1, which is the correct value, for the function with a print statement H_sigma_1(). However, for H_sigma_2(), the function without a print statement, the function returns 2. I've tested this on Ubuntu 16.04.1 as well as CentOS 6.8. Though, I was not able to reproduce this bug on Windows 10. Therefore, this seems to be a Linux issue.

Rcpp Code

library(Rcpp)

cppFunction ( 
  "double H_sigma_1(IntegerVector sigma, NumericMatrix J, NumericVector h)
  {
    double first_sum, second_sum = 0;
    int n = sigma.size();

    for(int i = 0; i < n; i++)
    {
      for(int j = 0; j < n; j++)
      {
        // skip inside loop if i >= j to stop double counting
        if(i >= j) {continue;}
        first_sum += J(i, j) * sigma[i] * sigma[j];
        Rcout << first_sum << std::endl;
      }
      second_sum += h[i] * sigma[i];
    }
    return(-1.0 * first_sum - second_sum);
  }"
)

cppFunction (
  "double H_sigma_2(IntegerVector sigma, NumericMatrix J, NumericVector h)
  {
    double first_sum, second_sum = 0;
    int n = sigma.size();

    for(int i = 0; i < n; i++)
    {
      for(int j = 0; j < n; j++)
      {
        // skip inside loop if i >= j to stop double counting
        if(i >= j) {continue;}
        first_sum += J(i, j) * sigma[i] * sigma[j];
        // Rcout << first_sum << std::endl;
      }
      second_sum += h[i] * sigma[i];
    }
    return(-1.0 * first_sum - second_sum);
  }"
)

Example call

Setup:

n = 2
params = rep(1, n) 
h = rep(params[1], n)
J = toeplitz(c(0, params[2], rep(0, n - 2)))

Test 1:

H_sigma_1(c(-1, -1), J, h)

Output:

1
[1] 1

Test 2

H_sigma_2(c(-1, -1), J, h)

Output

[1] 2

Answer 1

The problem that you ran into is not explicitly declaring starting values for sums before using them.

double first_sum, second_sum = 0;

is not the same as:

double first_sum = 0, second_sum = 0;

See compiler warning flags:

file11d36f564b15.cpp:17:5: warning: variable 'first_sum' is uninitialized when used here [-Wuninitialized]
    first_sum += J(i, j) * sigma[i] * sigma[j];
    ^~~~~~~~~
file11d36f564b15.cpp:8:21: note: initialize the variable 'first_sum' to silence this warning
    double first_sum, second_sum = 0;

---

You immediately use:

first_sum += J(i, j) * sigma[i] * sigma[j];

without setting a first_sum = ...;

---

Furthermore, another issue is:

second_sum = 0;

initializes a double with an integer value. While this issue is minor in scope, to correct this, all one has to do is use 0.0 instead of 0.

second_sum = 0.0;

This also applies for first_sum.

---

Code with above fixes:

library(Rcpp)

cppFunction ( 
    "double H_sigma_1(IntegerVector sigma, NumericMatrix J, NumericVector h)
    {
    double first_sum = 0.0, second_sum = 0.0;
    int n = sigma.size();

    for(int i = 0; i < n; i++) {
      for(int j = 0; j < n; j++) {
        // skip inside loop if i >= j to stop double counting
        if(i >= j) {continue;}
        first_sum += J(i, j) * sigma[i] * sigma[j];
        Rcout << first_sum << std::endl;
      }
      second_sum += h[i] * sigma[i];
    }
    return(-1.0 * first_sum - second_sum);
    }"
)

cppFunction ( 
    "double H_sigma_2(IntegerVector sigma, NumericMatrix J, NumericVector h)
    {
    double first_sum = 0.0, second_sum = 0.0;
    int n = sigma.size();

    for(int i = 0; i < n; i++) {
      for(int j = 0; j < n; j++) {
        // skip inside loop if i >= j to stop double counting
        if(i >= j) {continue;}
        first_sum += J(i, j) * sigma[i] * sigma[j];
      }
      second_sum += h[i] * sigma[i];
    }
    return(-1.0 * first_sum - second_sum);
    }"
)

Test:

n = 2
params = rep(1, n) 
h = rep(params[1], n)
J = toeplitz(c(0, params[2], rep(0, n - 2)))

H_sigma_1(c(-1, -1), J, h)

Output:

1
[1] 1

Test 2:

H_sigma_2(c(-1, -1), J, h)

Output:

[1] 1