Hereâs a function that computes unsigned long square roots using Newtonâs
method.
I believe that I have tested it quite carefully and am quite confident that it
works correctly.
It returns the unsigned long exact square root value if there is one.
Otherwise, it returns the closest unsigned long int which is less than the
exact real square root value.
For example, it would return
0 as the square root of 0
1 as the square root of 1, 2 and 3
2 as the square root of 4, 5, 6, 7 and 8
3 as the square root of 9, 10 ⦠15
4 as the square root of 16 ⦠24
etc.
I believe that it properly handles all unsigned long values of n from 0L to
((unsigned long)~0L) == 2^64 - 1.
There are probably ways of making it faster.
I am not including my test software because I would be much more confident
that the code works if it was independently tested.
-Danny
#define DELTA(a,b) ( (a) > (b) ? (a) - (b) : (b) - (a) )
unsigned long
ulong_sqrt( unsigned long n )
{
if ( n == 0 || n == 1 ) {
return n;
}
unsigned long guess = n >> 1;
unsigned long quotient = n / guess;
while ( 1 ) {
unsigned long gap = DELTA( guess, quotient );
if ( gap < 2 ) {
break;
}
guess = ( quotient + guess ) >> 1;
quotient = n / guess;
}
guess = guess > quotient ? quotient : guess;
return guess;
}
> On Sep 8, 2015, at 00:03 , Otto Moerbeek <[email protected]> wrote:
>
> On Mon, Sep 07, 2015 at 06:20:02PM -0400, Raul Miller wrote:
>
>> On Mon, Sep 7, 2015 at 6:03 PM, Ted Unangst <[email protected]> wrote:
>>> Michael Warmuth-Uhl wrote:
>>>> On 2015-09-07, Otto Moerbeek wrote:
>>>> It seems to fix the issues on amd64, but I'm not sure if the accuracy of
>>>> long double and sqrtl is guaranteed to be enough in general.
>>>
>>> using floating point seems suspect. i think what's needed is an integer
sqrt
>>> function.
>>
>> Floating point works fine for integers as long as your integers fit
>> within the mantissa component of the floating point representation.
>>
>> Assuming the usual IEEE-754 implementation:
>>
>> For double, this means integers cannot exceed 2^53 (9007199254740992)
>> without losing least significant bits (2^53-1 can be represented, 2^53
>> can be represented but has lost the least significant bit so typically
>> that means it's assumed to be zero - which is fine, but that means
>> that 1+2^53 cannot be represented).
>>
>> For long double, this means integers cannot exceed 2^113
>> (10384593717069655257060992658440192) without losing least significant
>> bits.
>>
>> That said, you can implement integer square root using binary search.
>
> My bet would be Newton iteration, which should converge faster than
> binary search. ping(8) has an implementation of integer sqrt.
>
> BTW, long double won't fix the problem on al platforms since it is not
> guaranteed to be longer than double.
>
> -Otto
