Unexpected behavior of roundTo:

CONTENTS DELETED

CONTENTS DELETED

CONTENTS DELETED

Once more, Floats are inexact.
    0.01 ~= 0.01s2
The name #absPrintExactlyOn:base: is lying, it does not print exactly,
but it prints the shortest decimal representation than will be rounded
to the same Float when read back.

To print it exactly, you need to use #printShowingDecimalPlaces: indeed.
As every finite Float is a represented internally as a Fraction with a
denominator being a power of 2, every finite Float has a decimal
representation with a finite number of decimals digits (just multiply
numerator and denominator with adequate power of 5, and you'll get the
digits).

So try:
0.01 printShowingDecimalPlaces: 59
-> 0.01000000000000000020816681711721685132943093776702880859375

You see that even if you try to execute the operation without rounding
error, then convert it back to Float, you get the error:
 (2.8011416510246336 asTrueFraction roundTo: 0.01 asTrueFraction) asFloat
->  2.8000000000000003

When you perform the #roundTo: operations in Float inexact arithmetic,
you may accumulate more rounding errors, so the result may vary.

If you want to round to an exact hundredth, then use exact arithmetic and try:
   2.8011416510246336 roundTo: 0.01s2

Hope it helps to consider Float for what they are, inexact but fast.
Don't put too much expectations on them.

Nicolas

2011/4/8 Hilaire Fernandes <[hidden email]>:

+1.  This seems like an extension of "never test equality on floats."



________________________________________
From: [hidden email] [[hidden email]] On Behalf Of Nicolas Cellier [[hidden email]]
Sent: Friday, April 08, 2011 1:14 PM
To: [hidden email]
Subject: Re: [Pharo-project] Unexpected behavior of roundTo:

Once more, Floats are inexact.
    0.01 ~= 0.01s2
The name #absPrintExactlyOn:base: is lying, it does not print exactly,
but it prints the shortest decimal representation than will be rounded
to the same Float when read back.

To print it exactly, you need to use #printShowingDecimalPlaces: indeed.
As every finite Float is a represented internally as a Fraction with a
denominator being a power of 2, every finite Float has a decimal
representation with a finite number of decimals digits (just multiply
numerator and denominator with adequate power of 5, and you'll get the
digits).

So try:
0.01 printShowingDecimalPlaces: 59
-> 0.01000000000000000020816681711721685132943093776702880859375

You see that even if you try to execute the operation without rounding
error, then convert it back to Float, you get the error:
 (2.8011416510246336 asTrueFraction roundTo: 0.01 asTrueFraction) asFloat
->  2.8000000000000003

When you perform the #roundTo: operations in Float inexact arithmetic,
you may accumulate more rounding errors, so the result may vary.

If you want to round to an exact hundredth, then use exact arithmetic and try:
   2.8011416510246336 roundTo: 0.01s2

Hope it helps to consider Float for what they are, inexact but fast.
Don't put too much expectations on them.

Nicolas

2011/4/8 Hilaire Fernandes <[hidden email]>:

Yes nice example to exhibit in school

This again happens even if performed exactly (then rounded to nearest Float)
 (2.8 asTrueFraction roundTo: 0.01 asTrueFraction) asFloat

As soon as you write 0.01 instead of (1/100) or (0.01s2), the worm is
in the fruit.

Nicolas

2011/4/8 Hilaire Fernandes <[hidden email]>:

Pour enfoncer le clou, let's play a bit more with inexact representations:

{
((2.8 asTrueFraction roundTo: 0.01 asTrueFraction) - (2.8
predecessor)) abs -> -1.
((2.8 asTrueFraction roundTo: 0.01 asTrueFraction) - (2.8)) abs -> 0.
((2.8 asTrueFraction roundTo: 0.01 asTrueFraction) - (2.8 successor)) abs -> 1.
} detectMin: [:e | e key ]

you get 0.0->1, which mean that:
  (2.8 asTrueFraction roundTo: 0.01 asTrueFraction) asFloat = (2.8 successor)

But remember:
  (2.8 asTrueFraction roundTo: 0.01 asTrueFraction) ~= (2.8 successor)

It must be interpreted as: the nearest Float to (2.8 asTrueFraction
roundTo: 0.01 asTrueFraction) is
(2.8 successor).

If you want to know how far it is, then get an idea with:

((2.8 asTrueFraction roundTo: 0.01 asTrueFraction) - (2.8 successor)
asTrueFraction) asFloat
 -2.0816681711721685e-16

Nicolas

2011/4/8 Schwab,Wilhelm K <[hidden email]>:

2011/4/8 Nicolas Cellier <[hidden email]>:
Oh one thing I forgot, you wanna know how many digits you need after
decimal point to print exactly ?

