Sqrt() innacurate?

Hey guys,
I was trying to help someone with a problem involving finding if a number is a perfect square. I came to the conclusion that

Code:

int root;
int square; 
for(square = 0; square < 100; square++)
{
  root = sqrt(square);
  if(root * root == square)
  {
    cout<<square;
  }
}

Given that square can have any beginning point and ending point within the numeric_limits. His teacher then responded. Due to truncation of sqrt() when dealing with alot of digits you will "miss" some perfect squares. This didn't make sense to me. Any number that requires a truncation of sqrt() wouldn't be a perfect square. To try and demonstrate the error that the teacher said would occur I coded a small program.

Code:

#include <iostream>
#include <stdio.h>
#include <limits>
#include <math.h>
using namespace std;

double y;
int rootfound;
int main(void)
{
  // start at the max value of an integer and find the first square using my method
  for(int x = numeric_limits<int>::max(); x > 0; x--)
  {
    int z = (int)sqrt(x);
    if(z*z == x)
    {
      rootfound = z;
      cout<<"found a perfect square: "<<x<<endl;
      cout<<"its root is "<<z<<endl;
      x = 0;
    } 
  }
  // now test the next perfect square and see if it is over the numeric limits of
  // an integar
  rootfound++;
  cout<<"testing the next number up from this root: "<<endl;
  cout<<"The new root to work with is "<<rootfound<<endl;
  cout<<"Its perfect square is "<<(rootfound*rootfound)<<endl;
  if((rootfound*rootfound) < 0) 
  {
    cout<<"this square was over the numeric limits so we found the first square
                 root lower than limits"<<endl;
  }
  else
    cout<<"this square was under the number limits so the sqrt() function
                caused us to miss a perfect square"<<endl;
 
  system("pause");
}

As I expected the result still finds the first perfect square counting down from the numberic_limits of a int. Can anyone demonstrate how the method used would exclude a perfect square that is within the limits of an int? Further why would sqrt() be in a standard library math.h if it wasn't dependable?

I guess due to floating point inaccuracies, it might not be unheard of for the result being slightly lower, e.g sqrt(9) coming out as 2.999999981.

However, this test program shows that it doesn't happen (at least for me) for the int range.

Code:

#include <iostream>
#include <limits>
#include <math.h>
using namespace std;

typedef int IntType ;
int main(void)
{
    for (IntType i = 0; i <= (IntType)sqrt(std::numeric_limits<IntType>::max()); ++i) {
        IntType result = (IntType)sqrt(i * i);
        if (result != i) std::cout << i << ' ' << result << '\n';
    }
    std::cout << "All\n";
}

Exactly my point. I can't demonstrate it happening within the integar range on a number that is a perfect square. If it can't be demonstrated I dont' understand a teacher imposing the restrictions, but obviously if its his course you have to follow his rules to pass. It just made me interested in trying to find the error to which he is refering. Maybe its not an accurate error rate like only comes back wrong 1/ 1000 times or someting, but in that case doesn't every programming method have a "possible" error rate due to data corruption and the 1/1,000,000 chance that your computer will misfire.

Further any other check for a perfect square not using a sqrt() would seem to become time consuming on the CPU running extra loops ect. What if you were dealing with root numbers 0 - 4000. The only other method for approaching I can think of would be

Code:

for(square = 0; square < 4000; square++)
{
  if(IntToCompare == square * square)
  {
    cout<<IntToCompare<<" is a perfect square.";
  }
}

Just checking this would be cumbersome. Then add the possibility that you are checking if a list of 2000 integars are a perfect square it seems this would take WAY to long in comparison to the method mentioned before.

You could also write your own sqrt (approximation - int precision is enough) which is guaranteed to suit the algorithm.

Or perform a binary search on the range.

Actually, if you use floats then you can run into a problem even without using sqrt.
The following fails when i is 4097 because 16785409 can't be represented as a float.
It instead stores 16785408.f into f:

Code:

	for (unsigned int i=0; i<65536; ++i)
	{
		float f = (float)i*i;
		ASSERT((unsigned int)f == i*i);
	}

That is possibly the kind of thing the teacher was worried about.

Neverthelsss, he/she is perfectly correct that using sqrt is not the ideal way to do this, as that usually entails use of floating point data types, and the issues they bring (performance, accuracy, etc).
You've recognised the need for a cast of the return value of sqrt, but you also need to recognise that the warning it otherwise produces is not merely there to get in your way, and that it really is trying to tell you something that's often worth paying attention to.

What makes you think that a floating point can represent any arbitrary integer precisely?

A double is 64 bits but its range is much, much larger than the range of a 64-bit integer. Therefore, there must be some integers it can't represent.

Right?

Originally Posted by anon

You could also write your own sqrt (approximation - int precision is enough) which is guaranteed to suit the algorithm.

Or perform a binary search on the range.

Yep, I've found that method to be effective as well. Another thing you can do is optimize the search range, such that it is already quite close to the actual square root. This implementation calculates the largest power of two needed to contain the root, eg: if N = 1234 then H (high) = 128... since the actual root is 111, H is obviously a pretty close initial guess. It might be improved upon by working out the minimum lower-bound as well, although I'm not sure that it's really worth the effort! It seems to run faster than the Newton-Raphson, AFAICT.

Code:

size_t isqrt( size_t val )
{
	size_t const
		one = 1;
	size_t
		div,
		res,
		low = one,
		hig = one, 
		tmp = val;	
	if( !val )
		return val;
	while( tmp >>= one )
		++hig;
	hig = one << ( ( hig + ( hig & one ) ) >> one );
	for( ;; )
	{
		res = low + ( ( hig - low ) >> one );
		div = val / res;
		if( div < res )
		{
			if( hig == res )
				break;
			hig = res;
		}
		else if( div > res )
		{
			if( low == res )
				break;
			low = res;
		}
		else
			break;
	}
	return res;
}