|
|
Message-ID: <20240618122357.GL3766212@port70.net>
Date: Tue, 18 Jun 2024 14:23:57 +0200
From: Szabolcs Nagy <nsz@...t70.net>
To: Damian McGuckin <damianm@....com.au>
Cc: MUSL <musl@...ts.openwall.com>
Subject: Re: roundf() (and round(), and ...)
* Damian McGuckin <damianm@....com.au> [2024-06-17 11:48:38 +1000]:
> Before I submit any suggestion for a change to roundf/round/.., could I get
> a critique of the style to make sure I compy with the style guide.
>
> Also, I have not used anything except for the now default FLT_EVAL_METHOD
> so I need guidance as to where to use 'float_t'. However, the only place
> where such arithmetic (here) might be affected is the expression 'a + a'
> which is exact anyway so I think it is irrelevant.
note that it is unspecified if your algorithm raises inexact or not.
(iirc musl asm implementations don't raise inexact, the c code does,
c23 now requires no inexact which i guess is what you tried to follow,
but that is hard to do in c)
otherwise i think your style is fine, but i added some comments.
FLT_EVAL_METHOD is not relevant for this code.
benchmark data may be useful (or code size e.g. on a soft-float target)
because that may be a valid justification to use this implementation.
>
> Thanks - Damian
>
> Code follows:
>
> float roundf(float x)
> {
> static const int b = 0x7f;
> const float a = fabsf(x);
> union { float f; uint32_t _f; } r = { a }, _x = { x };
>
> if (((int) (r._f >> 23)) < b + 23) /* i.e. effectively |x| < 2^(p-1) */
> {
> /*
> * capture the sign of the argument
> */
> const uint32_t s = _x._f - r._f;
i used to use explicit unions then switched to helper function asuint(x)
because i found that a bit clearer and compiler optimizes it just fine.
and some ppl complained that memcpy is more correct than union, so it is
better hidden away behind a helper function if somebody wants to switch.
nowadays i'd probably write
if (asuint(a) >> 23 < asuint(0x1p23f) >> 23)
but e.g. i find it just as clear writing
const uint32_t abits = asuint(a);
if (abits >> 23 < 0x7f + 23)
> /*
> * this should achieve rounding to nearest with any
> * ties (half-way cases) being rounded away-from-zero.
> * (is it wise to use uint32_t instead of int32_t here?)
> */
> const uint32_t rf = ((uint32_t) (a + a)) - (uint32_t) a;
>
> x = (r.f = (float) rf, r._f |= s, r.f);
it is isa dependent if int32_t or uint32_t is better, but i'd
expect signed convert is more widely supported than unsigned.
so i'd write this as
const float r = (float)((int32_t)(a + a) - (int32_t)a);
x = asfloat(asuint(r) | s);
> }
> return x;
> }
>
> Note that as per the latest IEEE-754 standard, the above does not raise an
> exception in the event of the rounding being inexact. This is not backwards
it is unspecified if (uint32_t)a raises inexact exception,
so this will be target and compiler dependent in practice.
ieee-754-2008:
5.8 Details of conversions from floating-point to integer formats
...
A language standard that permits implicit conversions or expressions
involving mixed types should require that these be implemented with
the inexact-signaling conversion operations below.
c23:
F.4 Floating to integer conversion
...
... whether conversion of a non-integral floating value raises the
"inexact" floating-point exception is unspecified.[note]
...
note: In those cases where it matters, library functions can be used
to effect such conversions with or without raising the "inexact"
floating-point exception.
> compatible with the existing MUSL equivalents. While it is another issue
> altogether, this latest standard also drops the nearbyint() in favour of
> routines called (as per the C standard) 'roundevenf()/roundeven()/...'.
Powered by blists - more mailing lists
Confused about mailing lists and their use? Read about mailing lists on Wikipedia and check out these guidelines on proper formatting of your messages.