Values of large magnitude are clamped by pexp and result in values which should only occur outside that range:
x << -88.37, -88.38, 88.37, 88.38;
Array4f y = exp(x);
for(int i=0; i<4; ++i) z(i) = std::exp(x(i));
std::cout << "x = " << x.transpose() << "\nexp(x) = " << y.transpose() << "\nshould be " << z.transpose() << "\n";
x = -88.37 -88.38 88.37 88.38
exp(x) = 0 0 2.39114e+38 inf
should be 4.18211e-39 4.14052e-39 2.39114e+38 2.41516e+38
Not super-critical, and likely irrelevant in most cases, but getting different results in the scalar and packet implementation is not good.
4.18211e-39 4.14052e-39 are denormals, so not sure we should care at all.
for 88.38, maybe the upper bound can be very slightly increased, I did not checked.
As far as I understand the current implementation, the upper bound is chosen such that
round(x/log(2)) == 128,
i.e., when adding the exponential bias of 127, the result is infinity.
This could be handled either by doing
m = floor(x/log(2));
y = x-m*log(2); // y in [0,log(2)] --> exp(y) in [1,2]
or if m=round(x/log(2)), by directly adding m to the exponent of exp(y) (this would require that the upper bound is chosen such that exp(upper) == 2^128, otherwise NaNs would be generated. The first alternative would require a different polynomial for approximation.
Both solutions will likely shorten the valid range at the lower bound.
An alternative would be to split m, and multiply twice by 2^(m/2), or by calculating (exp(x/2))^2 (I assume that will slightly decrease the accuracy, though).
I did not yet check pexp<double> in detail, but I assume it has essentially the same problems.
-- GitLab Migration Automatic Message --
This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.
You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1693.