Easy, (0.01 asTrueFraction denominator highBit - 1).

Nicolas

CONTENTS DELETED

2011/4/8 Hilaire Fernandes <[hidden email]>:
Previous versions were lying too by not printing a number that would
be reinterpreted into itself.
But if it is what you want, then you should use the old print message,
#absPrintOn:base:
Maybe it could be wrapped into a printApproximatelyOn: to be
polymorphic with other numbers ?

Nicolas

> --
> Education 0.2 -- http://blog.ofset.org/hilaire
>
>
>

On Apr 8, 2011, at 7:48 PM, Hilaire Fernandes wrote:

> Well, whatever the underneath representation, one can expect roundedTo:
> 2 to return a float with two decimals.
>
Sometimes I think we should use the resources that these amazing machines
give us these days to move programming languages closer to humans...

There are better Float models than the ones that are implemented in hardware.
E.g. another INRIA Project is this GNU MPFR:

"The MPFR library is a C library for multiple-precision floating-point computations with correct rounding."

        http://www.mpfr.org/

Why don't we make our language better at "real" math? The power in the machine is definitly there...
I am quite convincd that if there are normal programming languages in 50 years, the math part of them
will be closer to Mathematica then to C...

        Marcus





--
Marcus Denker  -- http://www.marcusdenker.de
INRIA Lille -- Nord Europe. Team RMoD.

2011/4/8 Marcus Denker <[hidden email]>:
mpfr are still floating point with inexact operations, you just change
the precision to whatever number of bits you want.

I played with my own modest Smalltalk implementation relying on
LargeInteger arithmetic
http://www.squeaksource.com/ArbitraryPrecisionFl.html
It implements correct rounding to (if correct is round to nearest even).

Nicolas

>
>
>
> --
> Marcus Denker  -- http://www.marcusdenker.de
> INRIA Lille -- Nord Europe. Team RMoD.
>
>
>

On Apr 8, 2011, at 7:14 PM, Nicolas Cellier wrote:

> Once more, Floats are inexact.
>    0.01 ~= 0.01s2
> The name #absPrintExactlyOn:base: is lying, it does not print exactly,
> but it prints the shortest decimal representation than will be rounded
> to the same Float when read back.

we should update the comment if it does not reflect that.

Yes!

Nicolas

would not you like the write (I can help) a two three 5 pages document
for the next book on this topics? :)

I would be gorgeous and everybody would love that

Stef

On Apr 8, 2011, at 7:18 PM, Nicolas Cellier wrote:

>>
> I played with my own modest Smalltalk implementation relying on
> LargeInteger arithmetic
> http://www.squeaksource.com/ArbitraryPrecisionFl.html
> It implements correct rounding to (if correct is round to nearest even).

Nicolas

would it be a path for you and us so that we use your library?
Ok it will put pressure on you but like that we move one the next age.

Stef

At the end of the day, we still insist on expecting there won't be loss
of information just because we wrote a decimal print string for a binary
floating point number.  One can get offended or irritated all one wants,
but the reality of the situation won't change.

