Calculating the modulo of two intervals

Question

I want to understand how the modulus operator works when applied to two intervals. Adding, subtracting and multiplying two intervals is trivial to implement in code, but how do you do it for modulus?

I'd be happy if someone can show me the formula, sample code or a link which explains how it works.

Background info: You have two integers x_lo < x < x_hi and y_lo < y < y_hi. What is the the lower and upper bound for mod(x, y)?

Edit: I'm unsure if it is possible to come up with the minimal bounds in an efficient manner (without calculating the mod for all x or for all y). If so, then I'll accept an accurate but non-optimal answer for the bounds. Obviously, [-inf,+inf] is a correct answer then :) but I want a bound that is more limited in size.

Answer 1

Let see mod(x, y) = mod. In general 0 <= mod <= y. So it's always true: y_lo < mod < y_hi

But we can see some specific cases below:

- if: x_hi < y_lo then div(x, y) = 0, then x_low < mod < x_hi
- if: x_low > y_hi then div(x, y) > 0, then y_low < mod < y_hi
- if: x_low < y_low < y_hi < x_hi, then y_low < mod < y_hi
- if: x_low < y_low < x_hi < y_hi, then y_low < mod < x_hi
- if: y_low < x_low < y_hi < x_hi, then y_low < mod < y_hi

....

Answer 2

It turns out, this is an interesting problem. The assumption I make is that for integer intervals, modulo is defined with respect to truncated division (round towards 0).

As a consequence, mod(-a,m) == -mod(a,m) for all a, m. Moreover, sign(mod(a,m)) == sign(a).

Definitions, before we start

Closed interval from a to b: [a,b]
Empty interval: [] := [+Inf,-Inf]
Negation: -[a,b] := [-b,-a]
Union: [a,b] u [c,d] := [min(a,c),max(b,d)]
Absolute value: |m| := max(m,-m)

Simpler Case: Fixed modulus m

It is easier to start with a fixed m. We will later generalize this to the modulo of two intervals. The definition builds up recursively. It should be no problem to implement this in your favorite programming language. Pseudocode:

def mod1([a,b], m):
    // (1): empty interval
    if a > b || m == 0:
        return []
    // (2): compute modulo with positive interval and negate
    else if b < 0:
        return -mod1([-b,-a], m)
    // (3): split into negative and non-negative interval, compute and join 
    else if a < 0:
        return mod1([a,-1], m) u mod1([0,b], m)
    // (4): there is no k > 0 such that a < k*m <= b
    else if b-a < |m| && a % m <= b % m:
        return [a % m, b % m]
    // (5): we can't do better than that
    else
        return [0,|m|-1]

Up to this point, we can't do better than that. The resulting interval in (5) might be an over-approximation, but it is the best we can get. If we were allowed to return a set of intervals, we could be more precise.

General case

The same ideas apply to the case where our modulus is an interval itself. Here we go:

def mod2([a,b], [m,n]):
    // (1): empty interval
    if a > b || m > n:
        return []
    // (2): compute modulo with positive interval and negate
    else if b < 0:
        return -mod2([-b,-a], [m,n])
    // (3): split into negative and non-negative interval, compute, and join 
    else if a < 0:
        return mod2([a,-1], [m,n]) u mod2([0,b], [m,n])
    // (4): use the simpler function from before
    else if m == n:
        return mod1([a,b], m)
    // (5): use only non-negative m and n
    else if n <= 0:
        return mod2([a,b], [-n,-m])
    // (6): similar to (5), make modulus non-negative
    else if m <= 0:
        return mod2([a,b], [1, max(-m,n)])
    // (7): compare to (4) in mod1, check b-a < |modulus|
    else if b-a >= n:
        return [0,n-1]
    // (8): similar to (7), split interval, compute, and join
    else if b-a >= m:
        return [0, b-a-1] u mod2([a,b], [b-a+1,n])
    // (9): modulo has no effect
    else if m > b:
        return [a,b]
    // (10): there is some overlapping of [a,b] and [n,m]
    else if n > b:
        return [0,b]
    // (11): either compute all possibilities and join, or be imprecise
    else:
        return [0,n-1] // imprecise

