Ask Your Question
5

What are the differences between RealDoubleField() and RealField(53) ?

asked 2013-04-04 13:34:40 -0500

tmonteil gravatar image

updated 2014-05-22 13:22:44 -0500

Hi,

This question is related to question 2402 (still open!) and tries to collect differences between RDF=RealDoubleField() and RR=RealField(53). These are two floating point real number fields with both 53 bits of precision. The first one comes from the processor floating-point arithmetic, the second one is "emulated" by mpfr. They are assumed to follow the same rounding standards (to the nearest, according to the sagebook, but i may be wrong).

However, we can see some differences between them:

sage: RDF(1/10)*10 == RDF(1)
False
sage: RDF(1/10)*10 - RDF(1)
-1.11022302463e-16

sage: RR(1/10)*10 == RR(1)  
True
sage: sage: RR(1/10)*10 - RR(1) 
0.000000000000000

Could you explain that ? EDIT: this was a bug and it is now fixed, see trac ticket 14416.

There are also some specificities on which field should be used for some methods.

  • For example, it seems that the eignevalues are not well computed on RR, but are correctly computed on RDF (see trac #13660). What is the reason for that ?

  • Also, it seems that when dealing with huge matrices, the fast atlas library in only used when entries are in RDF, not in RR (see trac #10815).

Are there other difference that should be known between both implementations of floating point numbers in Sage ?

Thierry

edit retag flag offensive close merge delete

Comments

It is worth noticing that the exponent is bounded in RDF and not in RR : sage: RDF(1/2^10000) == 0 True sage: RR(1/2^10000) == 0 False

tmonteil gravatar imagetmonteil ( 2013-04-04 13:35:05 -0500 )edit

As for matrices, I think that for a lot of these methods we send things to numpy and friends, which only do things in double precision... something along those lines, I'm not precisely sure. Some methods are even only available for `RDF`.

kcrisman gravatar imagekcrisman ( 2013-04-04 14:48:38 -0500 )edit

4 answers

Sort by ยป oldest newest most voted
1

answered 2013-04-05 00:23:52 -0500

updated 2014-06-16 01:13:43 -0500

Rounding in RR and RDF is not the same, as shown below (or the precision is different, I do not know):

sage: a = RR(1/10)   
sage: a.sign_mantissa_exponent()
(1, 7205759403792794, -56)
sage: (a * RR(10)).sign_mantissa_exponent()
(1, 4503599627370496, -52)

sage: b = RDF(1/10)
sage: b.sign_mantissa_exponent()   # same floating point as above
(1, 7205759403792793, -56)
sage: (b * RDF(10)).sign_mantissa_exponent()  # different from above !!
(1, 9007199254740991, -53)

In particular, RR has the equality 1/10 * 10 = 1 but not RDF because

sage: RR(1).sign_mantissa_exponent()                       
(1, 4503599627370496, -52)
sage: RDF(1).sign_mantissa_exponent()
(1, 4503599627370496, -52)
edit flag offensive delete link more

Comments

1

Indeed, it seems that for RDF, the rounding is not to the nearest as claimed in the [sagebook](http://sagebook.gforge.inria.fr/) (Section 11.3.2). Il looks like more like a rounding toward zero. sage: print "RDF ", RDF(1/10)*10 - RDF(1), RDF(-1/10)*10 + RDF(1) RDF -1.11022302463e-16 1.11022302463e-16 sage: for rounding in ['RNDZ','RNDD','RNDU','RNDN']: ....: R = RealField(prec=53,rnd=rounding) ....: print rounding, R(1/10)*10 - R(1), R(-1/10)*10 + R(1) RNDZ -1.11022302462515e-16 1.11022302462515e-16 RNDD -1.11022302462516e-16 -2.22044604925032e-16 RNDU 2.22044604925032e-16 1.11022302462516e-16 RNDN 0.000000000000000 0.000000000000000

tmonteil gravatar imagetmonteil ( 2013-04-05 00:49:19 -0500 )edit

Great analysis. I wonder what we can do about this.

kcrisman gravatar imagekcrisman ( 2013-04-05 02:48:23 -0500 )edit

In Sage 5.11 I get: sage: print "RDF ", RDF(1/10)*10 - RDF(1), RDF(-1/10)*10 + RDF(1) RDF 0.0 0.0

Eviatar Bach gravatar imageEviatar Bach ( 2013-08-19 12:17:21 -0500 )edit
1

answered 2013-04-05 06:24:28 -0500

tmonteil gravatar image

updated 2013-04-05 07:21:24 -0500

Actually, it seems to be even worse! The sagebook is right in saying that RDF rounds its computations to the nearest:

sage: sage: print "RDF ", RealField(100)(RDF(1)/RDF(10)), RealField(100)(RDF(-1)/RDF(10))
RDF  0.10000000000000000555111512313 -0.10000000000000000555111512313
sage: for rounding in ['RNDZ','RNDD','RNDU','RNDN']:
....:     R = RealField(prec=53,rnd=rounding)
....:     print rounding, RealField(100)(R(1)/R(10)), RealField(100)(R(-1)/R(10))
RNDZ 0.099999999999999991673327315311 -0.099999999999999991673327315311
RNDD 0.099999999999999991673327315311 -0.10000000000000000555111512313
RNDU 0.10000000000000000555111512313 -0.099999999999999991673327315311
RNDN 0.10000000000000000555111512313 -0.10000000000000000555111512313

But the conversion from QQ to RDF rounds towards zero:

sage: print "RDF ", RealField(100)(RDF(1/10)), RealField(100)(RDF(-1/10))
RDF  0.099999999999999991673327315311 -0.099999999999999991673327315311
sage: for rounding in ['RNDZ','RNDD','RNDU','RNDN']:
....:     R = RealField(prec=53,rnd=rounding)
....:     print rounding, RealField(100)(R(1/10)), RealField(100)(R(-1/10))
....:     
RNDZ 0.099999999999999991673327315311 -0.099999999999999991673327315311
RNDD 0.099999999999999991673327315311 -0.10000000000000000555111512313
RNDU 0.10000000000000000555111512313 -0.099999999999999991673327315311
RNDN 0.10000000000000000555111512313 -0.10000000000000000555111512313

See trac ticket 14416.

edit flag offensive delete link more
1

answered 2013-04-06 01:43:32 -0500

Volker Braun gravatar image

Slightly simplified, your question is about the following:

sage: RealField(150)( RDF(1)/RDF(10) )       # case 1
0.10000000000000000555111512312578270211815834
sage: RealField(150)( RR(1)/RR(10) )         # case 2
0.10000000000000000555111512312578270211815834

vs.

sage: R = RealField(prec=53, rnd='RNDZ')
sage: RealField(150)( R(1)/R(10) )           # case 3
0.099999999999999991673327315311325946822762489
sage: RealField(150)( RDF(1/10) )            # case 4
0.099999999999999991673327315311325946822762489

The first and second cases make guarantees about IEEE floating point behavior and therefore must yield the same answer. The third case explicitly states that it is rounding towards zero, and its to be expected that you get a slightly different answer with different rounding. The fourth case punts to the __float__ method of rationals, which is

def __float__(self):
    """
    Return floating point approximation to ``self`` as a Python float.

    OUTPUT: float

    EXAMPLES::

        sage: (-4/17).__float__()
        -0.23529411764705882
        sage: float(-4/17)
        -0.23529411764705882
    """
    return mpq_get_d(self.value)

According to the GMP/MPIR docs, mpq_get_d rounds towards zero so you get the same answer as in the third case.

edit flag offensive delete link more
0

answered 2013-04-04 14:55:45 -0500

kcrisman gravatar image

Here's at least one difference.

sage: RDF(1).ulp()
2.22044604925e-16
sage: RR(1).ulp()
2.22044604925031e-16

sage: a = RDF(1/10)
sage: a.ulp()
1.38777878078e-17
sage: b = RR(1/10)
sage: b.ulp()
1.38777878078145e-17

This seems relevant too; even though it's just printing, I think it might be connected.

sage: b
0.100000000000000
sage: a
0.1
sage: b^10
1.00000000000000e-10
sage: a^10
1e-10
edit flag offensive delete link more

Comments

Indeed, those differences seems only about the string representation of the result, not the numbers themselves that seem equal: sage: RealField(150)(RDF(1).ulp()) 2.2204460492503130808472633361816406250000000e-16 sage: RealField(150)(RR(1).ulp()) 2.2204460492503130808472633361816406250000000e-16 sage: RealField(150)(RDF(1/10).ulp()) 1.3877787807814456755295395851135253906250000e-17 sage: RealField(150)(RR(1/10).ulp()) 1.3877787807814456755295395851135253906250000e-17 Whereas, the difference i mentioned seems to be related in a strange rounding somewhere sage: RealField(150)(RDF(1/10)*10 - RDF(1)) -1.1102230246251565404236316680908203125000000e-16 sage: RealField(150)(RR(1/10)*10 - RR(1)) 0.000000000000000000000000000000000000000

tmonteil gravatar imagetmonteil ( 2013-04-05 00:10:42 -0500 )edit

See this entry in the FAQ on the Sage wiki: [What exactly does Sage do when I type 0.6**2?](http://wiki.sagemath.org/faq#What_exactly_does_Sage_do_when_I_type_0.6.2A.2A2.3F)

slelievre gravatar imageslelievre ( 2013-04-24 00:46:32 -0500 )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: 2013-04-04 13:34:40 -0500

Seen: 1,524 times

Last updated: Jun 16 '14