# python – count points in octant of a square of size M closer in angle to the x-axis

In response to easy function from a pair of 32-bit ints to a single 64-bit int that preserves rotational order by Aubrey da Cunha, Mark Dickinson presented a pair of mutually recursive functions instrumental to that end.
The almost repetition irritates me: Below the almost original code and a stab at unification for comparative review (Reaction of Mark Dickinson pending).
Any insufficient understanding in “module” docstring and additional variable name is mine.

``````""" Limits on count of points closer in angle to the x-axis than (p, q).
Mark Dickinson in https://stackoverflow.com/a/66910195/ 2021/04/01
"""

def floor_sum(p, q, M):
"""
Sum of floor(q * x / p) for 0 <= x < M.

Assumes p positive, q and M nonnegative.
"""
a = q // p
r = q % p
scaledSlope = a * M * (M - 1) // 2
if r == 0:
return scaledSlope
N = (M * r + p - 1) // p
return scaledSlope + (N - 1) * M - ceil_sum(r, p, N)

def ceil_sum(p, q, M):
"""
Sum of ceil(q * x / p) for 0 <= x < M.

Assumes p positive, q and M nonnegative.
"""
a = q // p
r = q % p
#    M_ = M - 1
scaledSlope = a * M * (M - 1) // 2
if r == 0:
return scaledSlope
N = (M * r + p - 1) // p
return scaledSlope + N * (M - 1) - floor_sum(r, p, N) ~~~
``````

versus

``````""" Limits on count of points closer in angle to the x-axis than (p, q).
Mark Dickinson in https://stackoverflow.com/a/66910195/ 2021/04/01
"""

def _sum(p, q, M, sum_ceil):
"""
Sum of ceil(q * x / p) for 0 <= x < M
- or floor, if not bool sum_ceil.
Assumes p positive, q and M nonnegative.
"""
decrement_M, decrement_N = 1, 0 if sum_ceil else 0, 1
sum_ = 0
MM_ = M * (M - 1)
while True:
a, r = divmod(q, p)
sum += a * MM_ // 2
if r == 0:
return sum_
N = (M * r + p - 1) // p
sum += (N - decrement_N) * (M - decrement_M)
decrement_M, decrement_N = decrement_N, decrement_M

def floor_sum(p, q, M):
"""
Sum of floor(q * x / p) for 0 <= x < M.

Assumes p positive, q and M nonnegative.
"""
return _sum(p, q, M, False)

def ceil_sum(p, q, M):
"""
Sum of ceil(q * x / p) for 0 <= x < M.

Assumes p positive, q and M nonnegative.
"""
return _sum(p, q, M, ceil_sum)
``````

For the sake of completeness, the mapping function using both:

``````""" Bijection from (0, M) x (0, M) to (0, M^2),
preserving the 'angle' ordering.
Mark Dickinson in https://stackoverflow.com/a/66910195/ 2021/04/01
"""
from math import gcd

def point_to_line(p, q, M):
"""
Bijection from (0, M) x (0, M) to (0, M^2), preserving
the 'angle' ordering.
"""
if p == q == 0:
return 0
elif q <= p:
return ceil_sum(p, q, M) + gcd(p, q)
else:
return M * (M - 1) - floor_sum(q, p, M) + gcd(p, q)
``````