real literals and IEEE-754

[Zimmermann Paul]

Hi,

[since I am not subscribed to sage-devel, please keep me in cc]

the concept of real literals, which was intended (as far as I understand)
to keep exactly track of inputs like "1e-20", leads to the following:

sage: a=RealField(53)(1e-20)
sage: Reals(200)(a)
1.0000000000000000000000000000000000000000000000000000000000e-20
sage: Reals(200)(a.exact_rational())
9.9999999999999994515327145420957165172950370278739244710772e-21

I consider this to be a bug.

In the C language, when one writes "double a = 1e-20; printf ("%.20e\n", a);",
everybody who is fluent with IEEE-754 2008 knows that 10^-20 is not exactly
representable in the binary64 format, thus nobody is surprised to see the
output 9.99999999999999945153e-21.

Moreover real literals seem to be created only with 53-bit precision:

sage: Reals(200)(RealField(52)(1e-20))
1.0000000000000000956165483594623726717277804672348392078969e-20
sage: Reals(200)(RealField(53)(1e-20))
1.0000000000000000000000000000000000000000000000000000000000e-20
sage: Reals(200)(RealField(54)(1e-20))
1.0000000000000000203849099068359721617286420850111158275023e-20

What about getting rid of real literals?

Paul

[Nils Bruin]

On Monday, January 6, 2014 12:46:19 PM UTC-8, Zimmermann Paul wrote:

What about getting rid of real literals?

I think they exist mainly to let

RealField(200)(1e-20)

work in the first place: The parser needs to generate code to instantiate the constant 1e-20 somewhere and that constant should remember its original string representation to allow convenient creation of an element in RealField(200). Consider for instance

sage: RealField(200)(a+(1e-800))
9.9999999999999994515327145420957165172950370278739244710772e-21

