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)