-
Notifications
You must be signed in to change notification settings - Fork 33
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
replacing expm1
implementation
#131
Comments
good catch |
here is my quick fix ... it should help you while I review it
|
#135 has a much better fix. |
#135 has a much better fix. |
Shouldn't this remain open until #135 is merged? |
yes |
expm1
expm1
for tiny values
I do not see the described behavior
I am closing this -- if there is something I missed, let me know |
this is still an issue for expm1(.5) |
???
|
Sorry, I chose the wrong number.
so in this case, |
After spending some time with the implementation, it has become clear that the only fix is to replace I did generate factorial coeffs for Double64 in the event we want to run the series solution another time. 1! .. 36! will fit in a Double64 without error. Consider 42!
For the record, these are appropriate for use with Double64 smallerfactorials = [ largerfactorials = [ Double64 inv factorials are exactly recoverable through factorial(big(30))smaller_invfactorials = [ these below are approximations, above are exactly recoverable(31, Double64(1.216125041553518e-34, 5.586290567888806e-51)) invfactorials = [ |
expm1
for tiny valuesexpm1
implementation
I implemented something here which seems to work: Might be worth taking a look at .... |
That looks like a pretty good version, although it won't be the most accurate since it ends in 9 multiplications that are uncompensated, so I think the error will be around 10 ULP error. |
@oscardssmith I took the Anyway, is there any easy way to fix that uncompensated multiplication problem? |
You could compensate the multiplication, but it would likely be a noticeable slowdown. 10 ULP may just be the price here. (The other option would be a table-based method) |
Here's a C++ implementation of https://github.com/tuwien-cms/libxprec/blob/v0.4.2/src/exp.cxx The idea is similar to what you are doing already ... but we simply split off the one in the |
expm1(x::Double64)
loses all digits of accuracy belowabs(x) < 1e-163
and returns0.0
instead ofx
.This regresses even beyond Float64, which correctly returns
x
.The text was updated successfully, but these errors were encountered: