Skip to content

Commit

Permalink
The correction vector is now only used when n > 1, so only declare an…
Browse files Browse the repository at this point in the history
…d initialize it inside the n > 1 block.
  • Loading branch information
aadler committed Aug 23, 2023
1 parent bbb876e commit e75c412
Showing 1 changed file with 4 additions and 6 deletions.
10 changes: 4 additions & 6 deletions src/compHorner.c
Original file line number Diff line number Diff line change
Expand Up @@ -47,17 +47,13 @@ extern SEXP compHorner_c(SEXP x, SEXP a) {
SEXP ret = PROTECT(allocVector(REALSXP, m));
double *pret = REAL(ret);

// Correction vector for compensated Horner
double correction[m];

// Since we will return 0 if n = 0---which is technically legal if asking for
// the best constant estimate or using rational to test polynomial---we must
// initialize the ret vector and the correction vector to 0.
// initialize the ret vector to 0.
memset(pret, 0, m * sizeof(double));
memset(correction, 0, m * sizeof(double));

// If n is at least 1, initialize ret with the last value of a. If n == 1 then
// there is no need to calculate piM or sigM and correction is already 0.
// there is no need to calculate piM, sigM, or correction.
if (n > 0) {
for (int i = 0; i < m; ++i) {
pret[i] = pa[nm1];
Expand Down Expand Up @@ -87,8 +83,10 @@ extern SEXP compHorner_c(SEXP x, SEXP a) {
if (n > 1) {
// x element of twoProdFMA used more than once so define as variable
volatile double Ax;
double correction[m];
double piM[nm1][m];
double sigM[nm1][m];
memset(correction, 0, m * sizeof(double));
memset(piM, 0, m * nm1 * sizeof(double));
memset(sigM, 0, m * nm1 * sizeof(double));

Expand Down

0 comments on commit e75c412

Please sign in to comment.