Checking square root algorithm for portability/correctness

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

  1. Clint Olsen

    Clint Olsen Guest

    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;
    }
     
    Clint Olsen, Mar 21, 2005
    #1
    1. Advertising

  2. 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
     
    Peter Nilsson, Mar 22, 2005
    #2
    1. Advertising

  3. Clint Olsen

    Clint Olsen Guest

    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
    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
     
    Clint Olsen, Mar 22, 2005
    #3
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Luca
    Replies:
    1
    Views:
    1,042
    salman sheikh
    Apr 29, 2004
  2. Replies:
    0
    Views:
    1,240
  3. Christian

    Fix point square root

    Christian, Apr 25, 2005, in forum: VHDL
    Replies:
    5
    Views:
    6,667
    jeppe
    Mar 18, 2010
  4. Jeremy Watts

    'big square root' for BigDecimal

    Jeremy Watts, May 26, 2005, in forum: Java
    Replies:
    4
    Views:
    3,572
    Boudewijn Dijkstra
    May 26, 2005
  5. popeyerayaz
    Replies:
    3
    Views:
    2,486
    Arne Vajhøj
    Jul 13, 2008
Loading...

Share This Page