Opened 5 years ago

# (1/(1006987929*pi - 3163545880)).n() raises division by zero error

Reported by: Owned by: slabbe major sage-7.5 numerical paulmasson N/A

### Description

As reported on sage-devel, using 7.4.beta6, I get

```sage: (1/(1006987929*pi - 3163545880)).n()
Traceback (most recent call last)
...
ValueError: power::eval(): division by zero
```

But providing digits works:

```sage: (1/(1006987929*pi - 3163545880)).n(digits=50)
3.2389954278022058869923921595406901102762355827180e6
```

Maybe also related:

```sage: a = 1/(1006987929*pi - 3163545880)
sage: b = 1/3/(333362599*pi - 1047289492)
sage: bool(a<b)
Traceback (most recent call last)
...
RuntimeError: ECL says: Error executing code in Maxima: expt: undefined: 0 to a negative exponent.
```

### comment:1 Changed 5 years ago by vdelecroix

It works with digits provided the number is big enough

```sage: (1/(1006987929*pi - 3163545880)).n(digits=14)
Traceback (most recent call last):
...
ValueError: power::eval(): division by zero
```

The default precision is 253 which is roughly 16 digits in base 10.

### comment:2 follow-up: ↓ 6 Changed 5 years ago by vdelecroix

You can also see that

```sage: 1 / (1006987929*RR.pi() - 3163545880)
+infinity
```

### comment:3 Changed 5 years ago by vdelecroix

But

```sage: 1 / (1006987929*RealField(64).pi() - 3163545880)
3.23904019306184012e6
```

### comment:5 follow-up: ↓ 8 Changed 5 years ago by rws

```sage: (1006987929*pi - 3163545880).n(54)
4.76837158203125e-7
sage: (1006987929*pi - 3163545880).n(53)
0.000000000000000
```

This is just one of these border cases. The only thing you can do atm is raising precision, which just sets the bar higher.

The only real solution would be to recognize that pi is irrational and raise the precision automatically. Drawback is that it could become slow.

### comment:6 in reply to: ↑ 2 ; follow-up: ↓ 7 Changed 5 years ago by rws

You can also see that

```sage: 1 / (1006987929*RR.pi() - 3163545880)
+infinity
```

Do you want infinity returned in this case too? I'll open a ticket if so.

### comment:7 in reply to: ↑ 6 Changed 5 years ago by vdelecroix

You can also see that

```sage: 1 / (1006987929*RR.pi() - 3163545880)
+infinity
```

Do you want infinity returned in this case too?

This would be very wrong. In the case of floating point there is a dramatic cancelation.

### comment:8 in reply to: ↑ 5 Changed 5 years ago by vdelecroix

```sage: (1006987929*pi - 3163545880).n(54)
4.76837158203125e-7
sage: (1006987929*pi - 3163545880).n(53)
0.000000000000000
```

This is just one of these border cases. The only thing you can do atm is raising precision, which just sets the bar higher.

The only real solution would be to recognize that pi is irrational and raise the precision automatically. Drawback is that it could become slow.

+1. `pi` is even transcendental. Hence for polynomials in `pi` with rational coefficients we should just increase precision (knowing in advance that it can not be zero).

Last edited 5 years ago by vdelecroix (previous) (diff)

### comment:9 follow-up: ↓ 10 Changed 5 years ago by slabbe

Maybe I created that ticket too fast: I wrongly thought that inputs given to `numerical_approx` was about the number of digits (or bits) of precision of the output, not of the intermediate computations.

Maybe that ticket could be closed as won't fix. Or maybe this ticket goal could be to improve the documentation of `numerical_approx` to explain the meaning of the inputs, because I believe this is misleading:

```   Return a numerical approximation of "self" with "prec" bits (or
decimal "digits") of precision.
```
Last edited 5 years ago by slabbe (previous) (diff)

### comment:10 in reply to: ↑ 9 ; follow-up: ↓ 11 Changed 5 years ago by rws

Maybe I created that ticket too fast: I wrongly thought that inputs given to `numerical_approx` was about the number of digits (or bits) of precision of the output, not of the intermediate computations.

Correct me but there is no way to get more precision out of operations with numbers of smaller precision, so the precision value is needed from the start.

Maybe that ticket could be closed as won't fix. Or maybe this ticket goal could be to improve the documentation of `numerical_approx` to explain the meaning of the inputs, because I believe this is misleading:

```   Return a numerical approximation of "self" with "prec" bits (or
decimal "digits") of precision.
```

Yes, and I think the division by zero `ValueError` should be caught and rethrown with an additional hint to raise the precision.

### comment:11 in reply to: ↑ 10 Changed 5 years ago by slabbe

Correct me but there is no way to get more precision out of operations with numbers of smaller precision, so the precision value is needed from the start.

The only thing that I would correct in the last sentence is : "so some high enough precision is needed from the start".

Maybe this is not what we want for `numerical_approx` function because it involves more computations (involving real interval arithmetics?) to compute what high enough precision is needed for the computations to get the output with the desired precision. But to me this is what I expected `numerical_approx` to do after reading its documentation:

```sage: def high_enough_precision(X, desired_prec):
....:     prec = desired_prec
....:     while -log(RealIntervalField(prec)(X).diameter(),2) < desired_prec:
....:         prec +=1
....:     return prec
....:
sage: high_enough_precision(pi, 53)
54
sage: high_enough_precision(pi-3.14, 53)
65
sage: B = (1/(1006987929*pi - 3163545880))
sage: high_enough_precision(B, 53)
108
```

Then, using 108 bits of precision in the internal computations really gives 53 bits of precision for the output:

```sage: RealIntervalField(108)(B)
3.238995427802206?e6
sage: B.n(digits=50)
3.2389954278022058869923921595406901102762355827180e6
```

Yes, and I think the division by zero `ValueError` should be caught and rethrown with an additional hint to raise the precision.

Yes, it might be an intermediate solution. Because sometimes there is no `ValueError`, but the precision should still be raised like in the case of `pi-3.14`.

Last edited 5 years ago by slabbe (previous) (diff)
Note: See TracTickets for help on using tickets.