Ask Your Question
1

Incorrect result for comparison (precision issues?)

asked 2016-01-29 04:58:13 +0200

ShreevatsaR gravatar image

updated 2017-08-01 12:18:44 +0200

FrédéricC gravatar image

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 flag offensive 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?

ShreevatsaR gravatar imageShreevatsaR ( 2016-01-29 08:10:23 +0200 )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?

ShreevatsaR gravatar imageShreevatsaR ( 2016-01-29 08:24:29 +0200 )edit

1 Answer

Sort by » oldest newest most voted
1

answered 2016-01-29 14:56:09 +0200

tmonteil gravatar image

updated 2016-01-29 15:36:09 +0200

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.

edit flag offensive delete link 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.

ShreevatsaR gravatar imageShreevatsaR ( 2016-01-29 17:20:49 +0200 )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.

tmonteil gravatar imagetmonteil ( 2016-01-29 17:30:22 +0200 )edit

The title Exact Real Arithmetic is indeed very misleading !

tmonteil gravatar imagetmonteil ( 2016-01-29 17:31:35 +0200 )edit

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

tmonteil gravatar imagetmonteil ( 2016-01-29 17:34:24 +0200 )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!

ShreevatsaR gravatar imageShreevatsaR ( 2016-01-29 17:46:17 +0200 )edit

I guess it may be impossible to handle the case when a is actually equal to b, but we can raise an exception: see https://ask.sagemath.org/question/535...

ShreevatsaR gravatar imageShreevatsaR ( 2020-09-24 14:21:23 +0200 )edit

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

Asked: 2016-01-29 04:58:13 +0200

Seen: 1,340 times

Last updated: Jan 29 '16