Re: Checking square root algorithm for portability/correctness

From: Peter Nilsson (airia_at_acay.com.au)
Date: 03/22/05 

Date: 21 Mar 2005 20:46:47 -0800

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

Relevant Pages

  • Checking square root algorithm for portability/correctness
    ... I posted a thread on comp.programming awhile back asking about an algorithm ... I implemented on square root. ... prime number as a convenient way to get a pseudorandom number generator. ... iterations could be arbitrary. ...
    (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)