Opened 10 years ago

Closed 10 years ago

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

Reported by: Owned by: mariah AlexGhitza minor sage-5.1 basic arithmetic sd40.5 sage-5.1.beta3 Volker Braun Douglas McNeil N/A

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

### comment:1 Changed 10 years ago by kcrisman

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

### comment:2 Changed 10 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 10 years ago by dsm

Okay, candidate patch attached.

### comment:4 Changed 10 years ago by dsm

• Status changed from new to needs_review

### comment:5 Changed 10 years ago by vbraun

• 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 10 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 10 years ago by dsm

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

### comment:8 Changed 10 years ago by dsm

• Status changed from needs_work to needs_review

### comment:9 Changed 10 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 10 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 10 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 10 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 10 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 10 years ago by vbraun

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

### comment:15 Changed 10 years ago by vbraun

• Status changed from needs_work to needs_review

### comment:16 Changed 10 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 10 years ago by vbraun

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

Updated patch

### comment:18 Changed 10 years ago by jdemeyer

• Status changed from needs_review to positive_review

### comment:19 Changed 10 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 10 years ago by jdemeyer

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

### comment:21 Changed 10 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.