Square Root Function not working

[bludstayne]

Why isn't this working? I always get a value that's too high. 4 would give me the correct value, but 9 up always goes over. It's kind of weird too, every perfect square makes it go over a little less, you can check it out and see what I mean.

float sqrt(unsigned int tar, float accr)
{
    float i;
    for (i=0.0f;i<tar;i+=accr)
    {
        if (i*i >= tar)
            return i;
    }
}

[Paul]

http://www.chiefdelphi.com/forums/sh...618#post221618
I responded to your origonal post. See if it helps you out any.

[seanwitte]

This code is based on an algorithm I found on Microchip's web site. Its for 32 bit integers, but can be adapted for 24 and 16 bit ints by changing the mask and result values.

typedef unsigned long uint32;

uint32 sqrt32(uint32 x)
{	
	uint32 mask = 0x1000;
	uint32 result = 0x800;

	while (mask > 0)
	{
		if ((result * result) > x)
		{
			result &= (~mask);
		}
		mask >>= 1;
		result |= mask;
	}
	return result;
}

[Chris Hibner]

Originally Posted by bludstayne

Why isn't this working? I always get a value that's too high. 4 would give me the correct value, but 9 up always goes over. It's kind of weird too, every perfect square makes it go over a little less, you can check it out and see what I mean.

Do you really need a sqare root? We originally had a square root in our waypoint calculation (to determine the distance from us to the waypoint), but we determined the square root was unessecary. Why? Because we're using feedback control to get to the position AND x^2 + y^2 becomes zero at the same point as sqrt(x^2 + y^2). Anyway, it's just a thought - you may not need the square root afterall.

-Chris

[Chris Hibner]

Originally Posted by bludstayne

Why isn't this working? I always get a value that's too high. 4 would give me the correct value, but 9 up always goes over. It's kind of weird too, every perfect square makes it go over a little less, you can check it out and see what I mean.

float sqrt(unsigned int tar, float accr)
{
    float i;
    for (i=0.0f;i<tar;i+=accr)
    {
        if (i*i >= tar)
            return i;
    }
}

Now, to answer your original question:

It is always high due to the >= that you have in the 'if' statement. It won't return a value until it becomes equal or GREATER THAN your value. The odds of it being equal are quite low. Therefore, you're always going to get something that is higher. In fact, it is IMPOSSIBLE to get a value that is lower than the actual square root.

I see two ways to increase your accuracy (I could probably come up with more, but I just woke up).

The first is to make 'accr' smaller. This will give you more accuracy. However, this will also result in slower code.

Here is another way the should give you decent accuracy without much more computation: linearly interpolate between the two points surrounding the square root. Here is an example:

float sqrt(unsigned int tar, float accr)
{
    float i;
    float i_prev;
    for (i=0.0f;i<tar;i+=accr)
    {
        if (i*i >= tar)
            break;
        i_prev = i;
    }
    // return a linearly interpolated value
    return (tar*tar - i_prev*i_prev)/(i*i - i_prev*i_prev) * 
             (i - i_prev) + i_prev;
}

If you don't want to get that complicated, just average the two values surrounding your target:

float sqrt(unsigned int tar, float accr)
{
    float i;
    float i_prev;
    for (i=0.0f;i<tar;i+=accr)
    {
        if (i*i >= tar)
            break;
        i_prev = i;
    }
    // return the average of the point before and the point after the sqrt
    return (i + i_prev)/2;
}

[kiracofe8]

Originally Posted by bludstayne

Why isn't this working? I always get a value that's too high. 4 would give me the correct value, but 9 up always goes over. It's kind of weird too, every perfect square makes it go over a little less, you can check it out and see what I mean.

float sqrt(unsigned int tar, float accr)
{
    float i;
    for (i=0.0f;i<tar;i+=accr)
    {
        if (i*i >= tar)
            return i;
    }
}

Why this goes over has already been answered, but I should also point out that this algorithm will be incredibly slow, and not particularly accurate. Two things:

