This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 1693 - pexp is inaccurate outside the clamping area
Summary: pexp is inaccurate outside the clamping area
Status: CONFIRMED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - vectorization (show other bugs)
Version: 3.4 (development)
Hardware: All All
: Low Accuracy Problem
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on: 1687
Blocks:
  Show dependency treegraph
 
Reported: 2019-03-18 17:39 UTC by Christoph Hertzberg
Modified: 2019-12-04 18:33 UTC (History)
4 users (show)



Attachments

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

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.
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.
Comment 3 Nobody 2019-12-04 18:33:51 UTC
-- 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.

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