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

[seblabbe]

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.

CC: @paulmasson

Component: numerical

Issue created by migration from https://trac.sagemath.org/ticket/21788

[videlec]

comment:1

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 2⁵³ which is roughly 16 digits in base 10.

[videlec]

comment:2

You can also see that

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

[videlec]

comment:3

But

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

[rwst]

comment:5

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.

[rwst]

comment:6

Replying to @videlec:

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.

[videlec]

comment:7

Replying to @rwst:

Replying to @videlec:

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.

[videlec]

comment:8

Replying to @rwst:

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

[seblabbe]

comment:9

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.

[rwst]

comment:10

Replying to @seblabbe:

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.

[seblabbe]

comment:11

Replying to @rwst:

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.