mpn_sqrtrem{1,2}

[Marco Bodrato]

Ciao,

Il Ven, 24 Marzo 2017 11:54 pm, Adrien Prost-Boucle ha scritto:
> This page, at the end, tends to discourage playing with rounding mode too
> frequently:
> https://www.gnu.org/software/libc/manual/html_node/Rounding.html

This page is about writing portable code using libc, here we have a quite
specific hardware: CPUs supporting SSE/SSE2.

Il Ven, 24 Marzo 2017 10:07 pm, Adrien Prost-Boucle ha scritto:
> But the tradeoff may not be in our favor. I didn't test.

We should eventually, but at first we have to ensure correctness and to
correctly select when this code can be compiled and used.

> Now, why float with 23-bit mantissa is enough for 32-bit inputs:

> The worst corner case (hypothetically) would be that the root is first
> too low by 1 because of the mantissa too short, and then too low by an
> additional -1 because of rounding mode.

We round towards 0 (or -inf) too. Maybe the worst risk might arise for the
other direction. But I do not know exactly how rounding on floats works.

Let N = 2147580932 = 46342^2 - 1, and rounding be toward +infinity.
Will N be rounded up to 8388989*256 when converted to float?
Then its root 11863552.60.../256 rounded up to 11863553/256?
and this number converted back to integer rounding it up to 46343?
If so, we have a possible +2.

Best regards,
m

-- 
http://bodrato.it/papers/