Checking square root algorithm for portability/correctness

From: Clint Olsen (clint_at_0lsen.net)
Date: 03/21/05 

Date: Mon, 21 Mar 2005 12:32:59 -0600

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;
}

Relevant Pages

  • Re: Checking square root algorithm for portability/correctness
    ... > an algorithm I implemented on square root. ... > number of iterations to make. ... mask the highest eventh bit and shift the mask as you go. ... Assuming you have M value bits, your algorithm performs ...
    (comp.lang.c)
  • Re: Estimating errors
    ... N is the number of iterations and epsilon the precesion I think... ... djamel wrote: ... I have symmetric matrix call V and I compute its square root, ...
    (comp.soft-sys.matlab)
  • Inaccuracy of ANSI clock() as ?-RNG
    ... I tried to build a portable ANSI C stopwatch with sub-/millisecond ... accuracy, using the clock-function. ... Inaccuracy is reflected as a varying number of iterations which fit ... "Juuso's inaccurate clock" - algorithm. ...
    (sci.crypt)
  • Re: Finding the Formula...
    ... algorithm for generating the coefficients. ... how you said it took many many iterations with huge matrices to ... If my "hunch" is right and the polynomials are of degree n, ...
    (sci.math)
  • Re: Square root algorithms and complexity
    ... > I developed an algorithm which includes multiplications,addistions, ... > and the computation of a square root. ... > number of multiplications that the algorithm performs, ... >> computing the square root using a specified algorithm? ...
    (sci.math)