mindahl
mindahl

Reputation: 77

What is a faster solution for integer coordinates in circle of radius x

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

Answers (3)

B. M.
B. M.

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.

  • For algorithm, I suggest you a better way : enter image description here

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

willeM_ Van Onsem
willeM_ Van Onsem

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.

Performance

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

ShlomiF
ShlomiF

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

Related Questions