ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Fri, 29 Jan 2016 10:46:17 -0600Incorrect result for comparison (precision issues?)http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/ Consider this session:
sage: if log(10) * 248510777753 < log(1024) * 82553493450:
....: print 'Less'
....: else:
....: print 'Not less'
....:
Not less
or more simply:
sage: bool(log(10)*248510777753 < log(1024)*82553493450)
False
But this is wrong, as we can see with higher-precision arithmetic:
sage: import mpmath
sage: mpmath.mp.dps = 22 # or anything greater
sage: bool(mpmath.log(10)*248510777753 < mpmath.log(1024)*82553493450)
True
I guess this is happening because Sage is computing to some finite precision. But when writing some bigger program, it's scary that a condition involving variables, like say,
if m * q < n * p:
can without warning give the wrong result and take the wrong path. Is there a way to prevent this from happening, i.e. to make sure that in the program, comparisons are done using as many bits of precision as are necessary to evaluate them correctly, without us having to pre-specify a precision (which may be both too large and wasteful, or too small and give incorrect results)?Thu, 28 Jan 2016 21:58:13 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/Comment by ShreevatsaR for <p>Consider this session:</p>
<pre><code>sage: if log(10) * 248510777753 < log(1024) * 82553493450:
....: print 'Less'
....: else:
....: print 'Not less'
....:
Not less
</code></pre>
<p>or more simply:</p>
<pre><code>sage: bool(log(10)*248510777753 < log(1024)*82553493450)
False
</code></pre>
<p>But this is wrong, as we can see with higher-precision arithmetic:</p>
<pre><code>sage: import mpmath
sage: mpmath.mp.dps = 22 # or anything greater
sage: bool(mpmath.log(10)*248510777753 < mpmath.log(1024)*82553493450)
True
</code></pre>
<p>I guess this is happening because Sage is computing to some finite precision. But when writing some bigger program, it's scary that a condition involving variables, like say,</p>
<pre><code>if m * q < n * p:
</code></pre>
<p>can without warning give the wrong result and take the wrong path. Is there a way to prevent this from happening, i.e. to make sure that in the program, comparisons are done using as many bits of precision as are necessary to evaluate them correctly, without us having to pre-specify a precision (which may be both too large and wasteful, or too small and give incorrect results)?</p>
http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32373#post-id-32373This is extremely odd:
sage: (log(1024)/log(10) * 82553493450 - 248510777753).n()
-0.0000610351562500000
but the correct value is (with `mpmath.mp.dps = 25`) positive, not negative:
sage: mpmath.log10(1024) * 82553493450 - 248510777753
mpf('5.652367462971596978604793549e-12')
But the odd part is that if I try log(1024)/log(10) * 82553493450 - 248510777753 > 0 [on Wolfram Alpha](http://www.wolframalpha.com/input/?i=log%281024%29%2Flog%2810%29+*+82553493450+-+248510777753+%3E+0) it correctly says True, but shows the approximate form of the difference as exactly the same wrong value: -0.0000610352 !! How are both Sage and Wolfram Alpha getting to the same wrong answer?Fri, 29 Jan 2016 01:24:29 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32373#post-id-32373Comment by ShreevatsaR for <p>Consider this session:</p>
<pre><code>sage: if log(10) * 248510777753 < log(1024) * 82553493450:
....: print 'Less'
....: else:
....: print 'Not less'
....:
Not less
</code></pre>
<p>or more simply:</p>
<pre><code>sage: bool(log(10)*248510777753 < log(1024)*82553493450)
False
</code></pre>
<p>But this is wrong, as we can see with higher-precision arithmetic:</p>
<pre><code>sage: import mpmath
sage: mpmath.mp.dps = 22 # or anything greater
sage: bool(mpmath.log(10)*248510777753 < mpmath.log(1024)*82553493450)
True
</code></pre>
<p>I guess this is happening because Sage is computing to some finite precision. But when writing some bigger program, it's scary that a condition involving variables, like say,</p>
<pre><code>if m * q < n * p:
</code></pre>
<p>can without warning give the wrong result and take the wrong path. Is there a way to prevent this from happening, i.e. to make sure that in the program, comparisons are done using as many bits of precision as are necessary to evaluate them correctly, without us having to pre-specify a precision (which may be both too large and wasteful, or too small and give incorrect results)?</p>
http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32372#post-id-32372Trying "arbitrary precision" doesn't help either:
RIF(10).log() * 248510777753 < RIF(1024).log() * 82553493450
or:
log(RLF(1024), base=10) * 82553493450 > 248510777753
both return False, versus
mpmath.log(1024, b=10) * 82553493450 > 248510777753
which correctly returns True. It's beginning to look like a bug in Sage's `log` function?Fri, 29 Jan 2016 01:10:23 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32372#post-id-32372Answer by tmonteil for <p>Consider this session:</p>
<pre><code>sage: if log(10) * 248510777753 < log(1024) * 82553493450:
....: print 'Less'
....: else:
....: print 'Not less'
....:
Not less
</code></pre>
<p>or more simply:</p>
<pre><code>sage: bool(log(10)*248510777753 < log(1024)*82553493450)
False
</code></pre>
<p>But this is wrong, as we can see with higher-precision arithmetic:</p>
<pre><code>sage: import mpmath
sage: mpmath.mp.dps = 22 # or anything greater
sage: bool(mpmath.log(10)*248510777753 < mpmath.log(1024)*82553493450)
True
</code></pre>
<p>I guess this is happening because Sage is computing to some finite precision. But when writing some bigger program, it's scary that a condition involving variables, like say,</p>
<pre><code>if m * q < n * p:
</code></pre>
<p>can without warning give the wrong result and take the wrong path. Is there a way to prevent this from happening, i.e. to make sure that in the program, comparisons are done using as many bits of precision as are necessary to evaluate them correctly, without us having to pre-specify a precision (which may be both too large and wasteful, or too small and give incorrect results)?</p>
http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?answer=32376#post-id-32376Your attempt to use interval arithmetics was on the right way. However, you should know that `RIF(a) < RIF(b)` returns `False` when both intervals intersect. So what you have to do is to test first whether the two intervals intersect. If yes, then you have to improve the precision, otherwise you can compare them. Here is a possible implementation:
def smaller(a,b):
prec = 53
while True:
R = RealIntervalField(prec)
try:
R(a).intersection(R(b))
except ValueError:
return R(a) < R(b)
else:
print '{} bits of precision is not enough'.format(prec)
prec += 1
Note that i had to use the `try...except` statement since `A.intersection(B)` raises a `ValueError` when both interval are disjoint.
Now, you can do:
sage: a = log(10) * 248510777753
sage: b = log(1024) * 82553493450
sage: smaller(a,b)
53 bits of precision is not enough
54 bits of precision is not enough
55 bits of precision is not enough
56 bits of precision is not enough
57 bits of precision is not enough
58 bits of precision is not enough
59 bits of precision is not enough
60 bits of precision is not enough
61 bits of precision is not enough
62 bits of precision is not enough
63 bits of precision is not enough
64 bits of precision is not enough
65 bits of precision is not enough
66 bits of precision is not enough
67 bits of precision is not enough
68 bits of precision is not enough
69 bits of precision is not enough
70 bits of precision is not enough
71 bits of precision is not enough
72 bits of precision is not enough
73 bits of precision is not enough
74 bits of precision is not enough
75 bits of precision is not enough
76 bits of precision is not enough
77 bits of precision is not enough
True
Of course, you can multiply the precision by two instead of incrementing to be faster. Also, you can add a symbolic equality test, and define an upper bound on the precision in case both numbers are equal but cannot be proven equal, and so on.
Fri, 29 Jan 2016 07:56:09 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?answer=32376#post-id-32376Comment by ShreevatsaR for <p>Your attempt to use interval arithmetics was on the right way. However, you should know that <code>RIF(a) < RIF(b)</code> returns <code>False</code> when both intervals intersect. So what you have to do is to test first whether the two intervals intersect. If yes, then you have to improve the precision, otherwise you can compare them. Here is a possible implementation:</p>
<pre><code>def smaller(a,b):
prec = 53
while True:
R = RealIntervalField(prec)
try:
R(a).intersection(R(b))
except ValueError:
return R(a) < R(b)
else:
print '{} bits of precision is not enough'.format(prec)
prec += 1
</code></pre>
<p>Note that i had to use the <code>try...except</code> statement since <code>A.intersection(B)</code> raises a <code>ValueError</code> when both interval are disjoint.</p>
<p>Now, you can do:</p>
<pre><code>sage: a = log(10) * 248510777753
sage: b = log(1024) * 82553493450
sage: smaller(a,b)
53 bits of precision is not enough
54 bits of precision is not enough
55 bits of precision is not enough
56 bits of precision is not enough
57 bits of precision is not enough
58 bits of precision is not enough
59 bits of precision is not enough
60 bits of precision is not enough
61 bits of precision is not enough
62 bits of precision is not enough
63 bits of precision is not enough
64 bits of precision is not enough
65 bits of precision is not enough
66 bits of precision is not enough
67 bits of precision is not enough
68 bits of precision is not enough
69 bits of precision is not enough
70 bits of precision is not enough
71 bits of precision is not enough
72 bits of precision is not enough
73 bits of precision is not enough
74 bits of precision is not enough
75 bits of precision is not enough
76 bits of precision is not enough
77 bits of precision is not enough
True
</code></pre>
<p>Of course, you can multiply the precision by two instead of incrementing to be faster. Also, you can add a symbolic equality test, and define an upper bound on the precision in case both numbers are equal but cannot be proven equal, and so on.</p>
http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32382#post-id-32382I also tried this `(82553493450 * log(RIF(1024), base=10))._cmp_(248510777753)` and got `TypeError: Cannot convert sage.rings.integer.Integer to sage.rings.real_mpfi.RealIntervalFieldElement`: now I understand thanks to your comment, that a `RealIntervalField` is one with a certain fixed precision (though you can vary the precision and get different `RealIntervalField`s); it is not a single field with arbitrary intervals of variable precision. If that existed, then any integer could be trivially converted to a RealIntervalFieldElement because it's just an interval of size 0 around that integer, which can be represented exactly. Thanks again!Fri, 29 Jan 2016 10:46:17 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32382#post-id-32382Comment by tmonteil for <p>Your attempt to use interval arithmetics was on the right way. However, you should know that <code>RIF(a) < RIF(b)</code> returns <code>False</code> when both intervals intersect. So what you have to do is to test first whether the two intervals intersect. If yes, then you have to improve the precision, otherwise you can compare them. Here is a possible implementation:</p>
<pre><code>def smaller(a,b):
prec = 53
while True:
R = RealIntervalField(prec)
try:
R(a).intersection(R(b))
except ValueError:
return R(a) < R(b)
else:
print '{} bits of precision is not enough'.format(prec)
prec += 1
</code></pre>
<p>Note that i had to use the <code>try...except</code> statement since <code>A.intersection(B)</code> raises a <code>ValueError</code> when both interval are disjoint.</p>
<p>Now, you can do:</p>
<pre><code>sage: a = log(10) * 248510777753
sage: b = log(1024) * 82553493450
sage: smaller(a,b)
53 bits of precision is not enough
54 bits of precision is not enough
55 bits of precision is not enough
56 bits of precision is not enough
57 bits of precision is not enough
58 bits of precision is not enough
59 bits of precision is not enough
60 bits of precision is not enough
61 bits of precision is not enough
62 bits of precision is not enough
63 bits of precision is not enough
64 bits of precision is not enough
65 bits of precision is not enough
66 bits of precision is not enough
67 bits of precision is not enough
68 bits of precision is not enough
69 bits of precision is not enough
70 bits of precision is not enough
71 bits of precision is not enough
72 bits of precision is not enough
73 bits of precision is not enough
74 bits of precision is not enough
75 bits of precision is not enough
76 bits of precision is not enough
77 bits of precision is not enough
True
</code></pre>
<p>Of course, you can multiply the precision by two instead of incrementing to be faster. Also, you can add a symbolic equality test, and define an upper bound on the precision in case both numbers are equal but cannot be proven equal, and so on.</p>
http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32381#post-id-32381For `RLF`, the issue seems to come from tests within the symbolic ring.Fri, 29 Jan 2016 10:34:24 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32381#post-id-32381Comment by tmonteil for <p>Your attempt to use interval arithmetics was on the right way. However, you should know that <code>RIF(a) < RIF(b)</code> returns <code>False</code> when both intervals intersect. So what you have to do is to test first whether the two intervals intersect. If yes, then you have to improve the precision, otherwise you can compare them. Here is a possible implementation:</p>
<pre><code>def smaller(a,b):
prec = 53
while True:
R = RealIntervalField(prec)
try:
R(a).intersection(R(b))
except ValueError:
return R(a) < R(b)
else:
print '{} bits of precision is not enough'.format(prec)
prec += 1
</code></pre>
<p>Note that i had to use the <code>try...except</code> statement since <code>A.intersection(B)</code> raises a <code>ValueError</code> when both interval are disjoint.</p>
<p>Now, you can do:</p>
<pre><code>sage: a = log(10) * 248510777753
sage: b = log(1024) * 82553493450
sage: smaller(a,b)
53 bits of precision is not enough
54 bits of precision is not enough
55 bits of precision is not enough
56 bits of precision is not enough
57 bits of precision is not enough
58 bits of precision is not enough
59 bits of precision is not enough
60 bits of precision is not enough
61 bits of precision is not enough
62 bits of precision is not enough
63 bits of precision is not enough
64 bits of precision is not enough
65 bits of precision is not enough
66 bits of precision is not enough
67 bits of precision is not enough
68 bits of precision is not enough
69 bits of precision is not enough
70 bits of precision is not enough
71 bits of precision is not enough
72 bits of precision is not enough
73 bits of precision is not enough
74 bits of precision is not enough
75 bits of precision is not enough
76 bits of precision is not enough
77 bits of precision is not enough
True
</code></pre>
<p>Of course, you can multiply the precision by two instead of incrementing to be faster. Also, you can add a symbolic equality test, and define an upper bound on the precision in case both numbers are equal but cannot be proven equal, and so on.</p>
http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32380#post-id-32380The title `Exact Real Arithmetic` is indeed very misleading !Fri, 29 Jan 2016 10:31:35 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32380#post-id-32380Comment by tmonteil for <p>Your attempt to use interval arithmetics was on the right way. However, you should know that <code>RIF(a) < RIF(b)</code> returns <code>False</code> when both intervals intersect. So what you have to do is to test first whether the two intervals intersect. If yes, then you have to improve the precision, otherwise you can compare them. Here is a possible implementation:</p>
<pre><code>def smaller(a,b):
prec = 53
while True:
R = RealIntervalField(prec)
try:
R(a).intersection(R(b))
except ValueError:
return R(a) < R(b)
else:
print '{} bits of precision is not enough'.format(prec)
prec += 1
</code></pre>
<p>Note that i had to use the <code>try...except</code> statement since <code>A.intersection(B)</code> raises a <code>ValueError</code> when both interval are disjoint.</p>
<p>Now, you can do:</p>
<pre><code>sage: a = log(10) * 248510777753
sage: b = log(1024) * 82553493450
sage: smaller(a,b)
53 bits of precision is not enough
54 bits of precision is not enough
55 bits of precision is not enough
56 bits of precision is not enough
57 bits of precision is not enough
58 bits of precision is not enough
59 bits of precision is not enough
60 bits of precision is not enough
61 bits of precision is not enough
62 bits of precision is not enough
63 bits of precision is not enough
64 bits of precision is not enough
65 bits of precision is not enough
66 bits of precision is not enough
67 bits of precision is not enough
68 bits of precision is not enough
69 bits of precision is not enough
70 bits of precision is not enough
71 bits of precision is not enough
72 bits of precision is not enough
73 bits of precision is not enough
74 bits of precision is not enough
75 bits of precision is not enough
76 bits of precision is not enough
77 bits of precision is not enough
True
</code></pre>
<p>Of course, you can multiply the precision by two instead of incrementing to be faster. Also, you can add a symbolic equality test, and define an upper bound on the precision in case both numbers are equal but cannot be proven equal, and so on.</p>
http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32379#post-id-32379No, `RIF` is the `RealIntervalField` with 53 bits of precision, so there is no adaptation. For `RLF`, it is a bit harder to explain, it maintains its origin (in the present case, an element of the symbolic ring) as a cache and asks to this cache first.Fri, 29 Jan 2016 10:30:22 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32379#post-id-32379Comment by ShreevatsaR for <p>Your attempt to use interval arithmetics was on the right way. However, you should know that <code>RIF(a) < RIF(b)</code> returns <code>False</code> when both intervals intersect. So what you have to do is to test first whether the two intervals intersect. If yes, then you have to improve the precision, otherwise you can compare them. Here is a possible implementation:</p>
<pre><code>def smaller(a,b):
prec = 53
while True:
R = RealIntervalField(prec)
try:
R(a).intersection(R(b))
except ValueError:
return R(a) < R(b)
else:
print '{} bits of precision is not enough'.format(prec)
prec += 1
</code></pre>
<p>Note that i had to use the <code>try...except</code> statement since <code>A.intersection(B)</code> raises a <code>ValueError</code> when both interval are disjoint.</p>
<p>Now, you can do:</p>
<pre><code>sage: a = log(10) * 248510777753
sage: b = log(1024) * 82553493450
sage: smaller(a,b)
53 bits of precision is not enough
54 bits of precision is not enough
55 bits of precision is not enough
56 bits of precision is not enough
57 bits of precision is not enough
58 bits of precision is not enough
59 bits of precision is not enough
60 bits of precision is not enough
61 bits of precision is not enough
62 bits of precision is not enough
63 bits of precision is not enough
64 bits of precision is not enough
65 bits of precision is not enough
66 bits of precision is not enough
67 bits of precision is not enough
68 bits of precision is not enough
69 bits of precision is not enough
70 bits of precision is not enough
71 bits of precision is not enough
72 bits of precision is not enough
73 bits of precision is not enough
74 bits of precision is not enough
75 bits of precision is not enough
76 bits of precision is not enough
77 bits of precision is not enough
True
</code></pre>
<p>Of course, you can multiply the precision by two instead of incrementing to be faster. Also, you can add a symbolic equality test, and define an upper bound on the precision in case both numbers are equal but cannot be proven equal, and so on.</p>
http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32378#post-id-32378Very useful; thank you very much! I'll use this (will also mark this answer as accepted). But I'm wondering why `log(RIF(1024), base=10)` (or `log(RLF(1024), base=10)`), which is printed as `3.010299956639812?`, does not automatically expand its precision when asked to be compared with `248510777753 / 82553493450` (which is about `3.0102999566398119521…`). I do notice now (thanks to you) that with `RIF` it returns False for both `<` and `>` comparison (whereas with `RLF` it returns True for `==`)… I would have expected Sage to have a true exact real arithmetic type (going by the title at http://doc.sagemath.org/html/en/reference/rings_numerical/index.html#exact-real-arithmetic) but it seems we need to write our own wrapper for this.Fri, 29 Jan 2016 10:20:49 -0600http://ask.sagemath.org/question/32371/incorrect-result-for-comparison-precision-issues/?comment=32378#post-id-32378