Using floating-point math to do integral sqrt

[wyn...@gmail.com]

I saw some code using floating-point math to perform integral square root.
Is it good? It seems if only the mantissa can store the integer.

#include <math.h>

uint16_t sqrt(uint16_t n) {
return ::sqrt(static_cast<double>(n));
};
uint32_t sqrt(uint32_t n) {
return ::sqrt(static_cast<double>(n));
};
uint64_t sqrt(uint64_t n) {
return ::sqrt(static_cast<long double>(n));
};

[Bonita Montero]

This is a down-rounded integral sqrt-function:

#include <cstdint>

uint32_t sqrt( std::uint32_t x )
{
using namespace std;
int n;
uint32_t squareRoot,
squareRootNew,
square,
squareNew;
for( int n = 15, squareRoot = 0, square = 0; n >= 0; n-- )
{
squareRootNew = squareRoot + ((uint32_t)1 << n);
squareNew = square + ((squareRoot + squareRootNew) << n);
if( squareNew <= x )
squareRoot = squareRootNew,
square = squareNew;
}
return squareRoot;
}

[Bonita Montero]

Am 02.08.2020 um 17:41 schrieb Bonita Montero:
> This is a down-rounded integral sqrt-function:
>
> #include <cstdint>
>
> uint32_t sqrt( std::uint32_t x )
> {
>     using namespace std;
>     int      n;
>     uint32_t squareRoot,
>              squareRootNew,
>              square,
>              squareNew;
>     for( int n = 15, squareRoot = 0, square = 0; n >= 0; n-- )

for( n = 15, squareRoot = 0, square = 0; n >= 0; n-- )

[Bo Persson]

Define "good". :-)

If the mantissa uses 52 bits, the first two would work ok.

The last one depends on what a long double really is, and/or if all 64
bits of n are actually used.


Bo Persson

[Mike Terry]

I have done this (or equivalent) in my own code, because it was the most
efficient algorithm I could find to achieve calculation of integral
square roots.

The reason is that floating point square roots have hardware support and
execute (at least in my environment) considerably quicker than anything
else I found.

The downside is that due to resolution issues, the floating point square
root is NOT always equal to the required integral square root, so
FURTHER PROCESSING IS REQUIRED TO TEST FOR THIS AND CORRECT IF REQUIRED.
[The floating result was sometimes out by one from the correct
result.] Even so, in timed tests this (with the follow-up corrections)
was clearly quicker than any other approach I looked at.

(So, it /might/ be good, but probably only as part of some larger
routine that handles any required corrections. Maybe on other
architectures or with other argument sizes the correction would not be
required. My environment was Windows 10 on a DELL laptop with Intel i7,
calculating square roots for 64-bit arguments. No, I don't think I have
the timing results any more...)

Mike.

[Jorgen Grahn]

Since the subject mentions sorting, you might as well sort on
sqrt(n)*sqrt(n), i.e. n.

Unless you need the loss of precision that makes e.g. integer sqrt(4)
and sqrt(5) compare equal.

/Jorgen

--
// Jorgen Grahn <grahn@ Oo o. . .
\X/ snipabacken.se> O o .

[Christian Gollwitzer]

Am 03.08.20 um 09:29 schrieb Jorgen Grahn:

> Since the subject mentions sorting, you might as well sort on
> sqrt(n)*sqrt(n), i.e. n.
>

Where do you see "sorting"? I see "store" and "sqrt" only.

Christian

[Öö Tiib]

Result can be one less. Those would be more round:

uint16_t sqrt(uint16_t n) {
return ::sqrt(static_cast<double>(n)) + 0.5;
};
uint32_t sqrt(uint32_t n) {
return ::sqrt(static_cast<double>(n)) + 0.5;
};
uint64_t sqrt(uint64_t n) {
return ::sqrt(static_cast<long double>(n)) + 0.5;
};

[Jorgen Grahn]

Oh, I misread. Thanks. I read "integral sort" when I browed past it
quickly yesterday, and for some reason that "sort" stuck in my head.