(adding this constant doesn't change the 53-bit representation of a but it does ensure the resulting number is not considered a RealLiteral any more).

For added convenience, real literals pretend to be floating point numbers themselves already (they have just some extra data hanging off them):

sage: parent(a)
Real Field with 53 bits of precision

I think this explains the discrepancy with bit lengths:

sage: type(RealField(53)(a))
<type 'sage.rings.real_mpfr.RealLiteral'>
sage: type(RealField(54)(a))
<type 'sage.rings.real_mpfr.RealNumber'>

Converting into a RealField doesn't particularly try to preserve the "RealLiteral" property. It's just that RealLiterals pretend to already be members of RealField(53) and hence the conversion code decides it can just return the input.

That also explains the answer you get from a.exact_rational(): There is no special exact_rational on RealLiterals, so you're just getting what applies to 53-bit real numbers.

So back to your question: I expect that real literals can't be removed because it would make creating high precision constants too much of a pain.

Here is something I think is really a bug:

sage: parent(RealField(200)(1) + 1e-20)
Real Field with 53 bits of precision

The whole point of a real literal is that it CAN coerce into the highest precision field, contrary to normal float coercion, which coerces into the lowest precision. So it seems to me that we should probably have a "RealLiteral field", which is going to be odd in the sense that normal arithmetic in it results in answers in RealField(53), but which coerces into any RealField(n).

[Zimmermann Paul]

Nils,

> Date: Mon, 6 Jan 2014 13:18:07 -0800 (PST)
> From: Nils Bruin <nbr...@sfu.ca>

>
> On Monday, January 6, 2014 12:46:19 PM UTC-8, Zimmermann Paul wrote:
>
> > What about getting rid of real literals?
> >
>
> I think they exist mainly to let
>
> RealField(200)(1e-20)
>
> work in the first place: The parser needs to generate code to instantiate
> the constant 1e-20 somewhere and that constant should remember its original
> string representation to allow convenient creation of an element in
> RealField(200).

but RealField(200)("1e-20") already does the job:

sage: RealField(200)("1e-20")
1.0000000000000000000000000000000000000000000000000000000000e-20

> Consider for instance
>
> sage: RealField(200)(a+(1e-800))
> 9.9999999999999994515327145420957165172950370278739244710772e-21
>
> (adding this constant doesn't change the 53-bit representation of a but it
> does ensure the resulting number is not considered a RealLiteral any more).

this demonstrates another inconsistency:

sage: a=1e-20
sage: RealField(200)(a)
1.0000000000000000000000000000000000000000000000000000000000e-20
sage: RealField(200)(a+0)
9.9999999999999994515327145420957165172950370278739244710772e-21
sage: RealField(200)(1*a)
9.9999999999999994515327145420957165172950370278739244710772e-21

> So back to your question: I expect that real literals can't be removed
> because it would make creating high precision constants too much of a pain.

I don't think so. We should educate the user:

* if she/he writes 1e-20, she/he should be aware that (like in other languages)
1e-20 is not exactly representable in binary, thus will be approximated to
the current precision

* if she/he wants to define a high precision constant, use "1e-20"

Paul

[Thierry]

Hi,

On Mon, Jan 06, 2014 at 01:18:07PM -0800, Nils Bruin wrote:
> Here is something I think is really a bug:
>
> sage: parent(RealField(200)(1) + 1e-20)
> Real Field with 53 bits of precision

On Mon, Jan 06, 2014 at 10:46:44PM +0100, Zimmermann Paul wrote:
> this demonstrates another inconsistency:
>
> sage: a=1e-20
> sage: RealField(200)(a)
> 1.0000000000000000000000000000000000000000000000000000000000e-20
> sage: RealField(200)(a+0)
> 9.9999999999999994515327145420957165172950370278739244710772e-21
> sage: RealField(200)(1*a)
> 9.9999999999999994515327145420957165172950370278739244710772e-21
>

What if, instead of inheriting from RealField(53), real literals inherit
from RealLazyField() ?

sage: a = RealLazyField()('1e-20')

sage: parent(RealField(200)(1) + a)
Real Field with 200 bits of precision

sage: RealField(200)(a)
1.0000000000000000000000000000000000000000000000000000000000e-20
sage: RealField(200)(a+0)

1.0000000000000000000000000000000000000000000000000000000000e-20
sage: RealField(200)(a*1)
1.0000000000000000000000000000000000000000000000000000000000e-20

Ciao,
Thierry

[Robert Bradshaw]

By that logic, rather than preparsing 1/2 we should force the user to
write ZZ(1) / ZZ(2). I find

sage: RealField(200)(1e-20)
9.9999999999999994515327145420957165172950370278739244710772e-21

surprising, as I explicitly asked for 200 digits of precision, this is
a natural way to write a "literal" of high precision. I'll admit

sage: Reals(200)(RealField(53)(1e-20))
1.0000000000000000000000000000000000000000000000000000000000e-20

is probably wrong. BTW, it was intended not just for real fields but
exact ones as well (thought that seems to have changed, as now

sage: QQ(1e-20)
1/99999999999999990439
sage: ZZ(1e30)
1000000000000000019884624838656

though this should be an easy fix. (Yes, I see why that's true, but
for someone who's not a numerical analyst it's quite unfortunate.)

(FWIW, the problem with putting literals in another field, like
RealLazyField, is what to do when doing arithmetic like 1.2 + 1.7 one
wants arithmetic with reals of unspecified parent to be "fast." The
literalness of floating point literals is supposed to just affect
constructors. It almost makes one want to encode the desired precision
in the literal itself....)

- Robert

[Marc Mezzarobba]

Nils Bruin wrote:
> Here is something I think is really a bug:
>
> sage: parent(RealField(200)(1) + 1e-20)
> Real Field with 53 bits of precision

Yes. And I'd say the behavior of RealField(200)(RR(1e-20))) is a bug as
well (the same bug in fact).

So I think there are two different issues here:

1. Real literals are elements of RR. IMO this is a bug that needs to be
fixed. However, Jeroen disagreed in the ticket where this discussion
originated[1], so I would like to hear other people's opinion. In my
view, literals such as 1e-20 could be preparsed to Python floating-
point decimals, rationals, elements of CLF... but certainly not
elements of RR.

[1] http://trac.sagemath.org/ticket/15542

2. Should real literals exist in the first place? I personally agree
with Paul that they do more harm than good (they add a gratuitous
difference between Sage and most languages, including Python, and
they do not fit too well in Sage's floating-point arithmetic model).
But having them in Sage is not an unreasonable design choice either.
And of course getting rid of them entirely would break compatibility
with earlier versions of Sage.

--
Marc

[Marc Mezzarobba]

Robert Bradshaw wrote:
> FWIW, the problem with putting literals in another field, like
> RealLazyField, is what to do when doing arithmetic like 1.2 + 1.7 one
> wants arithmetic with reals of unspecified parent to be "fast." The
> literalness of floating point literals is supposed to just affect
> constructors.

What about putting them in some ad hoc parent in which all arithmetic
operations (except perhaps unary -) would automatically convert their
operands to 53-bits floating-point numbers--and return elements of RR?

(I'm not saying I'm in favor of this solution: I would prefer real
literals to disappear entirely, or to become decimal floating point
numbers if really needs be.)

--
Marc

[Jeroen Demeyer]

The more I think about this, the more I agree with Paul. Real literals
were invented for one particular special case, namely
RealField(200)(1e-20). But from the moment you try something else with
them, stuff breaks. I don't mind RealField(200)(1e-20) to return
9.9999999999999994515327145420957165172950370278739244710772e-21
and hopefully it educates the user about floating-point precision.

I think the inconsistency noted above is quite bad because it is a
counterexample to a == b => f(a) == f(b).

[Volker Braun]

On Monday, January 6, 2014 6:50:05 PM UTC-10, Robert Bradshaw wrote:

By that logic, rather than preparsing 1/2 we should force the user to
write ZZ(1) / ZZ(2).

+1  It is important that the Sage REPL be easy to use with "mathematical" notation as far as possible. You should be able to input high precision floats in a natural way.

But there is also a fundamental difference to the real literals here, "1/2" becomes a honest element of the rational ring whereas real literals just masquerade as elements of RR (using it as their parent) while not actually being 53-bit floats. In particular, this makes RR(1e-20) again a real literal (since the parent is the same), which is quite unexpected.

Would it be possible to move real literals to their own parent? Arithmetic operations with other parents should always cast the literal to the other parent first, and arithmetic operations within the real literals would end up in RR. 

[Nils Bruin]

On Saturday, January 11, 2014 6:26:59 AM UTC-8, Volker Braun wrote:

Would it be possible to move real literals to their own parent? Arithmetic operations with other parents should always cast the literal to the other parent first, and arithmetic operations within the real literals would end up in RR. 

That should be pretty straightforward, but it has an ugly wart:

(hypothetical):

hypothetical_sage: a=RR(71)(1e20)
hypothetical_sage: b=1e0
hypothetical_sage: a+b+b
1.00000000000000000002e20
hypothetical_sage: b+b+a
 1.00000000000000e20

The fact that floating point operations aren't associative/commutative is one thing, but making it so that the resulting parent even changes is quite possibly a worse wart than having to write quotes around literals.

From an advertising point of view I think some form of RealLiteral is a great idea: It makes it really easy to demonstrate: look sage has multiprecision floats and they are "natural" objects for it (no quotes needed). It's just fundamentally underspecified if there is not a "preferred" precision out there. Most other computer algebra system solve that problem by having a preferred precision: a global setting. But it's nasty to have global settings, of course (especially if you want to turn code developed in an interactive session into library code)

In effect we already have a global setting for this, i.e., this already works:

sage: RealNumber=lambda s: RealField(71)(s)
sage: 1e0+1e0+1e20
1.00000000000000000002e20
sage: RealNumber=lambda s: RealField(53)(s)
sage: 1e20+1e0+1e0
1.00000000000000e20

and this would continue to work if we remove RealLiterals. So perhaps we just need better tools to (temporarily) set RealNumber and/or change the name that the preparser produces? We could then deprecate/remove the RealLiterals that we have now.

[Robert Bradshaw]

On Sat, Jan 11, 2014 at 9:11 AM, Nils Bruin <nbr...@sfu.ca> wrote:
> On Saturday, January 11, 2014 6:26:59 AM UTC-8, Volker Braun wrote:
>>
>> Would it be possible to move real literals to their own parent? Arithmetic
>> operations with other parents should always cast the literal to the other
>> parent first, and arithmetic operations within the real literals would end
>> up in RR.
> That should be pretty straightforward, but it has an ugly wart:
>
> (hypothetical):
>
> hypothetical_sage: a=RR(71)(1e20)
> hypothetical_sage: b=1e0
> hypothetical_sage: a+b+b
> 1.00000000000000000002e20
> hypothetical_sage: b+b+a
> 1.00000000000000e20

I think the primary wart here is that RR(1e20) is still a RealLiteral,
leading to inconsistent behavior. I would propose we could special
case this such that this cast actually changes it to a
RealNumber(prec=53) which wouldn't have the wart above or the
strangeness first mentioned in this ticket. (This could perhaps be
done by making its own Parent, but I don't think we should coerce to
the other on arithmetic--if the "true Parent" isn't specified quickly
and explicitly we should treat is as a RealNumber(prec=53).

> The fact that floating point operations aren't associative/commutative is
> one thing, but making it so that the resulting parent even changes is quite
> possibly a worse wart than having to write quotes around literals.
>
> From an advertising point of view I think some form of RealLiteral is a
> great idea: It makes it really easy to demonstrate: look sage has
> multiprecision floats and they are "natural" objects for it (no quotes
> needed). It's just fundamentally underspecified if there is not a
> "preferred" precision out there. Most other computer algebra system solve
> that problem by having a preferred precision: a global setting. But it's
> nasty to have global settings, of course (especially if you want to turn
> code developed in an interactive session into library code)
>
> In effect we already have a global setting for this, i.e., this already
> works:
>
> sage: RealNumber=lambda s: RealField(71)(s)
> sage: 1e0+1e0+1e20
> 1.00000000000000000002e20
> sage: RealNumber=lambda s: RealField(53)(s)
> sage: 1e20+1e0+1e0
> 1.00000000000000e20
>
> and this would continue to work if we remove RealLiterals. So perhaps we
> just need better tools to (temporarily) set RealNumber and/or change the
> name that the preparser produces? We could then deprecate/remove the
> RealLiterals that we have now.

I generally like to avoid modifying the global environment to produce
a local result. A with statement context could be a handy way to get
this though.

- Robert