Mailing List Archive: Float rounding

Float rounding

Aug 17, 2023, 5:44 AM

Post #1 of 4 (89 views)

As mentioned on Amsterdam IRC by Johan (summarized):

sciurius: 'printf "%.2f\n", -1.785;' prints "-1.78" for >= 5.30 and "-1.79" for <= 5.28

We (Johan, Abe and me) have done some test to check if this was a
correct statement and if so, what part of the toolchain were to blame.

I used my perl farm to conclude:

-1.78 voor 5.6.0 … 5.8.0, -1.79 voor 5.8.1 … 5.29.1, -1.78 voor 5.29.2 … 5.38.1

Then I suggested that glibc/libm might be debet, but none of us could
find proof that the version of glibc had influence on this outcome.

I *know* that floats are imprecise, but I promised I'd forward this to
the ML.

Abe then did a bisect to hit the culprit being

ce6f496d720f620645562842: PATCH: [perl #41202] text->float gives wrong answer

Again, *I* have no personal opinion or interest, but Rob and/or people
with knowledge in this area might have.

commit ce6f496d720f6206455628425320badd95b07372
Author: Sisyphus <sisyphus@cpan.org>
Date: Wed Aug 1 22:33:38 2018 +1000

PATCH: [perl #41202] text->float gives wrong answer

This changes to use Perl_strtod() when available, and that turns out to
be the key to fixing this bug.

S_mulexp10() is removed from embed.fnc to avoid repeating the
complicated prerequisites for defining Perl_strtod(). This works
because this static function already was defined before use in
numeric.c, and always called in full form without using a macro.

James Keenan fixed a file permissions problem originally introduced by
this commit, but the fix has been squashed into it.

--
H.Merijn Brand https://tux.nl Perl Monger http://amsterdam.pm.org/
using perl5.00307 .. 5.37 porting perl5 on HP-UX, AIX, and Linux
https://tux.nl/email.html http://qa.perl.org https://www.test-smoke.org

Re: Float rounding [ In reply to ]

sisyphus359 at gmail

Aug 17, 2023, 6:30 AM

Post #2 of 4 (89 views)

Permalink

On Thu, Aug 17, 2023 at 10:45?PM <perl5@tux.freedom.nl> wrote:

> As mentioned on Amsterdam IRC by Johan (summarized):
>
> sciurius: 'printf "%.2f\n", -1.785;' prints "-1.78" for >= 5.30 and
> "-1.79" for <= 5.28
>

I see the same on 64-bit builds of perl on Windows.
According to my checks against the mpfr library, the result, as given by
5.30 is correct, and the result as given by 5.28 is incorrect.

I'm finding that on perl-5.28, the value -1.785 is being assigned
incorrectly:

>perl -wle "printf '%a', -1.785;"
-0x1.c8f5c28f5c29p+0

"%.2f" formatting correctly outputs that incorrect value
(-1.7850000000000001) as -1.79.

whereas perl-5.30 assigns the value correctly:

>perl -wle "printf '%a', -1.785;"
-0x1.c8f5c28f5c28fp+0

"%.2f" formatting correctly outputs that correct value
(-1.7849999999999999) as -1.78.

I expect that "%a" formatting will reveal the same anomaly on the system
Johan was using.

Cheers,
Rob

Re: Float rounding [ In reply to ]

pjacklam at gmail

Aug 19, 2023, 2:05 AM

Post #3 of 4 (87 views)

Permalink

tor. 17. aug. 2023 kl. 15:30 skrev sisyphus <sisyphus359@gmail.com>:
>
> According to my checks against the mpfr library, the result, as
> given by 5.30 is correct, and the result as given by 5.28 is
> incorrect.

What do you get when you use different rounding modes?

When converting the value 1.785 to the "binary64" format, ignoring
the sign, we first strip the leading bit, since this is implicit
for this value in the "binary64" format. We are then left with the
significand 0.785. Since 0.5 <= 0.785 < 1, the exponent will be
zero, so there is no need to adjust for that. The last thing to do
is to represent 0.785 using 52 bits. One way to do this is to
multiply by 2**52 and round the result to an integer.

0.785 * 2**52 = 3535325707485839.36

If we round away from zero, we get the exact value

1 + 3535325707485840 / 2**52
= 1.7850000000000001421085471520200371742249

which can also be written as

0x1.c8f5c28f5c29p+0

If we round towards zero or round to nearest, we get the exact
value

1 + 3535325707485839 / 2**52
= 1.7849999999999999200639422269887290894985

which can also be written as

0x1.c8f5c28f5c28fp+0

So when the internal representation of 1.785 changes from
0x1.c8f5c28f5c29p+0 to 0x1.c8f5c28f5c28fp+0, I suspect this is due
to some change in the rounding done when 1.785 is converted to the
"binary64" format.

Peter

Re: Float rounding [ In reply to ]

sisyphus359 at gmail

Aug 19, 2023, 6:21 AM

Post #4 of 4 (87 views)

Permalink

On Sat, Aug 19, 2023 at 7:06?PM Peter John Acklam <pjacklam@gmail.com>
wrote:

> tor. 17. aug. 2023 kl. 15:30 skrev sisyphus <sisyphus359@gmail.com>:
> >
> > According to my checks against the mpfr library, the result, as
> > given by 5.30 is correct, and the result as given by 5.28 is
> > incorrect.
>
> What do you get when you use different rounding modes?
>

Firstly, I should point out that I was assuming that the perl in question
was one whose nvtype was 'double'. (I think that assumption was correct.)
I also assumed that, being perl, the rounding mode is "to nearest, with
ties to even".

That being so, doing $x=-1.785 should assign the value
-0x1.c8f5c28f5c28fp+0 to $x.
Assigning anything else is simply incorrect.

When converting the value 1.785 to the "binary64" format, ignoring
> the sign, we first strip the leading bit, since this is implicit
> for this value in the "binary64" format. We are then left with the
> significand 0.785. Since 0.5 <= 0.785 < 1, the exponent will be
> zero, so there is no need to adjust for that. The last thing to do
> is to represent 0.785 using 52 bits. One way to do this is to
> multiply by 2**52 and round the result to an integer.
>
> 0.785 * 2**52 = 3535325707485839.36
>
> If we round away from zero, we get the exact value
>
> 1 + 3535325707485840 / 2**52
> = 1.7850000000000001421085471520200371742249
>
> which can also be written as
>
> 0x1.c8f5c28f5c29p+0
>
> If we round towards zero or round to nearest, we get the exact
> value
>
> 1 + 3535325707485839 / 2**52
> = 1.7849999999999999200639422269887290894985
>
> which can also be written as
>
> 0x1.c8f5c28f5c28fp+0
>
> So when the internal representation of 1.785 changes from
> 0x1.c8f5c28f5c29p+0 to 0x1.c8f5c28f5c28fp+0, I suspect this is due
> to some change in the rounding done when 1.785 is converted to the
> "binary64" format.
>

You could well be right.
If that's the sort of procedure that perl was following in perl-5.28 and
earlier, then that may well explain why those perls assigned so many values
incorrectly.

Such random changes of rounding mode could account for discrepancies of
1ULP (as is the case in this particular example), though I've got a feeling
that sometimes the discrepancy would be *more than* 1ULP.
(It's a while ago, and I'm not sure about that.)

Cheers,
Rob