Author |
Topic: Integer Square Root (Read 1998 times) |
|
John_Gaughan
Uberpuzzler
Behold, the power of cheese!
    
 
Gender: 
Posts: 767
|
 |
Integer Square Root
« on: Apr 27th, 2004, 7:49pm » |
 Quote  Modify
|
This is a problem that has a real application for a program I am writing. I have not been able to come up with a good solution yet.
Given a C++ class:
Code:struct Point
{
int x;
int y;
...
// Calculate the distance between *this point
// and another one, discarding everything after
// the decimal (i.e. floor)
unsigned int distance (const Point &p) const
{
...
}
}; |
|
Space considerations aside (as long as they are not too bad), what is the fastest algorithm to calculate the distance between two Point objects? Keep in mind that coordinates are signed, distance is unsigned (think about it for a minute), and rounded down.
Standard sqrt() functions work on doubles, which is potentially less precise because of significant digits and floating point precision. This function needs to be precisely accurate in all cases.
So far I cannot think of anything that does not use log functions, which are probably not much faster than sqrt().
|
|
 IP Logged |
x = (0x2B | ~0x2B)
x == the_question
|
|
|
Barukh
Uberpuzzler
Gender:
Posts: 2276
|
|
Re: Integer Square Root
« Reply #1 on: Apr 28th, 2004, 12:03am » |
Quote Modify
|
John, are you assuming any constraints on the coordinates of the points?
If not, please note that the distance() ought to return long long, not unsigned int. This may affect the selection of the best algorithm.
|
|
IP Logged |
|
|
|
towr
wu::riddles Moderator
Uberpuzzler
Some people are average, some are just mean.
Gender:
Posts: 13730
|
|
Re: Integer Square Root
« Reply #2 on: Apr 28th, 2004, 1:29am » |
Quote Modify
|
I'd say unsigned long long then..
You could do a binary search for the right number d
( d^2 = (x1-x2)^2+(y1-y2)^2 )
You should get the number in log2(d) steps
There are probably some ways to optimise it.
|
|
IP Logged |
Wikipedia, Google, Mathworld, Integer sequence DB
|
|
|
John_Gaughan
Uberpuzzler
Behold, the power of cheese!
Gender:
Posts: 767
|
|
Re: Integer Square Root
« Reply #3 on: Apr 28th, 2004, 9:56am » |
Quote Modify
|
Yes, you are right, I would need to go up to the next size. The real application uses 16 bit integers so I could go with an unsigned 32 bit integer for the distance. I was thinking of one dimension but in a plane the distance can be greater.
The binary search idea sounds good. This problem essentially is "how do I approximate an integral square root quickly." There are potentially hundreds of points that will need these calculations, and coupled with calculations for other data structures, I want to find a solution that may not be ideal but is good. I have not profiled the standard C sqrt() function yet, but just looking at the code, I can probably come up with something slightly faster given the constraints of the problem domain.
Most points will be relatively close to the origin (within 1,000 units), so I can most likely optimize it by separating the sample space into disjoint chunks and filtering [smiley=vartriangle.gif]x2 [smiley=vartriangle.gif]y2 into the proper one. Even if I separated it into 1024 unit wide partitions that would take at most 10 comparisons.
|
|
IP Logged |
x = (0x2B | ~0x2B)
x == the_question
|
|
|
towr
wu::riddles Moderator
Uberpuzzler
Some people are average, some are just mean.
Gender:
Posts: 13730
|
|
Re: Integer Square Root
« Reply #5 on: Apr 28th, 2004, 10:47am » |
Quote Modify
|
given d^2 = (x1-x2)^2+(y1-y2)^2
You can get a first guess for d by finding the position of the most significant bit, if it's at position i, take a=2^(i/2)
(takes 5 or 6 steps, depending on wether you have 32 or 64 bits)
then iterate
a = a+(d^2-a*a)/(2a)
untill (a+1)^2 > d^2 >= a^2
It should converge faster than binary search..
|
| « Last Edit: Apr 28th, 2004, 12:25pm by towr » |
IP Logged |
Wikipedia, Google, Mathworld, Integer sequence DB
|
|
|
towr
wu::riddles Moderator
Uberpuzzler
Some people are average, some are just mean.
Gender:
Posts: 13730
|
|
Re: Integer Square Root
« Reply #6 on: Apr 28th, 2004, 1:23pm » |
Quote Modify
|
The following seems to work (I've had the computer test it for every integer up to 2^31-1)
Code:int isqrt(int square)
{
int a=square, c, m=16;
if (a >> m) m+=8; else m-=8;
if (a >> m) m+=4; else m-=4;
if (a >> m) m+=2; else m-=2;
if (a >> m) m+=2;
a = 1 << (m/2);
while(c = (square-a*a)/(2*a))
a+=c;
return a-1;
}
|
|
|
| « Last Edit: Apr 28th, 2004, 1:29pm by towr » |
IP Logged |
Wikipedia, Google, Mathworld, Integer sequence DB
|
|
|
John_Gaughan
Uberpuzzler
Behold, the power of cheese!
Gender:
Posts: 767
|
|
Re: Integer Square Root
« Reply #7 on: Apr 28th, 2004, 5:42pm » |
Quote Modify
|
Thanks towr, that helps. Now I need to do some testing and integrate it into my class.
|
|
IP Logged |
x = (0x2B | ~0x2B)
x == the_question
|
|
|
Barukh
Uberpuzzler
Gender:
Posts: 2276
|
|
Re: Integer Square Root
« Reply #8 on: Apr 29th, 2004, 12:17am » |
Quote Modify
|
I am totally confused…  It seems you are overcomplicating things, John. Here are my considerations:
1. If the coordinates are 16 bits, then double precision (53 bits) is more than enough. So, calculating the square of the distance, converting it to double, taking a square and converting it back to integer will do the work.
2. If the coordinates are 32 bits, then you may use long double precision (64 bits – but you need to check your compiler supports this level of accuracy) and use sqrtl() function that works on long doubles. Then, the result you get may differ from the true one by at most 1.
3. As for the speed, I compared the sqrt() implementation with the towr’s isqrt() implementation (which, by the way, I found very elegant), and got the following numbers (running 100,000,000 iterations on Pentium4 machine, gcc –O3 compilation):
sqrt() : 1.84 sec
isqrt(): 28.59 sec
|
|
IP Logged |
|
|
|
John_Gaughan
Uberpuzzler
Behold, the power of cheese!
Gender:
Posts: 767
|
|
Re: Integer Square Root
« Reply #9 on: Apr 29th, 2004, 6:05am » |
Quote Modify
|
on Apr 29th, 2004, 12:17am, Barukh wrote:| 1. If the coordinates are 16 bits, then double precision (53 bits) is more than enough. So, calculating the square of the distance, converting it to double, taking a square and converting it back to integer will do the work. |
|
The exact size of the integers involved is not all that important. The types will be 16 bit integers, and I believe that all CPUs supporting MMX or SSE can use those registers for 64 bit integer arithmetic. That would be most CPUs in the last decade.
Quote:| 2. If the coordinates are 32 bits, then you may use long double precision (64 bits – but you need to check your compiler supports this level of accuracy) and use sqrtl() function that works on long doubles. Then, the result you get may differ from the true one by at most 1. |
|
Long double functions are non-standard and generally inefficient on older hardware. Not everyone has a Pentium 4, and my program needs to run well even on older hardware.
Quote:3. As for the speed, I compared the sqrt() implementation with the towr’s isqrt() implementation (which, by the way, I found very elegant), and got the following numbers (running 100,000,000 iterations on Pentium4 machine, gcc –O3 compilation):
sqrt() : 1.84 sec
isqrt(): 28.59 sec
|
|
I find this odd. I do not have the time now (or yesterday) to test it out myself. It may take some tweaking by hand. I will also use gcc (g++ actually) to compile, so the libc should be the same or close enough to it. I imagine sqrt() is not something that changes much anymore between libc versions. All I remember from looking at libc's implementation of sqrt() is that I could probably do it faster given the (integer) constraints of the problem domain.
The timing looks decent, given that I will probably only need to do 1,000 calculations at most. More are possible, but extremely rare.
The fact that you run that on a Pentium 4 may be a factor too. It has very good floating point arithmetic, better than its predecessors. Anyway, I need to do some testing this weekend when I have time and I will share what I come up with.
If nothing else, it was a good exercise in CS 
|
|
IP Logged |
x = (0x2B | ~0x2B)
x == the_question
|
|
|
Grimbal
wu::riddles Moderator
Uberpuzzler
Gender:
Posts: 7527
|
|
Re: Integer Square Root
« Reply #11 on: May 6th, 2004, 1:36pm » |
Quote Modify
|
Here is an integer square root that should be fast without being too complicated. Faster method would involve lookup tables.
unsigned usqrt32c(unsigned long n)
{
int i;
unsigned d = 1U<<15;
unsigned long d2 = 1UL<<30;
unsigned result;
unsigned long result2, trial;
result = 0;
result2 = 0;
for (i = 16; i != 0; i--) {
trial = result2 + ((unsigned long)result<<i) + d2;
if (trial <= n) {
result += d;
result2 = trial;
}
d >>= 1;
d2 >>= 2;
}
return result;
}
|
| « Last Edit: Jan 23rd, 2005, 3:48pm by Grimbal » |
IP Logged |
|
|
|
ChipBizGuy
Guest
|
I don't know if you still care, but I believe there is a better solution than just coding isqrt because you have additional information. Let Dx = abs(x1-x2) and Dy = abs(y1-y2). The minimum possible distance is max(Dx, Dy). Then maximum is Dx+Dy. Use these as the bounds of the binary search. This reduces complexity to log(min(Dx,Dy)) compares.
Mark
|
|
IP Logged |
|
|
|
|
RIDDLES SITE
WRITE MATH!
Home 
Help 
Search 
Members 
Login 
Register
Reply
Notify of replies
Send Topic
Print
Remove