Opened 7 years ago

Closed 7 years ago

Last modified 7 years ago

#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)

trac_12557_fix_small_RDF_logs_v4.patch (2.6 KB) - added by vbraun 7 years ago.
Updated patch

Download all attachments as: .zip

Change History (22)

comment:1 Changed 7 years ago by kcrisman

I can confirm this exact behavior on Mac OS X 10.6 and on sage.math.

comment:2 Changed 7 years ago by dsm

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 dsm

Okay, candidate patch attached.

comment:4 Changed 7 years ago by dsm

  • Status changed from new to needs_review

comment:5 Changed 7 years ago by vbraun

  • 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 dsm

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 dsm

[Not sure why two copies were uploaded; they should be identical.]

comment:8 Changed 7 years ago by dsm

  • Status changed from needs_work to needs_review

comment:9 Changed 7 years ago by vbraun

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 vbraun

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 dsm

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 dsm

  • Authors set to Volker Braun
  • Reviewers changed from Volker Braun to Douglas McNeil
  • Status changed from needs_review to positive_review

comment:13 Changed 7 years ago by jdemeyer

  • 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 vbraun

Fixed in the updated patch. Only the documentation of the log() method changed.

comment:15 Changed 7 years ago by vbraun

  • Status changed from needs_work to needs_review

comment:16 Changed 7 years ago by kcrisman

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 vbraun

Sorry, I uploaded the wrong patch. Thats sleep deprivation for you ;) Fixed now.

Changed 7 years ago by vbraun

Updated patch

comment:18 Changed 7 years ago by jdemeyer

  • Status changed from needs_review to positive_review

comment:19 Changed 7 years ago by jdemeyer

  • 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 jdemeyer

See #13134 for a Solaris SPARC failure related to this ticket.

comment:21 Changed 7 years ago by jdemeyer

The problem at #13134 is actually due to a bug in ulp(). That ticket is now ready for review.

Note: See TracTickets for help on using tickets.