Have fun! :)

Answer 3

I'll assume both intervals are positive, in other words 0 <= x_lo and 0 < y_lo — only because negative numbers make things more complicated.

Let's look at the definition of mod:

x = k*y + mod(x, y)
k = floor(x/y)

Division isn't too complicated when applied to non-negative intervals: floor(x_lo / y_hi) <= k <= floor(x_hi / y_lo). Call those bounds k_lo and k_hi. If we calculate that k_lo = k_hi = k, then we can substitute back into the equation above and solve for mod(x, y).

mod(x, y) = x - k*y

So x_lo - k*y_hi <= mod(x, y) <= x_hi - k*y_lo when k_lo = k_hi.

---

If you look at a 3d graph of mod(x, y) alongside k = floor(x/y), then you'll see large regions that share a value of k, which look linear (flat planes) in x and y. The calculation above is for finding the bounds when the result is confined to a single one of those regions, where k has a particular value. The linearity makes it easy to calculate the bounds. These 3d regions remind me of how a 2d graph of mod(x, c) for a constant c has a sawtooth pattern, divided into linear regions. Next we'll look at when the modulo result can span between multiple regions.

In all other cases where k_lo < k_hi, the possible results of span across multiple values of k. Each of the regions is a plane/linear, so if we're looking for upper and lower bounds, they'll be along the boundaries of regions or the edges of the intervals. (I'm sure there's a justification in terms of derivatives and discontinuities.) Take everything that follows with a grain of salt.

Those boundaries between values of k = floor(x/y) are when x/y is an integer, that is whenever x is a multiple of y. If x and y were real numbers, the lines where x is a multiple of y would be where mod(x, y) approaches 0 from above, and y from below. We know the intervals for x and y contain a value on both sides of one of these boundaries since k_lo < k_hi, hence lower_bound = 0.

I think there are two ways for a point to be the maximum output of mod(x, y) within a region. A k region boundary intersected with either the x = x_hi plane or the y_hi plane. (If no k region boundary intersected either of them, then k_lo would equal k_hi.)

Setting x = x_hi as constant, then the maximum value of mod(x_hi, y) is at the smallest y such that k = floor(y / x_hi) is at its maximum (i.e. k = k_hi). I'm going by the graph here. So one possible upper bound is mod(x_hi, clamp(k_hi*x_hi, y_lo, y_hi)). Unfortunately, my calculations don't lead to the right answer here, so you might need to figure this part out and test.

Setting y = y_hi as constant, then the maximum will depend on whether the plane spans over a k region boundary, i.e. whether or not floor(x_lo / y_hi) < floor(x_hi / y_hi).

• If so, then the maximum for mod(x, y_hi) is just y_hi - 1. (Justification: if k = floor(x/y_hi) for some x, but k-1 = floor((x-1)/y_hi), then x is congruent to 0 mod y_hi. x-1 is congruent to 0-1 mod y_hi, in other words mod(x-1, y_hi) = y_hi - 1 which gives the upper bound.)
• If the plane doesn't span a k region boundary, then the maximum will just be mod(x_hi, y_hi) (simple because within each k region it's linear in x). I think this case will never matter anyway, because the maximum will be smaller than the result of setting x = x_hi anyway, but I'm not certain enough to leave it out.

That gives us the other upper bound.

If my maths is correct, we can just take the max of these two upper bounds for the case when k_lo < k_hi.

---

So all in all, it would look something like:

def mod_intervals(xl, xh, yl, yh):
    kl = xl // yh
    kh = xh // yl
    if kl == kh:
        k = kl
        return (xl - k*yh, xh - k*yl)
    else:
        bound_at_xh = xh % min(yh, max(yl, kh*xh)) # this line is likely to be incorrect
        if xl // yh < xh // yh:
            bound_at_yh = yh - 1
        else:
            bound_at_yh = xh % yh
        return (0, max(bound_at_xh, bound_at_yh))