Checking square root algorithm for portability/correctness

[Clint Olsen]

Hello:

I posted a thread on comp.programming awhile back asking about 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.
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.

About the code:

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

I'd appreciate any comments about correctness and efficiency if you have
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;
}

 

[Peter Nilsson]

I posted a thread on comp.programming awhile back asking about
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.

I'd appreciate any comments about correctness and
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.

 

[Clint Olsen]

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
radicand (a) for each iteration, and since you have to start with the most
significant bits, you can't just keep a running register and shift off the
bits to the right.

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.

Thanks for your comments.

-Clint