First off, it is a cardinal sin to have a floating point value as a loop index. The problem is that floating point numbers always have some round off error. And, each time you add, the round off error increases. So, for example, if tar = 1, and accr = 0.1, after 10 iterations, you may get i = 0.9999999 or you may get i = 1.0000001. If you get the former, you'll end up with the square of 1 = 1.1 b/c you'll go one step too far. If you get the later, you'll get the right answer. But there is no way to know which you'll get ahead of time.

The problem of floating point loop indexes can cause such subtle bugs that some languages do not allow it at all (e.g. recent versions of FORTRAN). The following is the recommended replacement:

float sqrt(unsigned int tar, float accr)
{
    int i;
    float j;

    j = 0;

    for (i=0;j<tar;i++)
    {
        j = i * accr;
        if (j*j >= tar)
            return j;
    }
    return -1; //signal an error.  should never happen, but it is bad form to
                  //simply not return any value.    
}

Secondly: the algorithm itself converges slowly. E.g. if you try to take the square root of, say, 100 with accr = 0.1, it will take a LONG time (1000 iterations) to get a correct answer. Try a different method. The most common method to numerically compute square roots is called Newton's method. Do a search on google for "Newton's method square root" and you should find some good example code. Newton's method will converge to a reasonable value after only a few iterations.

[Paul]

Originally Posted by kiracofe8

Why this goes over has already been answered, but I should also point out that this algorithm will be incredibly slow, and not particularly accurate. Two things:

First off, it is a cardinal sin to have a floating point value as a loop index. The problem is that floating point numbers always have some round off error. And, each time you add, the round off error increases. So, for example, if tar = 1, and accr = 0.1, after 10 iterations, you may get i = 0.9999999 or you may get i = 1.0000001. If you get the former, you'll end up with the square of 1 = 1.1 b/c you'll go one step too far. If you get the later, you'll get the right answer. But there is no way to know which you'll get ahead of time.

The problem of floating point loop indexes can cause such subtle bugs that some languages do not allow it at all (e.g. recent versions of FORTRAN). The following is the recommended replacement:

float sqrt(unsigned int tar, float accr)
{
    int i;
    float j;

    j = 0;

    for (i=0;j<tar;i++)
    {
        j = i * accr;
        if (j*j >= tar)
            return j;
    }
    return -1; //signal an error.  should never happen, but it is bad form to
                  //simply not return any value.    
}

Secondly: the algorithm itself converges slowly. E.g. if you try to take the square root of, say, 100 with accr = 0.1, it will take a LONG time (1000 iterations) to get a correct answer. Try a different method. The most common method to numerically compute square roots is called Newton's method. Do a search on google for "Newton's method square root" and you should find some good example code. Newton's method will converge to a reasonable value after only a few iterations.

I've used Newton's system of finding square roots, and it is pretty fast. Take a look at this code, it only works for ints but it should be good enough for what you need. It takes a a number of interations equal to the square root. For instance, to get the sqr of 16 it would take 4 interations, 100 would be 10 and so on and so forth. See how you like it.

inline int sroot(int n)
{
int c = 0; int a = 1;
do {n -= a; a += 2; ++c;}
while (n > 0);
return c;
}

[mtrawls]

Originally Posted by Chris Hibner

Do you really need a sqare root? We originally had a square root in our waypoint calculation (to determine the distance from us to the waypoint), but we determined the square root was unessecary. Why? Because we're using feedback control to get to the position AND x^2 + y^2 becomes zero at the same point as sqrt(x^2 + y^2). Anyway, it's just a thought - you may not need the square root afterall.

-Chris

If you do find you need to compute sqrt(x^2+y^2), might I suggest the CORDIC algorithm (mentioned here before for sin/cos by Kevin Watson, I believe ... I found it very useful). Here is a very good file explaining it: http://www.andraka.com/files/crdcsrvy.pdf. The nice thing about it, is that it is a "shift and add" algorithm ... which is very fast. Maybe I'm just wierd, but I thought it was cool, and was disappointed when I realized we wouldn't be using it ... like Chris said he didn't need to compute the sqrt, we figured that we might not need to compute trig functions ourselves either.