New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 1693 - pexp is inaccurate outside the clamping area
pexp is inaccurate outside the clamping area
 Status: CONFIRMED None Eigen Unclassified Core - vectorization (show other bugs) 3.4 (development) All All Low Accuracy Problem Nobody 1687 Show dependency tree / graph

 Reported: 2019-03-18 17:39 UTC by Christoph Hertzberg 2019-03-18 23:45 UTC (History) 4 users (show) chtz gael.guennebaud jacob.benoit.1 markos

Attachments

 Christoph Hertzberg 2019-03-18 17:39:47 UTC ```Values of large magnitude are clamped by pexp and result in values which should only occur outside that range: Array4f x(4); x << -88.37, -88.38, 88.37, 88.38; Array4f y = exp(x); Array4f z; 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"; Output: 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.``` Gael Guennebaud 2019-03-18 20:59:22 UTC ```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.``` Christoph Hertzberg 2019-03-18 23:45:37 UTC ```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 in detail, but I assume it has essentially the same problems.```

 Note You need to log in before you can comment on or make changes to this bug.