Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C Programming > Checking square root algorithm for portability/correctness

Reply
Thread Tools

Checking square root algorithm for portability/correctness

 
 
Clint Olsen
Guest
Posts: n/a
 
      03-21-2005
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;
}
 
Reply With Quote
 
 
 
 
Peter Nilsson
Guest
Posts: n/a
 
      03-22-2005
Clint Olsen wrote:
> 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.

--
Peter

 
Reply With Quote
 
 
 
 
Clint Olsen
Guest
Posts: n/a
 
      03-22-2005
On 2005-03-22, Peter Nilsson <(E-Mail Removed)> 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
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.

>> 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.


Thanks for your comments.

-Clint
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Fix point square root Christian VHDL 5 03-18-2010 02:14 PM
square root sqrt algorithm and flowchart popeyerayaz Java 3 07-13-2008 10:14 PM
'big square root' for BigDecimal Jeremy Watts Java 4 05-26-2005 09:01 PM
SRT DIvision, Square root and reciprocal square root alghazo@siu.edu VHDL 0 05-27-2004 06:23 AM
Square Root of floating point number Luca VHDL 1 04-29-2004 02:51 PM



Advertisments