Checking square root algorithm for portability/correctness

Discussion in 'C Programming' started by Clint Olsen, Mar 21, 2005.

1. Clint OlsenGuest

Hello:

I implemented on square root. The idea was to use the square root of a
prime number as a convenient way to get a pseudorandom number generator.
So, rather than calculate the square root to try to get a particular
precision answer, you would calculate it to x number of bits. This
particular function takes a radicand and the number of iterations as
arguments. So, for 32-bit unsigned integers, you'd want to iterate 16
times. But for the purposes of pseudorandom generation, the number of
iterations could be arbitrary.

I think I've identified one issue with this technique. Some
experimentatoin seems to indicate it is possible to overflow the difference
operand, which probably makes this algorithm non-sensical after a certain
number of iterations.

The notes for compiling are below. If you have time to take a look, that
would be great.

Thanks,

-Clint

Notes:

Compiling:

% cc -o sqrt -Wall -pedantic -g sqrt.c

If you don't want to see the debug messages, use -DNDEBUG as well.

The code is written so it doesn't make any assumptions about the width of
an unsigned long. It estimates this by taking the width of character and
then multiplying by the sizeof the parameterized type: sqrt_t. I made no
attempt to calculate the value bits of an unsigned long. In real life we'd
use the value bits to denote the proper number of iterations to make. The
odd initialization of shift was to account for some oddball machine that
had an odd number of bits in an unsigned long.

I took the algorithm from "Algorithms for Extracting Square Roots and Cube
Roots", Peng, 1981.

Running:

% sqrt 9 16

Usually you want to perform 16 iterations to produce a correct result for a
machine with 32-bit unsigned longs. However, for the purposes of
generating psuedo-random sequences on primes, I was interested in the
behavior if you specified something more arbitrary, like:

% sqrt 2 40

any.

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>

typedef unsigned long sqrt_t;

int main(int argc, char *argv[])
{
sqrt_t a; /* target of sqrt */
sqrt_t root = 0; /* root */
sqrt_t diff; /* running difference */
sqrt_t a_i, b_i;
unsigned width = CHAR_BIT * sizeof(sqrt_t);
int shift = width & 1 ? width - 1 : width - 2;
int i;

if (argc == 3) {
a = strtoul(argv[1], NULL, 10);
i = strtoul(argv[2], NULL, 10);

/* pass 1 */
a_i = a >> shift & 3;

diff = a_i - 1;

if (diff < a_i) /* overflow check */
root = 1;

for (shift -= 2, --i; i > 0; i--, shift -= 2) {

if (shift >= 0) {
a_i = a >> shift & 3;
} else {
a_i = 0;
}

b_i = diff << 2 | a_i;

if (root & 1)
diff = b_i - (root << 2 | 1);
else
diff = b_i + (root << 2 | 3);

root <<= 1;

if (diff < b_i)
root |= 1;

#ifndef NDEBUG
printf("root: %lu, d: %lu, shift: %d\n", root, diff, shift);
#endif
}
printf("%lu\n", root);
} else {
fprintf(stderr, "usage: %s <target> <iterations>\n", argv[0]);
}

return EXIT_SUCCESS;
}

Clint Olsen, Mar 21, 2005

2. Peter NilssonGuest

Clint Olsen wrote:
> an algorithm I implemented on square root. The idea was to use
> the square root of a prime number as a convenient way to get a
> pseudorandom number generator.
> ...
> The code is written so it doesn't make any assumptions about
> the width of an unsigned long.

But it does assume that there are no padding bits.

> It estimates this by taking the width of character and then
> multiplying by the sizeof the parameterized type: sqrt_t.
> I made no attempt to calculate the value bits of an unsigned
> long.

Why not? It's trivial...

int value_bits = 0;
untype_t x;
for (x = -1, x = x/2+1; x; x >>= 1) value_bits++;

> In real life we'd use the value bits to denote the proper
> number of iterations to make.

Or, as in the sample I gave in comp.programming you can simply
mask the highest eventh bit and shift the mask as you go.

> ...
> Usually you want to perform 16 iterations to produce a
> correct result for a machine with 32-bit unsigned longs.
> However, for the purposes of generating psuedo-random
> sequences on primes, I was interested in the behavior
> if you specified something more arbitrary,

Although, that's not something clc delves into.

> efficiency if you have any.

Assuming you have M value bits, your algorithm performs
M-2 + M-4 + M-6 ... shifts. You could reduce this to
precisely M.

My only other style comment is...

> sqrt_t root = 0; /* root */
> sqrt_t diff; /* running difference */
> sqrt_t a_i, b_i;
> ...
> /* pass 1 */
> a_i = a >> shift & 3;
> diff = a_i - 1;
>
> if (diff < a_i) /* overflow check */
> root = 1;

This last if-statement is an obfuscated way of writing...

root = (a_i != 0);

In other words, your 'overflow' occurs if and only if a_i is 0.

--
Peter

Peter Nilsson, Mar 22, 2005

3. Clint OlsenGuest

On 2005-03-22, Peter Nilsson <> wrote:
> But it does assume that there are no padding bits.

Yes.

> Why not? It's trivial...
>
> int value_bits = 0;
> untype_t x;
> for (x = -1, x = x/2+1; x; x >>= 1) value_bits++;

Yeah, I know it's trivial. I just didn't do it

> Or, as in the sample I gave in comp.programming you can simply
> mask the highest eventh bit and shift the mask as you go.

I'll have to look for your post on this.

> Assuming you have M value bits, your algorithm performs
> M-2 + M-4 + M-6 ... shifts. You could reduce this to
> precisely M.

I'm curious how you'd get around this? You need to extract 2 bits from the
significant bits, you can't just keep a running register and shift off the
bits to the right.

>> if (diff < a_i) /* overflow check */
>> root = 1;

>
> This last if-statement is an obfuscated way of writing...
>
> root = (a_i != 0);
>
> In other words, your 'overflow' occurs if and only if a_i is 0.