Opened 5 years ago
Last modified 5 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 5 years ago by
comment:2 follow-up: ↓ 6 Changed 5 years ago by
You can also see that
sage: 1 / (1006987929*RR.pi() - 3163545880) +infinity
comment:3 Changed 5 years ago by
But
sage: 1 / (1006987929*RealField(64).pi() - 3163545880) 3.23904019306184012e6
comment:4 Changed 5 years ago by
- Cc paulmasson added
comment:5 follow-up: ↓ 8 Changed 5 years ago by
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
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 5 years ago by
Replying to rws:
Replying to vdelecroix:
You can also see that
sage: 1 / (1006987929*RR.pi() - 3163545880) +infinityDo 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
Replying to rws:
sage: (1006987929*pi - 3163545880).n(54) 4.76837158203125e-7 sage: (1006987929*pi - 3163545880).n(53) 0.000000000000000This 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).
comment:9 follow-up: ↓ 10 Changed 5 years ago by
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.
comment:10 in reply to: ↑ 9 ; follow-up: ↓ 11 Changed 5 years ago by
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 5 years ago by
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
.
It works with digits provided the number is big enough
The default precision is 2^{53} which is roughly 16 digits in base 10.