RDF(1e-17).log() gives NaN

[mariah]

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

[kcrisman]

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

[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.

[dsm]

Okay, candidate patch attached.

[vbraun]

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()?

[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. ;-)

[dsm]

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

[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...

[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!

[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.

[jdemeyer]

The documentation needs work:

docstring of sage.rings.real_double.RealDoubleElement.log:22: WARNING: Literal block expected; none found.

[vbraun]

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

[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?

[vbraun]

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

[vbraun]

Updated patch

[jdemeyer]

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

[jdemeyer]

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