avoid some index arithmetic

[richfitz]

p_G_D[, , ] <- if (as.integer(step) >= n_p_G_D_steps)
  p_G_D_step[n_p_G_D_steps, i] * rel_p_G_D[i, j, k] else
    p_G_D_step[step + 1, i] * rel_p_G_D[i, j, k]

We currently generate code that does a lot of work to create the index into p_G_D and from rel_p_G_D

    for (int i = 1; i <= shared->dim_p_G_D_1; ++i) {
      for (int j = 1; j <= shared->dim_p_G_D_2; ++j) {
        for (int k = 1; k <= shared->dim_p_G_D_3; ++k) {
          internal.p_G_D[i - 1 + shared->dim_p_G_D_1 * (j - 1) + shared->dim_p_G_D_12 * (k - 1)] = (static_cast<int>(step) >= shared->n_p_G_D_steps ? shared->p_G_D_step[shared->dim_p_G_D_step_1 * (i - 1) + shared->n_p_G_D_steps - 1] * shared->rel_p_G_D[shared->dim_rel_p_G_D_12 * (k - 1) + shared->dim_rel_p_G_D_1 * (j - 1) + i - 1] : shared->p_G_D_step[shared->dim_p_G_D_step_1 * (i - 1) + step + 1 - 1] * shared->rel_p_G_D[shared->dim_rel_p_G_D_12 * (k - 1) + shared->dim_rel_p_G_D_1 * (j - 1) + i - 1]);
        }
      }
    }

This might be simpifiable to:

    for (int i = 1, idx = 0; i <= shared->dim_p_G_D_1; ++i) {
      for (int j = 1; j <= shared->dim_p_G_D_2; ++j) {
        for (int k = 1; k <= shared->dim_p_G_D_3; ++k, ++idx) {
          internal.p_G_D[idx] = (static_cast<int>(step) >= shared->n_p_G_D_steps ? shared->p_G_D_step[shared->dim_p_G_D_step_1 * (i - 1) + shared->n_p_G_D_steps - 1] * shared->rel_p_G_D[idx] : shared->p_G_D_step[shared->dim_p_G_D_step_1 * (i - 1) + step + 1 - 1] * shared->rel_p_G_D[idx]);
        }
      }
    }

[richfitz]

Looking at this, I think that loops run the wrong way around to support this, we should be doing k, j, i not i, j, k so that i moves fastest - this might be causing some access issues at the moment actually.