Ancestro común más bajo - algoritmo de Farach-Colton y Bender

Sea $G$ un árbol. Para cada consulta de la forma $(u, v)$ queremos encontrar el ancestro común más bajo de los nodos $u$ y $v$, es decir, queremos encontrar un nodo $w$ que yace en el camino de $u$ al nodo raíz, que yace en el camino de $v$ al nodo raíz, y si hay varios nodos elegimos el que está más lejos del nodo raíz. En otras palabras, el nodo deseado $w$ es el ancestro más bajo de $u$ y $v$. En particular, si $u$ es un ancestro de $v$, entonces $u$ es su ancestro común más bajo.

El algoritmo que se describe en este artículo fue desarrollado por Farach-Colton y Bender. Es asintóticamente óptimo.

Algoritmo

Usamos la reducción clásica del problema de LCA al problema de RMQ. Recorremos todos los nodos del árbol con DFS y mantenemos un arreglo con todos los nodos visitados y las alturas de estos nodos. El LCA de dos nodos $u$ y $v$ es el nodo entre las ocurrencias de $u$ y $v$ en el tour que tiene la menor altura.

En la siguiente imagen se puede ver un posible tour de Euler de un grafo y en la lista de abajo se pueden ver los nodos visitados y sus alturas.

$$\overline{)\begin{array}{lccccccccccccc}\text{Nodos:}& 1& 2& 5& 2& 6& 2& 1& 3& 1& 4& 7& 4& 1\\ \text{Alturas:}& 1& 2& 3& 2& 3& 2& 1& 2& 1& 2& 3& 2& 1\end{array}}$$$$\begin{array}{|l|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline \text{Nodos:} & 1 & 2 & 5 & 2 & 6 & 2 & 1 & 3 & 1 & 4 & 7 & 4 & 1 \ \hline \text{Alturas:} & 1 & 2 & 3 & 2 & 3 & 2 & 1 & 2 & 1 & 2 & 3 & 2 & 1 \ \hline \end{array}$$

Se puede leer más sobre esta reducción en el artículo Ancestro común más bajo. En ese artículo el mínimo de un rango se encontraba o bien por descomposición por raíz cuadrada en $O(\sqrt{N})$ o bien en $O(\log N)$ usando un Árbol de Segmentos. En este artículo vemos cómo podemos resolver las consultas de mínimo en un rango dadas en tiempo $O(1)$, tomando todavía solo $O(N)$ de tiempo para el preprocesamiento.

Nótese que el problema de RMQ reducido es muy específico: cualquier par de elementos adyacentes del arreglo difieren exactamente en uno (ya que los elementos del arreglo no son más que las alturas de los nodos visitados en el orden del recorrido, y o bien vamos a un descendiente, en cuyo caso el siguiente elemento es uno mayor, o bien volvemos al ancestro, en cuyo caso el siguiente elemento es uno menor). El algoritmo de Farach-Colton y Bender describe una solución exactamente para este problema de RMQ especializado.

Denotemos con $A$ el arreglo sobre el que queremos realizar las consultas de mínimo en un rango. Y $N$ será el tamaño de $A$.

Hay una estructura de datos sencilla que podemos usar para resolver el problema de RMQ con preprocesamiento $O(N \log N)$ y $O(1)$ por cada consulta: la Tabla Dispersa. Creamos una tabla $T$ donde cada elemento $T[i][j]$ es igual al mínimo de $A$ en el intervalo $[i, i + 2^j - 1]$. Obviamente $0 \leq j \leq \lceil \log N \rceil$, y por tanto el tamaño de la Tabla Dispersa será $O(N \log N)$. Se puede construir la tabla fácilmente en $O(N \log N)$ observando que $T[i][j] = \min(T[i][j-1], T[i+2^{j-1}][j-1])$.

¿Cómo podemos responder una consulta RMQ en $O(1)$ usando esta estructura de datos? Sea la consulta recibida $[l, r]$, entonces la respuesta es $\min(T[l][\text{sz}], T[r-2^{\text{sz}}+1][\text{sz}])$, donde $\text{sz}$ es el mayor exponente tal que $2^{\text{sz}}$ no es mayor que la longitud del rango $r-l+1$. En efecto, podemos tomar el rango $[l, r]$ y cubrirlo con dos segmentos de longitud $2^{\text{sz}}$: uno que empieza en $l$ y el otro que termina en $r$. Estos segmentos se solapan, pero esto no interfiere con nuestro cálculo. Para lograr realmente la complejidad temporal de $O(1)$ por consulta, necesitamos conocer los valores de $\text{sz}$ para todas las longitudes posibles de $1$ a $N$. Pero esto se puede precomputar fácilmente.

Ahora queremos mejorar la complejidad del preprocesamiento hasta $O(N)$.

Dividimos el arreglo $A$ en bloques de tamaño $K = 0.5 \log N$, donde $\log$ es el logaritmo en base 2. Para cada bloque calculamos el elemento mínimo y lo guardamos en un arreglo $B$. $B$ tiene tamaño $\frac{N}{K}$. Construimos una tabla dispersa a partir del arreglo $B$. Su tamaño y su complejidad temporal serán:

$$\frac{N}{K}\log\left(\frac{N}{K}\right) = \frac{2N}{\log(N)} \log\left(\frac{2N}{\log(N)}\right) =$$

$$= \frac{2N}{\log(N)} \left(1 + \log\left(\frac{N}{\log(N)}\right)\right) \leq \frac{2N}{\log(N)} + 2N = O(N)$$

Ahora solo tenemos que aprender a responder rápidamente consultas de mínimo en un rango dentro de cada bloque. De hecho, si la consulta de mínimo en un rango recibida es $[l, r]$ y $l$ y $r$ están en bloques distintos, entonces la respuesta es el mínimo de los siguientes tres valores: el mínimo del sufijo del bloque de $l$ que empieza en $l$, el mínimo del prefijo del bloque de $r$ que termina en $r$, y el mínimo de los bloques entre ellos. El mínimo de los bloques intermedios se puede responder en $O(1)$ usando la Tabla Dispersa. Así que solo nos quedan las consultas de mínimo en un rango dentro de los bloques.

Aquí aprovecharemos la propiedad del arreglo. Recordemos que los valores del arreglo — que no son más que valores de altura en el árbol — siempre diferirán en uno. Si quitamos el primer elemento de un bloque y lo restamos de cada otro elemento del bloque, cada bloque se puede identificar por una secuencia de longitud $K - 1$ formada por los números $+1$ y $-1$. Como estos bloques son tan pequeños, solo hay unas pocas secuencias distintas que pueden ocurrir. El número de secuencias posibles es:

$$2^{K-1} = 2^{0.5 \log(N) - 1} = 0.5 \left(2^{\log(N)}\right)^{0.5} = 0.5 \sqrt{N}$$

Así, el número de bloques distintos es $O(\sqrt{N})$, y por tanto podemos precomputar los resultados de las consultas de mínimo en un rango dentro de todos los bloques distintos en tiempo $O(\sqrt{N} K^2) = O(\sqrt{N} \log^2(N)) = O(N)$. Para la implementación podemos caracterizar un bloque por una máscara de bits de longitud $K-1$ (que cabrá en un int estándar) y guardar el índice del mínimo en un arreglo $\text{block}[\text{mask}][l][r]$ de tamaño $O(\sqrt{N} \log^2(N))$.

Así aprendimos a precomputar consultas de mínimo en un rango dentro de cada bloque, así como consultas de mínimo en un rango sobre un rango de bloques, todo en $O(N)$. Con estas precomputaciones podemos responder cada consulta en $O(1)$, usando a lo sumo cuatro valores precomputados: el mínimo del bloque que contiene l, el mínimo del bloque que contiene r, y los dos mínimos de los segmentos solapados de los bloques entre ellos.

Implementación

int n;
vector<vector<int>> adj;
 
int block_size, block_cnt;
vector<int> first_visit;
vector<int> euler_tour;
vector<int> height;
vector<int> log_2;
vector<vector<int>> st;
vector<vector<vector<int>>> blocks;
vector<int> block_mask;
 
void dfs(int v, int p, int h) {
    first_visit[v] = euler_tour.size();
    euler_tour.push_back(v);
    height[v] = h;
    
    for (int u : adj[v]) {
        if (u == p)
            continue;
        dfs(u, v, h + 1);
        euler_tour.push_back(v);
    }
}
 
int min_by_h(int i, int j) {
    return height[euler_tour[i]] < height[euler_tour[j]] ? i : j;
}
 
void precompute_lca(int root) {
    // obtener el tour de Euler y los índices de las primeras ocurrencias
    first_visit.assign(n, -1);
    height.assign(n, 0);
    euler_tour.reserve(2 * n);
    dfs(root, -1, 0);
 
    // precomputar todos los valores de log
    int m = euler_tour.size();
    log_2.reserve(m + 1);
    log_2.push_back(-1);
    for (int i = 1; i <= m; i++)
        log_2.push_back(log_2[i / 2] + 1);
 
    block_size = max(1, log_2[m] / 2);
    block_cnt = (m + block_size - 1) / block_size;
 
    // precomputar el mínimo de cada bloque y construir la tabla dispersa
    st.assign(block_cnt, vector<int>(log_2[block_cnt] + 1));
    for (int i = 0, j = 0, b = 0; i < m; i++, j++) {
        if (j == block_size)
            j = 0, b++;
        if (j == 0 || min_by_h(i, st[b][0]) == i)
            st[b][0] = i;
    }
    for (int l = 1; l <= log_2[block_cnt]; l++) {
        for (int i = 0; i < block_cnt; i++) {
            int ni = i + (1 << (l - 1));
            if (ni >= block_cnt)
                st[i][l] = st[i][l-1];
            else
                st[i][l] = min_by_h(st[i][l-1], st[ni][l-1]);
        }
    }
 
    // precomputar la máscara de cada bloque
    block_mask.assign(block_cnt, 0);
    for (int i = 0, j = 0, b = 0; i < m; i++, j++) {
        if (j == block_size)
            j = 0, b++;
        if (j > 0 && (i >= m || min_by_h(i - 1, i) == i - 1))
            block_mask[b] += 1 << (j - 1);
    }
 
    // precomputar RMQ para cada bloque único
    int possibilities = 1 << (block_size - 1);
    blocks.resize(possibilities);
    for (int b = 0; b < block_cnt; b++) {
        int mask = block_mask[b];
        if (!blocks[mask].empty())
            continue;
        blocks[mask].assign(block_size, vector<int>(block_size));
        for (int l = 0; l < block_size; l++) {
            blocks[mask][l][l] = l;
            for (int r = l + 1; r < block_size; r++) {
                blocks[mask][l][r] = blocks[mask][l][r - 1];
                if (b * block_size + r < m)
                    blocks[mask][l][r] = min_by_h(b * block_size + blocks[mask][l][r], 
                            b * block_size + r) - b * block_size;
            }
        }
    }
}
 
int lca_in_block(int b, int l, int r) {
    return blocks[block_mask[b]][l][r] + b * block_size;
}
 
int lca(int v, int u) {
    int l = first_visit[v];
    int r = first_visit[u];
    if (l > r)
        swap(l, r);
    int bl = l / block_size;
    int br = r / block_size;
    if (bl == br)
        return euler_tour[lca_in_block(bl, l % block_size, r % block_size)];
    int ans1 = lca_in_block(bl, l % block_size, block_size - 1);
    int ans2 = lca_in_block(br, 0, r % block_size);
    int ans = min_by_h(ans1, ans2);
    if (bl + 1 < br) {
        int l = log_2[br - bl - 1];
        int ans3 = st[bl+1][l];
        int ans4 = st[br - (1 << l)][l];
        ans = min_by_h(ans, min_by_h(ans3, ans4));
    }
    return euler_tour[ans];
}