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.