Reputation: 77
I'm trying to create a method that returns the count of integer coordinates within a circle of radius rad but I believe my current solution is too slow. What would you recommend for a better implementation? I'd like to code the solution myself as I'm still a beginner.
This is my current solution:
def points(rad):
possiblePoints = []
negX = -rad
negY = -rad
#creating a list of all possible points from -rad,-rad to rad, rad
while(negX <= rad):
while(negY <= rad):
possiblePoints.append((negX,negY))
negY += 1
negY = -rad
negX += 1
count = 0
index=0
negX = -rad
negY = -rad
distance = 0
#counting the possible points for distances within rad
for i in possiblePoints:
for j in range(0,2):
distance += (possiblePoints[index][j])**2
distance = distance**(0.5)
if(distance<=rad):
count+=1
distance = 0
index += 1
return count
Upvotes: 1
Views: 1051
Reputation: 18628
For faster, the same logic but a quick numpy solution :
def count_in(radius): r=radius+1 X,Y=np.mgrid[-r:r,-r:r] return (X2 + Y2 <= radius**2).sum()
You will speed your algorithm by ~x30 generally. Detail explanations later.
Just follow the blue line, evaluating pink area. Have fun.
Explanations for radius=3 :
In [203]: X
Out[203]:
array([[-4, -4, -4, -4, -4, -4, -4, -4],
[-3, -3, -3, -3, -3, -3, -3, -3],
[-2, -2, -2, -2, -2, -2, -2, -2],
[-1, -1, -1, -1, -1, -1, -1, -1],
[ 0, 0, 0, 0, 0, 0, 0, 0],
[ 1, 1, 1, 1, 1, 1, 1, 1],
[ 2, 2, 2, 2, 2, 2, 2, 2],
[ 3, 3, 3, 3, 3, 3, 3, 3]])
In [204]: Y
Out[204]:
array([[-4, -3, -2, -1, 0, 1, 2, 3],
[-4, -3, -2, -1, 0, 1, 2, 3],
[-4, -3, -2, -1, 0, 1, 2, 3],
[-4, -3, -2, -1, 0, 1, 2, 3],
[-4, -3, -2, -1, 0, 1, 2, 3],
[-4, -3, -2, -1, 0, 1, 2, 3],
[-4, -3, -2, -1, 0, 1, 2, 3],
[-4, -3, -2, -1, 0, 1, 2, 3]])
In [205]: X*X+Y*Y
Out[205]:
array([[32, 25, 20, 17, 16, 17, 20, 25],
[25, 18, 13, 10, 9, 10, 13, 18],
[20, 13, 8, 5, 4, 5, 8, 13],
[17, 10, 5, 2, 1, 2, 5, 10],
[16, 9, 4, 1, 0, 1, 4, 9],
[17, 10, 5, 2, 1, 2, 5, 10],
[20, 13, 8, 5, 4, 5, 8, 13],
[25, 18, 13, 10, 9, 10, 13, 18]])
And finally :
In [207]: (X**2 + Y**2 <= 3**2)
Out[207]:
array([[False, False, False, False, False, False, False, False],
[False, False, False, False, True, False, False, False],
[False, False, True, True, True, True, True, False],
[False, False, True, True, True, True, True, False],
[False, True, True, True, True, True, True, True],
[False, False, True, True, True, True, True, False],
[False, False, True, True, True, True, True, False],
[False, False, False, False, True, False, False, False]], dtype=bool)
In [208]: (X**2 + Y**2 <= 3**2).sum() # True is one
Out[208]: 29
Upvotes: 2
Reputation: 476574
A point (x, y) is within the radius given x2+y2≤r2. We can calculate the number of coordinates for a fixed y-value as:
ny=1 + 2×⌊√(r2-y2)⌋
If we then slice y between 0 and r, and we count each ny for y>0 twice, we obtain the total number, so we can calculate this as:
import numpy as np
from math import floor
def total(radius):
r0 = floor(radius)
ys = np.arange(1, r0+1)
halves = 2 * np.sum(2*np.floor(np.sqrt(radius * radius - ys * ys)) + 1)
return halves + 1 + 2 * r0
We thus calculate for each "layer" the number of integral coordinates, and we double each layer, since there is a "co-layer" with the same number of integral coordinates in the opposite direction. We then add the number of coordinates with y=0 which is 1+2×⌊r⌋.
The above thus works in O(r) with r the radius of the circle.
Sample output:
>>> total(0)
1.0
>>> total(0.1)
1.0
>>> total(1)
5.0
>>> total(1.1)
5.0
>>> total(1.42)
9.0
>>> total(2)
13.0
>>> total(3)
29.0
A slower alternative in O(r2) is to generate a grid, and then perform a comparison in bulk on it, like:
import numpy as np
from math import floor
def total_slow(radius):
r0 = floor(radius)
xs = np.arange(-r0, r0+1)
ys = xs[:, None]
return np.sum(xs*xs+ys*ys <= radius*radius)
The above is however useful if we want to verify that our solution works. For example:
>>> total_slow(0)
1
>>> total_slow(1)
5
>>> total_slow(2)
13
>>> total_slow(3)
29
produces the same results, but we can use this to verify that it always produces the same result for a large number of radiuses, for example:
>>> [total(r) == total_slow(r) for r in np.arange(0, 10, 0.03)]
[True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True]
This thus means that we verified that the two generated the same result for 0
, 0.03
, 0.06
, etc. up to 10
. The above is of course not a formal proof of the correctness, but it gives us some empricial evidence.
We can test performance by using the timeit
module, and we tested the algorithms with three radiusses, each experiment was repeated 10'000 times:
>>> timeit(partial(total, 1.23), number=10000)
0.11901686200144468
>>> timeit(partial(total, 12.3), number=10000)
0.1255619800031127
>>> timeit(partial(total, 123), number=10000)
0.1318465179974737
>>> timeit(partial(total_slow, 1.23), number=10000)
0.11859546599953319
>>> timeit(partial(total_slow, 12.3), number=10000)
0.15540562200112618
>>> timeit(partial(total_slow, 123), number=10000)
1.3335393390007084
>>> timeit(partial(total_slow, 1.23), number=10000)
0.11859546599953319
>>> timeit(partial(total_slow, 12.3), number=10000)
0.15540562200112618
>>> timeit(partial(total_slow, 123), number=10000)
1.3335393390007084
>>> timeit(partial(points, 1), number=10000)
0.1152820099996461
>>> timeit(partial(points, 12), number=10000)
3.7115225509987795
>>> timeit(partial(points, 123), number=10000)
433.3227958699972
With points
the method of @ShlomiF, we here use integral numbers, since for a radius with a floating point number, this method will produce incorrect results.
We thus obtain the following table:
radius total total_slow points
1.23 0.119 0.119 0.115*
12.3 0.125 0.155 3.711*
123. 0.131 1.334 433.323*
* with integral radiuses
This is expected behavior if we take time complexity into account: eventually a linear approach will outperform a quadratic approach.
Upvotes: 3
Reputation: 2895
Things can definitely be made faster than using lots of unnecessary loops.
Why not "build" and "check" in the same loop?
Also, you can make this (almost) a one-liner with the help of product
from itertools:
import numpy as np
from itertools import product
def points(rad):
x = np.arange(-rad, rad + 1)
return len([(p1, p2) for p1, p2 in product(x, x) if p1 ** 2 + p2 ** 2 <= rad ** 2])
# Examples
print(points(5))
>> 81
print(points(1))
>> 5
print(points(2))
>> 13
The list comprehension gives the actual points, if that's of any use. Good luck!
Some explanations in case needed:
x
is just a list of integers between (and including) -rad
and rad
.
product(x, x)
is the Cartesian product of x and itself, thus all the points in the square defined by both coordinates being between -rad
and rad
.
The list comprehension returns only points with a distance from the origin smaller or equal to rad
.
The function returns the length of this list.
Upvotes: 1