#12557 closed defect (fixed)
RDF(1e-17).log() gives NaN
Reported by: | mariah | Owned by: | AlexGhitza |
---|---|---|---|
Priority: | minor | Milestone: | sage-5.1 |
Component: | basic arithmetic | Keywords: | sd40.5 |
Cc: | Merged in: | sage-5.1.beta3 | |
Authors: | Volker Braun | Reviewers: | Douglas McNeil |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description
Using sage-4.8 on skynet/eno
sage: RDF(1e-16).log() -36.7368005697 sage: RDF(1e-17).log() NaN sage: log(1e-17) -39.1439465808988
Attachments (1)
Change History (22)
comment:1 Changed 8 years ago by
comment:2 Changed 8 years ago by
The problem is this branch:
if self._value < 2: if self._value == 0: return RDF(-1)/RDF(0) if self._value < 0: return RDF(0)/RDF(0) sig_on() a = self._new_c(gsl_sf_log_1plusx(self._value - 1) / log_of_base) sig_off() return a
When self._value is small enough, self._value-1 becomes -1 and we get the log of 0. Frankly, even before this happens we lose precision unnecessarily:
sage: [(log(x.n(1000))-RDF(x).log().n()).n(53) for x in [1/10**i for i in [0..17]]] [0.000000000000000, 0.000000000000000, -8.88178419700125e-16, -8.88178419700125e-16, 1.10134124042816e-13, 4.55102622254344e-12, -2.87556645162113e-11, 5.26355847796367e-10, -5.02475927532942e-9, 2.82819314634253e-8, -8.27403674463767e-8, -8.27403674463767e-8, 0.0000221219648111060, -0.000310896853829234, 0.000799597430194865, 0.000799597430194865, -0.104560918227634, NaN]
Probably the simplest fix is to use the 1+x trick only when the argument is actually close to 1. Simply changing the condition to
if 0.5 <= self._value < 2:
will change the above results to the less exciting but more useful
[0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, .000000000000000, 0.000000000000000, 0.000000000000000]
Any objections? Performance seems similar for values which previously would have used gsl_sf_log_1plusx but now simply use gsl_sf_log.
comment:3 Changed 7 years ago by
Okay, candidate patch attached.
comment:4 Changed 7 years ago by
- Status changed from new to needs_review
comment:5 Changed 7 years ago by
- Keywords sd40.5 added
- Reviewers set to Volker Braun
- Status changed from needs_review to needs_work
Fixes the bug!
I think a better doctest would be to check that it is below one ulp:
sage: def check_error(x): ... log_RDF = CDF(x).log().real() ... log_exact = log(x.n(1000)) ... return abs(log_RDF-log_exact) <= log_RDF.ulp() sage: all( check_error( RDF(2^x) ) for x in range(-100,100) ) True
Also, instead of dividing 0/0
why not return RDF.NaN()
?
comment:6 Changed 7 years ago by
Sure, that looks good. (I'll keep the near-1 test in: want to make sure nothing's broken near the boundary values.) RDF(0)/RDF(0) was there before, so I won't take responsibility for that one. ;-)
comment:7 Changed 7 years ago by
[Not sure why two copies were uploaded; they should be identical.]
comment:8 Changed 7 years ago by
- Status changed from needs_work to needs_review
comment:9 Changed 7 years ago by
I tried CDF(x).log().real()
for fun but of course want to test RDF(x).log()
, thats my fault.
The whole branch calling gsl_sf_log_1plusx(self._value - 1)
for values close to 1 nonsense. Sure the gsl_sf_log_1plusx
function is more accurate, but now all the error comes from the numerically unstable self._value - 1
.
Also the doctests fail for me, and I need to set the bound at 2*ulp for it to pass.
Finally, you must doctest check_error(x)
and not check_error(RDF(x))
. Otherwise, the exact log value (for the comparison) is computed at the wrong value:
sage: x = 0.999 sage: x.parent() Real Field with 53 bits of precision sage: x.n(100) # as expected 0.99900000000000000000000000000 sage: RDF(x).n(100) # not the intended value of x.n(100) 0.99899999999999999911182158030
Updated patch follows...
comment:10 Changed 7 years ago by
Too early in the morning ;-)
We have to use check_error(RDF(x))
and not check_error(x)
so that the arbitrary-precision floating point value is the same as the double precision floating point number, not the other way round. I fixed this in the updated patch.
Now also the doctest works with <= 1*ulp
precision!
comment:11 Changed 7 years ago by
Hey, I'm in no position to object: I was trying to preserve the "1+1e-100" case when we're already in RDF and there *is* no "1+1e-100", only 1, which doesn't make a great deal of sense. :-)
Works for me. Positive review.
comment:12 Changed 7 years ago by
- Reviewers changed from Volker Braun to Douglas McNeil
- Status changed from needs_review to positive_review
comment:13 Changed 7 years ago by
- Status changed from positive_review to needs_work
The documentation needs work:
docstring of sage.rings.real_double.RealDoubleElement.log:22: WARNING: Literal block expected; none found.
comment:14 Changed 7 years ago by
Fixed in the updated patch. Only the documentation of the log()
method changed.
comment:15 Changed 7 years ago by
- Status changed from needs_work to needs_review
comment:16 Changed 7 years ago by
Still not right, Volker - we want a single colon after "TESTS" and a double colon after the Trac mention. Or did you upload the wrong one?
comment:17 Changed 7 years ago by
Sorry, I uploaded the wrong patch. Thats sleep deprivation for you ;) Fixed now.
comment:18 Changed 7 years ago by
- Status changed from needs_review to positive_review
comment:19 Changed 7 years ago by
- Merged in set to sage-5.1.beta3
- Resolution set to fixed
- Status changed from positive_review to closed
comment:20 Changed 7 years ago by
See #13134 for a Solaris SPARC failure related to this ticket.
comment:21 Changed 7 years ago by
The problem at #13134 is actually due to a bug in ulp()
. That ticket is now ready for review.
I can confirm this exact behavior on Mac OS X 10.6 and on sage.math.