The general idea of the algorithm is pretty straightforward, it basically boils down to the dichotomic search idea by cutting your problem in half at each step. But the parity of the length of the arrays introduces some noisy intricacies.
Notations
$X$ and $Y$ two sorted arrays of the same size (their starting index is 0).
$n = |X| = |Y|$ the arrays' size
$m = \bigg\lfloor \dfrac{n}{2} \bigg\rfloor$ the middle index
$x = X[m]$ and $y = Y[m]$ the middle values
Note that, since the two arrays are the same length, there is a unique pair of median values.
The basic idea
At each step of the algorithm, we are discarding the same number of values to the left and right of the merged sorted arrays : the medians are preserved. For example :
The mathematical proof would be cumbersome to write.
The algorithm
Let's assume that $n$ is odd. If :
$x = y$ then they must be our two medians. If we sorted the merged array, those values would necessarily appear next to each other in the middle of it.
$x < y$ the problem can be reduced. It is the same as to find the two medians of the arrays $X[m:n]$ and $Y[0:m + 1]$ (I use the Python's slice notation here, endpoint is excluded).
$x > y$ is the symmetrical case, we should consider the outer part $X[0:m+1]$ and $Y[m:n]$.
Note that at each step, $X$ and $Y$ are always of the same length. As you can see, we can recursively apply this process until the two arrays are of total size of 4 or 2, they are the base cases.
If $n$ is pair, we have to be careful about the next endpoints we consider as the two medians, as they aren't a true middle in the two arrays. Let be
$m^- = m - 1$, $~x^- = X[ m^- ]$ and $~y^- = Y[ m^- ]$. If :
$x = y$ then the two medians must be the two medians of $[x^-, x],~[y^-, y]$ : it is one of our base case.
$x < y$ then we should consider $X[m - 1:n]$ and $Y[0:m + 1]$
$x < y$ we should consider $X[0: m + 1]$ and $Y[m - 1: n]$
As you can see, we carry the two middle elements to the next iteration !
Naive first implementation
If we implement the above algorithm in Python :
def is_odd(n: int) -> bool:
return n & 1
def base_case_medians_2(a: int, b: int) -> tuple[int, int]:
""" Returns the two medians (ordered) of the 2 given values. """
return (a, b) if a < b else (b, a)
def base_case_medians_4(a: int, b: int, c: int, d: int) -> tuple[int, int]:
""" Returns the two medians (ordered) of the 4 given values. """
return ((c, d) if b >= d else (b, c) if b <= c else (c, b)) if a <= c else ((a, b) if b <= d else (a, d) if a <= d else (d, a))
def medians_recursive_naive(X: list, Y: list) -> tuple[int, int]:
""" Returns the two medians of the merged (assumed sorted and of the same size) given arrays.
Runs in quasilinear time as unecessary copies of the arrays are made. """
n = len(X)
# Recursion's base cases
if n == 1:
return base_case_medians_2(X[0], Y[0])
elif n == 2:
return base_case_medians_4(X[0], X[1], Y[0], Y[1])
m = n // 2
x, y = X[m], Y[m]
# Branch on the parity of the arrays' size
if is_odd(n):
if x == y:
return (x, y)
elif x < y:
return medians_recursive_naive(X[m:], Y[:m + 1])
else:
return medians_recursive_naive(X[:m + 1], Y[m:])
else:
if x == y:
return base_case_medians_4(X[m - 1], x, Y[m - 1], y)
elif x < y:
return medians_recursive_naive(X[m - 1:], Y[: m + 1])
else:
return medians_recursive_naive(X[: m + 1], Y[m - 1:])
As the code shown above, an iteration takes $O(n)$ as copies of the arrays are made when forwarded to the next recursion call, but copies are unecessary if we instead keep track of the start and end indices of the two subarrays we consider.
Achieving logarithmic complexity
So, if we use instead indices, our function signature would look like :
def medians_recursive(x_start, x_end, y_start, y_end):
But, this recursive algorithm could be simplified into an iterative algorithm :
def medians_iterative(X: list, Y: list) -> tuple[int, int]:
""" Returns the two medians of the merged (assumed sorted and of the same size) given arrays.
Runs in logarithmic time. """
x_start, x_end = 0, len(X) # note that the `end` is excluded
y_start, y_end = 0, len(Y)
# Algorithm's main loop, cut in half (almost) the arrays at each iteration
while x_end - x_start > 2:
x_mid = (x_start + x_end) // 2
y_mid = (y_start + y_end) // 2
x, y = X[x_mid], Y[y_mid]
odd = is_odd(x_end - x_start)
if x == y:
return (x, x) if odd else base_case_medians_4(X[x_mid - 1], x, Y[y_mid - 1], y)
elif x < y:
x_start = x_mid + odd - 1
y_end = y_mid + 1
else:
x_end = x_mid + 1
y_start = y_mid + odd - 1
# Base cases
if x_end - x_start == 1:
return base_case_medians_2(X[x_start], Y[y_start])
else:
return base_case_medians_4(X[x_start], X[x_start + 1], Y[y_start], Y[y_start + 1])
