🚀 go-pugleaf

In <1993Jan6.093937.7848@canon.co.uk> ads@canon.co.uk (Adam Billyard) writes:

>murrayk@prism.CS.ORST.EDU (the Dodger) writes:


>> Well? Does anyone know of a quick and dirty way to get a reasonably accurate
>> square root of a fixed point integer number? My numbers are stored in a 16.16

>It depends on what "reasonably accurate" means.  I used a Newton-Raphson
>approximation in a 3D pool game I did a while back, and I got enough accuracy 
>and quickly enough for my needs..

>	Adam
>	
>NB If you're looking to do vector normalization then you probably need quite a 
>lot of accuracy.
Here is how I do it. This code is from a flight simulator. I use 16 bit
ints but need to take sqrt of 32bit numbers (the result of squaring the
ints).

The program uses a table for the starting value. In practice I write the
tables to a file which I later include into the source.

#define FASTCALL _fastcall		/* microsoft style */

typedef unsigned char	Uchar;
typedef unsigned int	Uint;
typedef unsigned short	Ushort;
typedef unsigned long	Ulong;

static Uint	far sqrtab0[256] = {0};
static Uint	far sqrtab1[256] = {0};
static Uint	far sqrtab2[256] = {0};
static Uint	far sqrtab3[256] = {0};

static int far
lsqrt (Ulong x)				/* used for initialization only */
{
	long	e;
	Ulong	r, t;

	if (x & 0xffff0000L)
		r = 662 + x / 17916;
	else if (x & 0x0000ff00L)
		r = 3 + x / 70;
	else
		r = 2 + x / 11;

	do {
		t = x / r;
		e = (long)(r - t) / 2;
		r = (r + t) / 2;
	} while (e);
	return ((int)r);
}

static void far
set_lsqrt (void)			/* initialise tables */
{
	int	i;

	for (i = 0; i < 256; ++i) {
		sqrtab3[i] = lsqrt (i*256L*256L*256L);
		sqrtab2[i] = sqrtab3[i] >> 4;
		sqrtab1[i] = sqrtab2[i] >> 4;
		sqrtab0[i] = sqrtab1[i] >> 4;
	}
}

static Uint far FASTCALL
my_sqrt (Ulong x)
{
	register Uint	r, t;
	register int	e, ind;

	ind = 0;
	if (x & 0xff000000L) {
		if (x & 0xc0000000L) {
			ind = 1;
			x /= 4;
		}
		r = sqrtab3[x>>24];
	} else if (x & 0x7fff0000L)
		r = sqrtab2[x>>16];
	else if (x & 0x007fff00L)
		r = sqrtab1[x>>8];
	else
		return (sqrtab0[x]);

	do {
		t = (Uint)(x / r);
		e = (int)(r - t) / 2;
		r -= e;
	} while (e);

	if (ind)
		return (t+r);
	else
		return ((t+r)/2);
}

main ()
{
	set_lsqrt ();

	return (0);
}