Ask Your Question
1

adding real literal and real number of high precision

asked 2012-07-27 13:05:52 +0100

Daniel Friedan gravatar image

updated 2012-07-27 13:14:18 +0100

When Sage is adding a real literal to a real number of high precision, shouldn't it calculate the sum in the high precision ring? Instead, Sage seems to calculate in double precision:

RF=RealField(150); RF

Real Field with 150 bits of precision

RF(0.9 + RF(1e-18))

0.90000000000000002220446049250313080847263336

RF(1.0+ RF(1e-18))

1.0000000000000000000000000000000000000000000

RF(1+ RF(1e-18))

1.0000000000000000010000000000000000000000000

I'm trying to use high precision arithmetic (2658 bits) in Sage to verify some results produced by the high precision semidefinite program solver sdpa_gmp. Sage's treatment of real literals in these calculations has made me anxious about the possibility that I'm overlooking other ways in which the calculations might be unreliable.

Is there anywhere an explanation of Sage's treatment of real literals in high precision arithmetic?

Added: Immediately after posting this question, the list of Related Questions in the sidebar pointed me to question/327/set-global-precision-for-reals where I learned that 'RealNumber = RF' would make all real literals lie in the high precision ring. Still, I wonder why the default behavior is to discard precision that is present in the original real literal.

thanks, Daniel Friedan

thanks for

edit retag flag offensive close merge delete

Comments

1

I'm sure there are technical answers, but I think the basic answer is just that we would go to the *lesser* precision, to be sure we are right, when we combine two numbers of different precision. That makes sense to me. `0.9` and friends are default 53-bit precision, while `1` is an `Integer` and so has arbitrarily high precision and preserves what you want.

kcrisman gravatar imagekcrisman ( 2012-07-27 13:23:49 +0100 )edit

There seem to be inconsistencies: RF(.5 + .4 ) 0.90000000000000002220446049250313080847263336 RF(.3+.7) 1.0000000000000000000000000000000000000000000 If .4 and .5 are added as double precision, why not .3 and .7 ?

Daniel Friedan gravatar imageDaniel Friedan ( 2012-07-27 15:07:08 +0100 )edit

EDIT: my previous comment missed the point of your question. The computation is actually done in double precision both times, but RF(RDF(1.0)) = 1.0000000000000000000000000000000000000000000 while RF(RDF(0.9)) = 0.90000000000000002220446049250313080847263336. Which makes the second addition look like it was done in RF all along.

daniels gravatar imagedaniels ( 2012-07-27 16:07:49 +0100 )edit

But RF(RDF(0.3))= 0.29999999999999998889776975374843459576368332 and RF(RDF(0.7))= 0.69999999999999995559107901499373838305473328 so how can RF(0.3+0.7)= 1.0000000000000000000000000000000000000000000 if the addition is being done in double precision?

Daniel Friedan gravatar imageDaniel Friedan ( 2012-07-27 17:34:13 +0100 )edit

1 Answer

Sort by ยป oldest newest most voted
2

answered 2012-07-27 19:31:31 +0100

daniels gravatar image

[This is not really a full answer to your original question but the comment size limit prevents me from posting it as such.]

I don't quite see where the problem should be with the last example: The addition is performed in RDF, giving RDF(1.0) and only after that the conversion of RDF(1.0) to RF happens.

I'll try to explain what's happening as I understand it. First as a remark: The conversion of an element from RDF (53-bits of precision) to RF (RealField with 150-bits of precision) is of course non-canonical. In Sage terms: There is no coercion from RDF to RF, Sage will only automatically coerce from higher precision to lower.

By writing RF(RDF(0.7)) we explicitly ask Sage to convert a lower precision element of RDF to RF anyway, filling up the remaining digits in whatever way it sees fit. I guess this might be the confusing part, because this means that it is not (necessarily) true that RF(RDF(0.3) + RDF(0.7)) = RF(RDF(0.3)) + RF(RDF(0.7)).

