New user self-registration is disabled due to spam. Please email eigen-core-team @ 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
Summary: pexp is inaccurate outside the clamping area
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - vectorization (show other bugs)
Version: 3.4 (development)
Hardware: All All
: Low Accuracy Problem
Assignee: Nobody
Depends on: 1687
  Show dependency treegraph
Reported: 2019-03-18 17:39 UTC by Christoph Hertzberg
Modified: 2019-03-18 23:45 UTC (History)
4 users (show)


Description 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";

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.
Comment 1 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.
Comment 2 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<double> 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.