Ask Your Question

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

edit retag close merge delete

Comments

Trying "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?

( 2016-01-29 01:10:23 -0600 )edit

This 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 http://www.wolframalpha.com/input/?i=log%281024%29%2Flog%2810%29+*+82553493450+-+248510777753+%3E+0 (on Wolfram Alpha) 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?

( 2016-01-29 01:24:29 -0600 )edit

1 answer

Sort by » oldest newest most voted

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

more

Comments

Very 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/refer...) but it seems we need to write our own wrapper for this.

( 2016-01-29 10:20:49 -0600 )edit

No, 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.

( 2016-01-29 10:30:22 -0600 )edit

The title Exact Real Arithmetic is indeed very misleading !

( 2016-01-29 10:31:35 -0600 )edit

For RLF, the issue seems to come from tests within the symbolic ring.

( 2016-01-29 10:34:24 -0600 )edit

I 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 RealIntervalFields); 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!

( 2016-01-29 10:46:17 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Stats

Asked: 2016-01-28 21:58:13 -0600

Seen: 1,045 times

Last updated: Jan 29 '16