M
x N
.X
, Y
]r
We can only move horizontally or vertically by 1 space (we always have to move). The goal is to find the total number of possibilities where we can be.
In this example we have M=5
, N=4
, [X=2
, Y=1
] and r=3
.
We can see, that r
is an odd number, so we're searching for odd numbers.
The result is obviously (after summing all the odd numbers) 12
I first tried brute force but it was too slow (O(n^2)
).
So I tried it more mathematically. I tried the Von Neumann Neighborhood algorithm, but I got stuck there totally. Finally, I'm trying it by calculating the vertical and horizontal counts in the row corresponding to x
and y
(count in x
is 3
, count in y
is 3
).
My next step is to find the corners and using Taxicab geometry calculate how many r
s it would take to go there.
Then I divide the image into 4 parts (quadrants). When the path from [X, Y]
to the corner of a given quadrant is less than r
, just divide the number of squares by 2. When the path is larger, I create r_
, which is determined by the difference of path to corner with r → r_ = taxicab - r
. Then just subtract r_
from the rectangle's contents and divide by 2 again.
As we can see, the result also comes out correctly, 3 + 3 + 1 + 1 + 2 + 2 = 12
.
BUT
Here I ask:
We are still working with even and odd numbers so that rounding errors often occur after dividing an odd number by 2. (for example for M=2, N=4, X=0, Y=0, r=4
- see picture. There are 8 fields - 3 blank. But 5/2 is 2.5 => why take 3 and not 2? I tried adding different rules, but it varied from example to example.)
Is there a more efficient algorithm for such a calculation that would be similarly fast and less error prone?
We can count the squares which lie outside of our grid, and subtract those from the total number of odd/even squares in range (regardless of in/out of bounds).
This is actually easier to do than we may think. First, we want to be able to count the number of odd/even squares in range, with no regards for boundaries. Using the formula for the sum of natural numbers up to N
- S = N * (N + 1) / 2
we can derive a couple of simple equations for this:
def count_all(size):
if size % 2: # odd size
return (size + 1) ** 2
return size * (size + 2) + 1
It may be a good exercise to try deriving these yourself, or at least verifying with a few examples that they are in fact correct.
Moving on- we can eliminate the points that fall out of bounds from above by "shrinking" our radius. This is very visual, so let me give you a diagram. Imagine only the points where the range is some fixed number, say range=13
, and our center is something in the lower-right quadrant, say (17, 5)
. If we plot these points, connecting them with lines, it creates a diamond:
| /\
| / \
| / \
-+-----------
| / \
| \ /
| \ /
| \ /
| \/
If we're only concerned with counting the points above the axis, we could equivalently just count the points above the axis of a smaller diamond that's shifted upwards accordingly. Example:
| /\
| / \
| / \
-+-----------
| \ /
| \ /
| \/
Now, this is very convenient to work with, because exactly half of the diamond is above, half is below. We do half to be careful though- there are points that fall on the axis, and both need to either be considered in bounds or out of bounds equivalently, but we can easily account for that.
Using this insight, we can count the number of points that fall out of bounds across an axis by shrinking the range and shifting the center point, and counting points on half of this new plot. The counting code:
def count_side(size):
if size % 2:
return (size + 1) // 2
return size // 2 + 1
def count_half(size):
if size < 0:
return 0
return count_all(size) // 2 + count_side(size)
Note that we have to be careful for even ranges, since we need to count the center (range 0) exactly once.
We aren't done yet though- if we just subtract out the number of points that are out of bounds above and then to the left independently, we're overcount the number of points to remove, since we count points in the top-left quadrant twice. To handle this, we use the same trick. We'll shrink + shift the diamond across the x-axis first, and then we'll do it again on this new diamond, but across the y-axis. Note that this new diamond will end up centered on the origin. I'd recommend that you pause a moment at this point to visualize this and convince yourself that this is fine, and will in fact give us a new diamond for any specific range, containing 4x the number of points which fall in the top-left quadrant.
Using this, we count the number of points in the top-left quadrant, re-add them to the total. We then repeat this same process for the right side, bottom, and other three corners, to obtain the overall total. The entire solution below:
from itertools import product
def count_all(size):
if size % 2:
return (size + 1) ** 2
return size * (size + 2) + 1
def count_side(size):
if size % 2:
return (size + 1) // 2
return size // 2 + 1
def count_half(size):
if size < 0:
return 0
return count_all(size) // 2 + count_side(size)
def count_quarter(size):
if size < 0:
return 0
return count_all(size) // 4 + count_side(size)
def get_deltas(pos, limit):
return -(pos + 1), pos - (limit + 1)
def count_inside(c, r, x, y, s):
total = count_all(s)
vertical_deltas = get_deltas(x, c)
horizontal_deltas = get_deltas(y, r)
out_sides = sum(count_half(s + delta)
for delta in horizontal_deltas + vertical_deltas)
out_corners = sum(count_quarter(s + delta_vert + delta_horiz)
for delta_vert, delta_horiz in product(vertical_deltas, horizontal_deltas))
inside = total - out_sides + out_corners
return inside
A couple of examples:
>>> print(count_inside(5, 4, 2, 1, 3))
12
>>> print(count_inside(5, 4, 2, 1, 4))
14
>>> print(count_inside(5, 4, 2, 1, 5))
15
>>> print(count_inside(10, 6, 3, 2, 8))
36