# When algorithm is placed inside loop it produces different results, C++

I create the following algorithm in Rcpp and compile it in R.

```// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]

arma::colvec Demo(arma::mat n, int K){

arma::colvec N(K);

for(int j=0; j<K; ++j){
for(int i=0; i<(K-j); ++i){
N[j] += accu(n.submat(i,0,i,j));
}
}
return N;
}

/***R
K = 4
n = cbind(c(1008, 5112, 1026, 25, 0), 0, 0, 0, 0)
Demo(n,K)

for(i in 1:3){
print(Demo(n,K))
print(K)
print(n)
}
*/
```

However, something really weird happens when I run it inside a loop.

For example, if I have

```> K = 4
> n
[,1] [,2] [,3] [,4] [,5]
[1,] 1008    0    0    0    0
[2,] 5112    0    0    0    0
[3,] 1026    0    0    0    0
[4,]   25    0    0    0    0
[5,]    0    0    0    0    0
```

Then if I run the algorithm `Demo` a single time I receive the correct result

```> Demo(n,K)
[,1]
[1,] 7171
[2,] 7146
[3,] 6120
[4,] 1008
```

However, if I run it multiple times inside a loop, it starts to behave weird

```for(i in 1:3){
print(Demo(n,K))
print(K)
print(n)
}
[,1]
[1,] 7171
[2,] 7146
[3,] 6120
[4,] 1008
[1] 4
[,1] [,2] [,3] [,4] [,5]
[1,] 1008    0    0    0    0
[2,] 5112    0    0    0    0
[3,] 1026    0    0    0    0
[4,]   25    0    0    0    0
[5,]    0    0    0    0    0
[,1]
[1,] 14342
[2,] 14292
[3,] 12240
[4,]  2016
[1] 4
[,1] [,2] [,3] [,4] [,5]
[1,] 1008    0    0    0    0
[2,] 5112    0    0    0    0
[3,] 1026    0    0    0    0
[4,]   25    0    0    0    0
[5,]    0    0    0    0    0
[,1]
[1,] 21513
[2,] 21438
[3,] 18360
[4,]  3024
[1] 4
[,1] [,2] [,3] [,4] [,5]
[1,] 1008    0    0    0    0
[2,] 5112    0    0    0    0
[3,] 1026    0    0    0    0
[4,]   25    0    0    0    0
[5,]    0    0    0    0    0
```

In the first run, it computes it correctly, then in the second run it gives the correct output multiplied by 2, and in the third run, it gives the correct output multiplied by 3. But based on the algorithm steps, I do not see an obvious step that produces this kind of behavior.

The correct output should have been

```for(i in 1:3){
print(Demo(n,K))
}
[,1]
[1,] 7171
[2,] 7146
[3,] 6120
[4,] 1008
[,1]
[1,] 7171
[2,] 7146
[3,] 6120
[4,] 1008
[,1]
[1,] 7171
[2,] 7146
[3,] 6120
[4,] 1008
```

You are incrementing `N` in place via `+=`.

Your function fails to ensure it is initialized at zero. `Rcpp` tends to do that by default (as I think it is prudent) — but this can be suppressed for speed if you know you are doing.

A minimally repaired version of your code (with the correct header, and a call to `.fill(0)`) follows.

```// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::colvec Demo(arma::mat n, int K){
arma::colvec N(K);
N.fill(0);   // important, or construct as N(k, arma::fill::zeros)
for(int j=0; j<K; ++j){
for(int i=0; i<(K-j); ++i){
N[j] += accu(n.submat(i,0,i,j));
}
}
return N;
}

/***R
K = 4
n = cbind(c(1008, 5112, 1026, 25, 0), 0, 0, 0, 0)
Demo(n,K)

for(i in 1:3) {
print(Demo(n,K))
print(K)
print(n)
}
*/
```

You could also call `.zeros()` (once constructed) or use `zeros(k)` (to construct) or … deploy a number of different ways to ensure your content is cleared before adding to it.

The shortest, after checking the documentation, may be `arma::colvec(N, arma::fill::zeros)`.