numericDeriv alters result of eval in R?

Question

I've run into something unexpected (for me, maybe). Take a look at the code bellow:

> env1 = new.env()
> env1$x = rnorm(10)
> env1$a = 3
> eval(quote(x^a), env1)
 [1]   0.26046508   1.37188255  -0.17077313   0.25959068  -3.38565765
 [6]   0.11537940  -0.06051991  -1.13542624  16.84291641 -10.60959041
> numericDeriv(quote(x^a), "a", env1)
Error in numericDeriv(quote(x^a), "a", env1) : 
  Missing value or an infinity produced when evaluating the model
> eval(quote(x^a), env1)
 [1]  0.2604651  1.3718826        NaN  0.2595907        NaN  0.1153794
 [7]        NaN        NaN 16.8429171        NaN
> env1$a
[1] 3
> eval(quote(x^a), env1)
 [1]  0.2604651  1.3718826        NaN  0.2595907        NaN  0.1153794
 [7]        NaN        NaN 16.8429171        NaN
> env1$a = 3
> eval(quote(x^a), env1)
 [1]   0.26046508   1.37188255  -0.17077313   0.25959068  -3.38565765
 [6]   0.11537940  -0.06051991  -1.13542624  16.84291641 -10.60959041

The error in numericDeriv is fine since derivative of x^a wrt a is x^a * log(x) and log is not defined for negative numbers. However, why does it changes the result of eval(quote(x^a), env1)?

Edit 1

As pointed by @Roland, the env1$a changes from 3 to 3.00000004470348. I've also checked that this doesn't happen when numericDeriv completes without errors

Edit 2: Session Info

As requested by @jay.sf, here is the output of sessionInfo():

> sessionInfo()
R version 4.0.1 (2020-06-06)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods  
[8] base     

loaded via a namespace (and not attached):
[1] compiler_4.0.1

Edit 3

I'm not sure if the numericDeriv function still uses this C code but if so, at the very bottom of the page we find

for(i = 0, start = 0; i < LENGTH(theta); i++) {
  for(j = 0; j < LENGTH(VECTOR_ELT(pars, i)); j++, start += LENGTH(ans)) {
    SEXP ans_del;
    double origPar, xx, delta;

    origPar = REAL(VECTOR_ELT(pars, i))[j];
    xx = fabs(origPar);
    delta = (xx == 0) ? eps : xx*eps;
    REAL(VECTOR_ELT(pars, i))[j] += rDir[i] * delta;
    PROTECT(ans_del = eval(expr, rho));
    if(!isReal(ans_del)) ans_del = coerceVector(ans_del, REALSXP);
    UNPROTECT(1);
    for(k = 0; k < LENGTH(ans); k++) {
      if (!R_FINITE(REAL(ans_del)[k]))
        error(_("Missing value or an infinity produced when evaluating the model"));
      REAL(gradient)[start + k] = rDir[i] * (REAL(ans_del)[k] - REAL(ans)[k])/delta;
    }
    REAL(VECTOR_ELT(pars, i))[j] = origPar;
  }
}

and we can see that this behavior, apparently, comes from the inner most for loop. Maybe one quick fix would be changing the if in the inner most to

if (!R_FINITE(REAL(ans_del)[k])) {
  REAL(VECTOR_ELT(pars, i))[j] = origPar;
  error(_("Missing value or an infinity produced when evaluating the model"));
}

I have reported to the author listed on ?numericDeriv.