[musl] Double Rounding in ....

[musl] Double Rounding in ....

[Damian McGuckin]

 	hypot, hypotf, hypotl
and
 	scalbn, scalbnf, scalbnl

Does anybody have a good/clear explanation for how double rounding occurs 
in the case of squaring-of/multiplication-by small numbers and how the 
action taken in the code addresses that? Not sure if it needs to mention

 	FLT_EVAL_METHOD

because certainly the float and double versions of those routines are
written in terms of float_t and double_t respectively.

Thanks - Damian

Pacific Engineering Systems International ..... 20D Grose St, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer

Re: [musl] Double Rounding in ....

[Szabolcs Nagy]

* Damian McGuckin <damianm@esi.com.au> [2023-11-11 19:11:41 +1100]:
> 
> 	hypot, hypotf, hypotl
> and
> 	scalbn, scalbnf, scalbnl
> 
> Does anybody have a good/clear explanation for how double rounding occurs in
> the case of squaring-of/multiplication-by small numbers and how the action
> taken in the code addresses that? Not sure if it needs to mention

n * 2^n is normally exact, but in the subnormal range there
are less precision bits so rounding can occur.

if you do the scaling in two steps (e.g. because the scale
factor has such a big exponent that it is not representible
in a single float) then you can get double rounding.

e.g.

scalbn(0x1.000000000000bp-1, -1024) == 0x1.0000000000008p-1025

normally the scaling is done in two steps: x*0x1p-1022*0x1p-2
because up to 0x1p-1022 the scaling is easy to construct from an int,
while 0x1p-1024 is subnormal.

with naive implementation the bit pattern at the end is
..001011 x
..00110  x*0x1p-1022 (rounded)
..010    x*0x1p-2 (rounded)

with the fix
..001011 x*0x1p-969 (969 = 1022-53)
..001    x*0x1p-55 (rounded)

so simply split the scale factor such that the first scaling can
only do rounding if the second scaling (<0x1p-53) completely rounds
the entire result away.

> 
> 	FLT_EVAL_METHOD
> 
> because certainly the float and double versions of those routines are
> written in terms of float_t and double_t respectively.
> 
> Thanks - Damian
> 
> Pacific Engineering Systems International ..... 20D Grose St, Glebe NSW 2037
> Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
> Views & opinions here are mine and not those of any past or present employer