Menor Ancestral Comum - Farach-Colton e Bender

Deixe $G$ ser uma árvore. Para cada consulta da forma $(u, v)$ queremos encontrar o menor ancestral comum dos nós $u$ e $v$, ou seja, queremos encontrar um nó $w$ que se encontra no caminho de $u$ para o nó raiz e que se encontra no caminho de $v$ para o nó raiz e, se houver vários nós, escolheremos o que estiver mais distante do nó raiz. Em outras palavras, o nó desejado $w$ é o menor ancestral de $u$ e $v$. Em particular, se $u$ for um ancestral de $v$, então $u$ é o menor ancestral comum.

O algoritmo que será descrito neste artigo foi desenvolvido por Farach-Colton e Bender. É assintoticamente otimizado.

Algoritmo

Usamos a redução clássica do problema do LCA para o do RMQ. Atravessamos todos os nós da árvore com uma DFS e mantemos uma array com todos os nós visitados e as alturas desses nós. O LCA de dois nós $u$ e $v$ é o nó entre as ocorrências de $u$ e $v$ no tour, que tem a menor altura.

Na figura a seguir, é possível ver um possível tour de Euler de um grafo e, na lista abaixo, os nós visitados e suas alturas.

LCA_Euler_Tour
$$\begin{array}{|l|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline \text{Nos:} & 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}$$

Você pode ler mais sobre essa redução no artigo Menor Ancestral Comum. Nesse artigo, o mínimo de um intervalo foi encontrado por "sqrt-decomposition" em $O(\sqrt{N})$ ou em $O(\log N)$ usando uma Árvore Segmentária. Neste artigo, veremos como podemos resolver as consultas de intervalo mínimo(RMQ) em $O(1)$, enquanto ainda gastamos apenas $O(N)$ para pré-processamento.

Observe que o problema da RMQ reduzido é muito específico: quaisquer dois elementos adjacentes na matriz diferem exatamente por um (uma vez que os elementos da array nada mais são do que as alturas dos nós visitados na ordem da travessia, e nós vamos para um descendente, nesse caso o próximo elemento é um maior ou então, retornamos ao ancestral, nesse caso, o próximo elemento será um menor). O algoritmo de Farach-Colton e Bender descreve uma solução exatamente para esse problema específico da RMQ.

Vamos definir $A$ como a array na qual queremos executar as consultas de intervalo mínimo. E $N$ será o tamanho de $A$.

Existe uma estrutura de dados que podemos usar para resolver o problema do RMQ com pré-processamento $O(N \log N)$ e $O(1)$ para cada consulta: Sparse Table. Criamos uma tabela $T$ em que cada elemento $T[i][j]$ é igual ao mínimo de $A$ no intervalo $[i, i + 2^j - 1]$. Claramente, $0 \leq j \leq \lceil \log N \rceil$, e, portanto, o tamanho da Sparse Table será $O(N \log N)$. Você pode criar a tabela em $O(N \log N)$ notando que: $T[i][j] = \min(T[i][j-1], T[i+2^{j-1}][j-1])$.

Como podemos responder a uma consulta RMQ em $O(1)$ usando essa estrutura de dados? Seja a consulta recebida: $[l,r]$, e a resposta é $\min(T[l][\text{sz}], T[r-2^{\text{sz}}+1][\text{sz}])$, onde $\text{sz}$ é o maior expoente em que $2^{\text{sz}}$ não seja maior que o comprimento do intervalo $r-l+1$. De fato, podemos pegar o intervalo $[l, r]$ e cobrir dois segmentos de comprimento $2^{\text{sz}}$ - um começando em $l$ e o outro terminando em $r$. Esses segmentos se sobrepõem, mas isso não interfere em nosso cálculo. Para realmente atingir a complexidade de tempo $O(1)$ por consulta, precisamos conhecer os valores de $\text{sz}$ para todos os comprimentos possíveis de $1$ até $N$. Isso pode ser pré-computado

Agora queremos melhorar a complexidade do pré-processamento para $O(N)$.

Dividimos a array $A$ em blocos de tamanho $K = 0.5 \log N$ com $\log$ sendo o logaritmo de base 2. Para cada bloco, calculamos o elemento mínimo e os armazenamos em uma array $B$. $B$ tem tamanho $\frac{N}{K}$. Construímos uma sparse table da array $B$. O tamanho e a complexidade do tempo serão: $$\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)$$

Agora só precisamos aprender a responder rapidamente às consultas de intervalo mínimo em cada bloco. De fato, se a consulta de intervalo mínimo recebido for $[l,r]$ e $l$ e $r$ estiverem em blocos diferentes, a resposta será o mínimo dos três valores a seguir: o mínimo do sufixo do bloco de $l$ começando em $l$, o mínimo do prefixo do bloco de $r$ terminando em $r$, e o mínimo dos blocos entre eles. O mínimo dos blocos intermediários pode ser respondido em $O(1)$ usando a Sparse Table. Portanto, isso nos deixa apenas com as consultas de intervalo mínimo dentro dos blocos.

Aqui vamos explorar a propriedade da array. Lembre-se de que os valores na array - que são apenas valores de alturas na árvore - sempre diferem em um. Se removermos o primeiro elemento de um bloco e subtraí-lo de todos os outros itens do bloco, cada bloco poderá ser identificado por uma sequência de comprimento $K - 1$ consistindo do número $+1$ e $-1$. Como esses blocos são muito pequenos, existem apenas algumas sequências diferentes que podem ocorrer. O número de sequências possíveis é: $$2^{K-1} = 2^{0.5 \log(N) - 1} = 0.5 \left(2^{\log(N)}\right)^{0.5} = 0.5 \sqrt{N}$$

Portanto, o número de blocos diferentes é $O(\sqrt{N})$, e, portanto, podemos pré-calcular os resultados de consultas de intervalo mínimo dentro de todos os blocos diferentes em $O(\sqrt{N} K^2) = O(\sqrt{N} \log^2(N)) = O(N)$. Para a implementação, podemos caracterizar um bloco com uma máscara de bits de comprimento $K-1$ (que caberá em um int padrão) e armazenar o índice do mínimo em uma array $\text{block}[\text{mask}][l][r]$ de tamanho $O(\sqrt{N} \log^2(N))$.

Assim, poderemos pré-calcular consultas de intervalo mínimo em cada bloco, bem como consultas de intervalo mínimo sobre um intervalo de blocos, tudo em $O(N)$. Com essas pré-computações, podemos responder a cada consulta em $O(1)$, usando no máximo quatro valores pré-computados: o mínimo do bloco contendo l, o mínimo do bloco contendo r, e os dois mínimos dos segmentos sobrepostos dos blocos entre eles.

Implementação

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) {
    // pegar o tour de euler & os índices das primeiras ocorrências
    first_visit.assign(n, -1);
    height.assign(n, 0);
    euler_tour.reserve(2 * n);
    dfs(root, -1, 0);

    // precomputar todos os 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 o mínimo de cada bloco e construir a sparse table
    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 a mask(máscara de bits) para cada bloco
    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 bloco ú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];
}