If we do RF(0.3+0.7) this is the same as RF(RDF(0.3)+RDF(0.7)), thus the two numbers are added in RDF, giving RDF(1.0), and then converted to RF. If we do on the other hand RF(0.3) + RF(0.7) then 0.3 and 0.7 are interpreted as 150-bit numbers and added in RF (note that they are in fact parsed with the higher precision, they are not first stored as double elements and then converted). Finally and still different, RF(RDF(0.3))+RF(RDF(0.7)) will create the elements with 53-bit precision, then convert them to RF and add them there.

The documentation explains quite nicely how such coercion is performed in general. In particular the explain feature could be interesting for you:

sage: RF = RealField(150)
sage: cm = sage.structure.element.get_coercion_model()
sage: RF(0.3+0.7)
1.0000000000000000000000000000000000000000000
sage: cm.explain(0.3,0.7,add)
Identical parents, arithmetic performed immediately.
Result lives in Real Field with 53 bits of precision
Real Field with 53 bits of precision
sage: RF(RDF(0.3)+RDF(0.7))
1.0000000000000000000000000000000000000000000
sage: cm.explain(RDF(0.3),RDF(0.7),add)
Identical parents, arithmetic performed immediately.
Result lives in Real Double Field
Real Double Field
sage: RF(RDF(0.3)) + RF(RDF(0.3))
0.59999999999999997779553950749686919152736664
sage: cm.explain(RF(RDF(0.3)), RF(RDF(0.7)), add)
Identical parents, arithmetic performed immediately.
Result lives in Real Field with 150 bits of precision
Real Field with 150 bits of precision
sage: RF(0.3) + RDF(0.7)            
1.0
sage: cm.explain(RF(0.3), RDF(0.7), add)
Coercion on left operand via
    Native morphism:
      From: Real Field with 150 bits of precision
      To:   Real Double Field
Arithmetic performed after coercions.
Result lives in Real Double Field
Real Double Field

It may ... (more)

edit flag offensive delete link more

Comments

Thank you very much for the explanations and for the pointers to the documentation. I gather that Sage's interpretation of a token such as '0.7' is context dependent. In RF(0.3+0.4), '0.3' and '0.4' are elements of RDF, and 0.3+0.4=0.7 in RDF, so RF(0.3+0.4)=RF(RDF(0.7)). But in RF(0.7), '0.7' is not an element of RDF, so RF(0.3+0.4) \ne RF(0.7). In order to avoid worrying about this context dependence, I can use 'RealNumber = RF' whenever I need to do high precision calculations where an accidental invocation of RDF would be disastrous. I'll be careful to do this in the future. I must say I find this behavior of Sage to be somewhat peculiar. I assumed that an expression of the form RF ...(more)

Daniel Friedan gravatar imageDaniel Friedan ( 2012-07-27 22:12:08 +0100 )edit

Sort of! This is handled by the preparser. I'm certainly not the right person to ask about that but digging around a bit in the source gave me the following, which at least explains how it works, if not why: When the preparser sees a floating point literal it calls RealNumber(s) with `s` the string it sees. This will give you a rings.real_mpfr.RealLiteral object with precision good enough to represent the literal. If you pass this to RF it will make the right element object. But if you start adding it, it has to to that in some field, and by default it takes RDF. So inputing floating point literals is ok, but as soon as you do mathematical ops with them you end up with an element in RDF instead of a RealLiteral. OTOH 1+pi is symbolic in Sage and can be coerced to any precision.

daniels gravatar imagedaniels ( 2012-07-28 04:47:12 +0100 )edit

You can see how it drops from a RealLiteral to an element of RDF once you do operations (the parent() decides where the operation takes place - the RealLiteral is actually stored with all the precision that your input string has): sage: x = 0.3 sage: type(x), x.parent() (<type 'sage.rings.real_mpfr.realliteral'="">, Real Field with 53 bits of precision) sage: y = x + 0.7 sage: type(y), y.parent() (<type 'sage.rings.real_mpfr.realnumber'="">, Real Field with 53 bits of precision)

daniels gravatar imagedaniels ( 2012-07-28 04:50:26 +0100 )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

Stats

Asked: 2012-07-27 13:05:52 +0100

Seen: 43,364 times

Last updated: Jul 27 '12