for (j = i + 1; j <= n; j++) {
//Rprintf("S[%i]=%f, S[%i]=%f, D[%i,%i]=%f, B=%f",i,S[i],j,S[j],i,j,D[give_index(i,j,n)],B);
A=S[i]+S[j]-B*D[give_index(i,j,n)];
for (j = i + 1; j <= n; j++) {
//Rprintf("S[%i]=%f, S[%i]=%f, D[%i,%i]=%f, B=%f",i,S[i],j,S[j],i,j,D[give_index(i,j,n)],B);
A=S[i]+S[j]-B*D[give_index(i,j,n)];