Opened 4 years ago

Last modified 4 years ago

#21788 new defect

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

Reported by: slabbe Owned by:
Priority: major Milestone: sage-7.5
Component: numerical Keywords:
Cc: paulmasson Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

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.

Change History (11)

comment:1 Changed 4 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: Changed 4 years ago by vdelecroix

You can also see that

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

comment:3 Changed 4 years ago by vdelecroix

But

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

comment:4 Changed 4 years ago by paulmasson

  • Cc paulmasson added

comment:5 follow-up: Changed 4 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: Changed 4 years ago by rws

Replying to vdelecroix:

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 4 years ago by vdelecroix

Replying to rws:

Replying to 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 4 years ago by vdelecroix

Replying to 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.

+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 4 years ago by vdelecroix (previous) (diff)

comment:9 follow-up: Changed 4 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 4 years ago by slabbe (previous) (diff)

comment:10 in reply to: ↑ 9 ; follow-up: Changed 4 years ago by rws

Replying to 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.

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 4 years ago by slabbe

Replying to rws:

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 4 years ago by slabbe (previous) (diff)
Note: See TracTickets for help on using tickets.