Maybe I need new glasses after all.

[wyn...@gmail.com]

Mike Terry於 2020年8月3日星期一 UTC+8上午8時58分09秒寫道：

Uh, it took me a while to understand what you were talking about was
representation issues!

[Bonita Montero]

> uint64_t sqrt(uint64_t n) {
> return ::sqrt(static_cast<long double>(n));
> };

Which compiler today has long double support > 64 bit ?
I think they're rather rare.

[jacobnavia]

Under the ARM64 bits, gcc has long double of 128 bits implemented in
software. I have tried to implement that also in my compiler but it
isn't ready yet.

My compiler lcc-win has 80 bit long double support under x86 (32 and 64
bits)

[Juha Nieminen]

Bonita Montero <Bonita....@gmail.com> wrote:
> Which compiler today has long double support > 64 bit ?
> I think they're rather rare.

For Intel targets gcc and clang (and probably icc) have generated 80-bit
floating point for long doubles since the dawn of time. (This is
probably configurable with a compiler option, though.)

I think ARM64 specifies support for 128-bit floating point, but
whether this is supported by the actual CPU hardware is up to the
CPU manufacturer (and, AFAIK, the vast majority don't). In that
case it's usually simulated in software.

[Manfred]

The documentation for MSVC still lists FP10.obj [*] as a link option to
change "the default precision control to 64 bits" which means 80-bit
intermediate values, however the linked to page "Floating-Point Support"
[**] has removed the FP10.obj option from recent versions of VS.

[*] https://docs.microsoft.com/en-us/cpp/c-runtime-library/link-options
[**]
https://docs.microsoft.com/en-us/cpp/c-runtime-library/floating-point-support

[Bonita Montero]

> I think ARM64 specifies support for 128-bit floating point, but
> whether this is supported by the actual CPU hardware is up to the

> CPU manufacturer (and, AFAIK, the vast majority don't). ...
According to the Wikipedia the POWER9-CPU is the only one with 128 bit
fp hardware-support. I think there's rare need for 128 bit fp, but no
need for 128 bit fp with hw-speed.

[Chris M. Thomasson]

Deep fractal zooms would benefit from fast 128-bit floats in the GPU.

[Bonita Montero]

>> According to the Wikipedia the POWER9-CPU is the only one with 128 bit
>> fp hardware-support. I think there's rare need for 128 bit fp, but no
>> need for 128 bit fp with hw-speed.

> Deep fractal zooms would benefit from fast 128-bit floats in the GPU.

There are even climate-predictions made with endless chain-calculations.
And no GPU would ever have 128 bit fp.

[Chris M. Thomasson]

Why? They have 64 bit...

https://docs.nvidia.com/cuda/floating-point/index.html

To get really deep fractal zooms, people create arbitrary precision in
the shader.

[Bonita Montero]

>> There are even climate-predictions made with endless chain-calculations.
>> And no GPU would ever have 128 bit fp.

> Why? They have 64 bit...

Because there's no need for it.

[Öö Tiib]

No one cares what you think they need. They just need.
They need 128 and even 256 bit floating points:
<https://link.springer.com/article/10.1007/s10915-012-9679-3>

[Bonita Montero]

>> Because there's no need for it.

> No one cares what you think they need. They just need.
> They need 128 and even 256 bit floating points:
> <https://link.springer.com/article/10.1007/s10915-012-9679-3>

If there would be need for >= 128 bit FP in GPUs there would
be >= 128 bit FP in GPUs.

[Öö Tiib]

Sure. And so your sentence that "no GPU would ever have 128 bit
fp." is false. As pointed article told, people did 128 and 256
bit fp on whole battery of GPUs 8 years ago.

[Bonita Montero]

> Sure. And so your sentence that "no GPU would ever have 128 bit
> fp." is false.

No, it's not false. 128-bit fp-supports means native support,
not emulated support.

[Chris M. Thomasson]

128 native should be a heck of a lot faster than emulating it using,
say, two 64 bit floats.