I think the real solution to these "problems" is to implement decimal
floating point as per IEEE-754/2008.  Then you can write something like
1.2345 and know for a fact the number represented is *exactly* 1.2345.
At least one VisualWorks platform already supports decimal floating
point in hardware (IBM's POWER line).  C99 extensions exist at least in
draft form to extend C so that it supports decimal floating point.  IBM
has released a GPL library in ANSI C that implements the feature.  I am
not sure what is the deal with the license and if you can e.g.:
interface to it without making your whole Smalltalk GPL.

Food for thought...

On 4/8/11 12:53 , Marcus Denker wrote:

On Apr 9, 2011, at 9:09 AM, Andres Valloud wrote:

> At the end of the day, we still insist on expecting there won't be loss of information just because we wrote a decimal print string for a binary floating point number.  One can get offended or irritated all one wants, but the reality of the situation won't change.
>
> I think the real solution to these "problems" is to implement decimal floating point as per IEEE-754/2008.  \\

The MPFR lib seems to me (as far as I understand) is inspired by IEEE 754

<SNIP>
MPFR is a portable library written in C for arbitrary precision arithmetic on floating-point numbers. It is based on the GNU MP library. It aims to provide a class of floating-point numbers with precise semantics. The main characteristics of MPFR, which make it differ from most arbitrary precision floating-point software tools, are:

        • the MPFR code is portable, i.e., the result of any operation does not depend on the machine word size mp_bits_per_limb (64 on most current processors);
        • the precision in bits can be set exactly to any valid value for each variable (including very small precision);
        • MPFR provides the four rounding modes from the IEEE 754-1985 standard, plus away-from-zero, as well as for basic operations as for other mathematical functions.

....


The following five rounding modes are supported:

        • MPFR_RNDN: round to nearest (roundTiesToEven in IEEE 754-2008),
        • MPFR_RNDZ: round toward zero (roundTowardZero in IEEE 754-2008),
        • MPFR_RNDU: round toward plus infinity (roundTowardPositive in IEEE 754-2008),
        • MPFR_RNDD: round toward minus infinity (roundTowardNegative in IEEE 754-2008),
        • MPFR_RNDA: round away from zero (experimental).

The ‘round to nearest’ mode works as in the IEEE 754 standard: in case the number to be rounded lies exactly in the middle of two representable numbers, it is rounded to the one with the least significant bit set to zero. For example, the number 2.5, which is represented by (10.1) in binary, is rounded to (10.0)=2 with a precision of two bits, and not to (11.0)=3. This rule avoids the drift phenomenon mentioned by Knuth in volume 2 of The Art of Computer Programming (Section 4.2.2).
</SNIP>


> Then you can write something like 1.2345 and know for a fact the number represented is *exactly* 1.2345. At least one VisualWorks platform already supports decimal floating point in hardware (IBM's POWER line).  C99 extensions exist at least in draft form to extend C so that it supports decimal floating point.  IBM has released a GPL library in ANSI C that implements the feature.  I am not sure what is the deal with the license and if you can e.g.: interface to it without making your whole Smalltalk GPL.
>

MPFR is LGPL, so linking that in the VM or using it via FFI is no problem license wise.

        Marcus


--
Marcus Denker  -- http://www.marcusdenker.de
INRIA Lille -- Nord Europe. Team RMoD.

On 4/9/11 0:28 , Marcus Denker wrote:
>
> On Apr 9, 2011, at 9:09 AM, Andres Valloud wrote:
>
>> I think the real solution to these "problems" is to implement
>> decimal floating point as per IEEE-754/2008.  \\
>
> The MPFR lib seems to me (as far as I understand) is inspired by IEEE
> 754
 >
 > <quote from the docs>

The docs seem to suggest MPRF implements variable precision binary
floating point numbers.  IEEE-754/2008, as opposed to IEEE-754/1985,
also specifies types for base 10 (i.e.: decimal) floating point numbers.
  In other words, 2.5 wouldn't be represented as (101)_2, but (25)_10.
Thus, going back to the original problem that caused this discussion,

> 2.8011416510246336 roundTo: 0.01
>
> returns
>
>  2.8000000000000003

the problem here is that there is no binary floating point number that
can represent the number 0.01.  The arguments why this is so are the
usual (and I just put them here for the sake of completeness): binary
floating point numbers represent fractions of the form

+/- m * 2^k

where m, k are integers such that 0 <= m < M, and e <= k <= K.  So let's
try to find which values of m and k match for 0.01:

1/100 = m * 2^k

Well, m is not zero, so it must be at least 1, and so k must be
negative.  Set k' = -k and rearrange:

2^k' = 100m

Now we invoke the help of the fundamental theorem of arithmetic that
says integers can be factored into primes in only one way.  On the left
we only have 2 as a prime factor.  On the right, there is a factor of 5.
  Therefore, we conclude the above equality is impossible.  So, m and k
do not exist, and so no binary floating point number is equal to 0.01.

To make the story short, you cannot assume lossless operation just
because you can see a base 10 print string.  If loss of precision is the
issue because we only have 2^k' on the left, let's try base 10 instead.
  In other words, floating point numbers in which the base is 10
represent values of the form

+/- m * 10^k

where m, k are integers such that 0 <= m < M, and e <= k <= K.  So let's
try to find which values of m and k match for 0.01:

1/100 = m * 10^k

So now we just match the fraction on the left with m=1 and k=-2.  Done.
  Of course, there are values that cannot be represented without loss of
precision, but the loss of precision involves only m and k.  For example,

67423523487542376543765376234576235726345623472 = m * 10^k

when 0 <= m < M and M = 9.  There are not enough digits in the mantissa,
so the result is going to be 7*10^k for some k.  But this issue happens
as well with binary floating point numbers when k is positive.  Try

67423523487542376543765376234576235726345623472 asFloat truncated

Increasing the precision of binary floating point numbers to arbitrary
precision help you if you need a bigger mantissa, but a bigger mantissa
does not solve the fundamental problem of binary floating point numbers
that causes imprecision whenever you use base 10 to write down values.

> MPFR is LGPL, so linking that in the VM or using it via FFI is no
> problem license wise.

I am not convinced that's what you need:

http://en.wikipedia.org/wiki/Decimal_floating_point

And it looks like Intel wrote a decimal floating point library, too...
now I need to check the license on this one as well:

http://software.intel.com/en-us/blogs/2008/03/06/intel-decimal-floating-point-math-library/

Andres.