Descomposición por raíz cuadrada (Sqrt Decomposition)

La descomposición por raíz cuadrada (Sqrt Decomposition) es un método (o una estructura de datos) que permite realizar algunas operaciones comunes (encontrar la suma de los elementos de un subarreglo, encontrar el elemento mínimo/máximo, etc.) en $O(\sqrt n)$ operaciones, lo cual es mucho más rápido que $O(n)$ del algoritmo trivial.

Primero describimos la estructura de datos para una de las aplicaciones más simples de esta idea, luego mostramos cómo generalizarla para resolver algunos otros problemas, y finalmente vemos un uso un poco distinto de esta idea: partir las consultas de entrada en bloques de tamaño raíz cuadrada.

Estructura de datos basada en descomposición por raíz cuadrada

Dado un arreglo $a[0 \dots n-1]$ , implementar una estructura de datos que permita encontrar la suma de los elementos $a[l \dots r]$ para $l$ y $r$ arbitrarios en $O(\sqrt n)$ operaciones.

Descripción

La idea básica de la descomposición por raíz cuadrada es el preprocesamiento. Dividiremos el arreglo $a$ en bloques de longitud aproximadamente $\sqrt n$ , y para cada bloque $i$ precalcularemos la suma de los elementos que contiene $b[i]$ .

Podemos asumir que tanto el tamaño del bloque como la cantidad de bloques son iguales a $\sqrt n$ redondeado hacia arriba:

$s = \lceil \sqrt n \rceil$

Entonces el arreglo $a$ se divide en bloques de la siguiente manera:

$\underbrace{a[0], a[1], \dots, a[s-1]}$

El último bloque puede tener menos elementos que los demás (si $n$ no es múltiplo de $s$ ); esto no es importante para la discusión (ya que se puede manejar fácilmente). Así, para cada bloque $k$ , conocemos la suma de los elementos en él $b[k]$ :

$b[k] = \sum\limits_{i=k\cdot s}^{\min {(n-1,(k+1)\cdot s - 1})} a[i]$

Así, hemos calculado los valores de $b[k]$ (esto requirió $O(n)$ operaciones). ¿Cómo nos ayudan a responder cada consulta $[l, r]$ ? Nótese que si el intervalo $[l, r]$ es suficientemente largo, contendrá varios bloques enteros, y para esos bloques podemos encontrar la suma de sus elementos en una sola operación. Como resultado, el intervalo $[l, r]$ contendrá partes de solo dos bloques, y tendremos que calcular la suma de los elementos en esas partes de forma trivial.

Así, para calcular la suma de los elementos en el intervalo $[l, r]$ solo necesitamos sumar los elementos de las dos “colas”: $[l\dots (k + 1)\cdot s-1]$ y $[p\cdot s\dots r]$ , y sumar los valores $b[i]$ en todos los bloques desde $k + 1$ hasta $p-1$ :

$\sum\limits_{i=l}^r a[i] = \sum\limits_{i=l}^{(k+1) \cdot s-1} a[i] + \sum\limits_{i=k+1}^{p-1} b[i] + \sum\limits_{i=p\cdot s}^r a[i]$

Nota: Cuando $k = p$ , es decir, $l$ y $r$ pertenecen al mismo bloque, la fórmula no se puede aplicar, y la suma se debe calcular de forma trivial.

Este enfoque nos permite reducir de forma significativa la cantidad de operaciones. En efecto, el tamaño de cada “cola” no supera la longitud del bloque $s$ , y la cantidad de bloques en la suma no supera $s$ . Como hemos elegido $s \approx \sqrt n$ , la cantidad total de operaciones necesarias para encontrar la suma de los elementos en el intervalo $[l, r]$ es $O(\sqrt n)$ .

Implementación

Empecemos con la implementación más simple:


// datos de entrada
int n;
vector<int> a (n);
 
// preprocesamiento
int len = (int) sqrt (n + .0) + 1; // tamaño del bloque y cantidad de bloques
vector<int> b (len);
for (int i=0; i<n; ++i)
    b[i / len] += a[i];
 
// responder las consultas
for (;;) {
    int l, r;
  // leer los datos de entrada de la siguiente consulta
    int sum = 0;
    for (int i=l; i<=r; )
        if (i % len == 0 && i + len - 1 <= r) {
            // si el bloque entero que empieza en i pertenece a [l, r]
            sum += b[i / len];
            i += len;
        }
        else {
            sum += a[i];
            ++i;
        }
}

Esta implementación tiene una cantidad excesiva de operaciones de división (que son mucho más lentas que otras operaciones aritméticas). En su lugar, podemos calcular los índices de los bloques $c_l$ y $c_r$ que contienen los índices $l$ y $r$ , e iterar sobre los bloques $c_l+1 \dots c_r-1$ con un procesamiento separado de las “colas” en los bloques $c_l$ y $c_r$ . Este enfoque corresponde a la última fórmula de la descripción, y convierte el caso $c_l = c_r$ en un caso especial.


int sum = 0;
int c_l = l / len,   c_r = r / len;
if (c_l == c_r)
    for (int i=l; i<=r; ++i)
        sum += a[i];
else {
    for (int i=l, end=(c_l+1)*len-1; i<=end; ++i)
        sum += a[i];
    for (int i=c_l+1; i<=c_r-1; ++i)
        sum += b[i];
    for (int i=c_r*len; i<=r; ++i)
        sum += a[i];
}

Otros problemas

Hasta ahora estábamos discutiendo el problema de encontrar la suma de los elementos de un subarreglo continuo. Este problema se puede extender para permitir actualizar elementos individuales del arreglo. Si un elemento $a[i]$ cambia, basta con actualizar el valor de $b[k]$ del bloque al que pertenece este elemento ( $k = i / s$ ) en una operación:

$b[k] += a_{new}[i] - a_{old}[i]$

Por otro lado, la tarea de encontrar la suma de los elementos se puede reemplazar por la de encontrar el elemento mínimo/máximo de un subarreglo. Si este problema también tiene que contemplar actualizaciones de elementos individuales, actualizar el valor de $b[k]$ también es posible, pero requerirá iterar sobre todos los valores del bloque $k$ en $O(s) = O(\sqrt{n})$ operaciones.

La descomposición por raíz cuadrada se puede aplicar de forma similar a toda una clase de otros problemas: encontrar la cantidad de elementos cero, encontrar el primer elemento distinto de cero, contar elementos que cumplen cierta propiedad, etc.

Otra clase de problemas aparece cuando hay que actualizar elementos del arreglo en intervalos: incrementar los elementos existentes o reemplazarlos por un valor dado.

Por ejemplo, digamos que podemos hacer dos tipos de operaciones sobre un arreglo: sumar un valor dado $\delta$ a todos los elementos del arreglo en el intervalo $[l, r]$ o consultar el valor del elemento $a[i]$ . Guardemos en $b[k]$ el valor que hay que sumar a todos los elementos del bloque $k$ (inicialmente todos los $b[k] = 0$ ). Durante cada operación “add” hay que sumar $\delta$ a $b[k]$ para todos los bloques que pertenecen al intervalo $[l, r]$ y sumar $\delta$ a $a[i]$ para todos los elementos que pertenecen a las “colas” del intervalo. La respuesta a la consulta $i$ es simplemente $a[i] + b[i/s]$ . De este modo la operación “add” tiene complejidad $O(\sqrt{n})$ , y responder una consulta tiene complejidad $O(1)$ .

Finalmente, esas dos clases de problemas se pueden combinar si la tarea requiere hacer ambas cosas: actualizaciones de elementos en un intervalo y consultas en un intervalo. Ambas operaciones se pueden hacer con complejidad $O(\sqrt{n})$ . Esto requerirá dos arreglos de bloques $b$ y $c$ : uno para llevar el seguimiento de las actualizaciones de elementos y otro para llevar el seguimiento de las respuestas a la consulta.

Existen otros problemas que se pueden resolver usando descomposición por raíz cuadrada; por ejemplo, un problema de mantener un conjunto de números que permita agregar/eliminar números, comprobar si un número pertenece al conjunto y encontrar el $k$ -ésimo número más grande. Para resolverlo hay que guardar los números en orden creciente, partidos en varios bloques con $\sqrt{n}$ números en cada uno. Cada vez que se agrega/elimina un número, hay que rebalancear los bloques moviendo números entre los comienzos y los finales de bloques adyacentes.

Algoritmo de Mo (Mo’s algorithm)

Una idea similar, basada en la descomposición por raíz cuadrada, se puede usar para responder consultas de rango ( $Q$ ) offline en $O((N+Q)\sqrt{N})$ . Esto puede sonar bastante peor que los métodos de la sección anterior, ya que es una complejidad un poco peor que la que teníamos antes y no se pueden actualizar valores entre dos consultas. Pero en muchas situaciones este método tiene ventajas. Durante una descomposición por raíz cuadrada normal, hay que precomputar las respuestas de cada bloque y fusionarlas al responder las consultas. En algunos problemas este paso de fusión puede ser bastante problemático. P. ej. cuando cada consulta pide encontrar la moda de su rango (el número que aparece con más frecuencia). Para esto cada bloque tendría que guardar la cantidad de cada número que contiene en algún tipo de estructura de datos, y ya no podemos realizar el paso de fusión lo suficientemente rápido. El algoritmo de Mo usa un enfoque completamente distinto, que puede responder este tipo de consultas de forma rápida, porque solo mantiene una estructura de datos, y las únicas operaciones sobre ella son fáciles y rápidas.

La idea es responder las consultas en un orden especial basado en los índices. Primero responderemos todas las consultas cuyo índice izquierdo está en el bloque 0, luego responderemos todas las consultas cuyo índice izquierdo está en el bloque 1, y así sucesivamente. Y además tendremos que responder las consultas de un bloque en un orden especial, concretamente ordenadas por el índice derecho de las consultas.

Como ya se dijo, usaremos una sola estructura de datos. Esta estructura de datos guardará información sobre el rango. Al principio este rango estará vacío. Cuando queremos responder la siguiente consulta (en el orden especial), simplemente extendemos o reducimos el rango, agregando/eliminando elementos en ambos extremos del rango actual, hasta transformarlo en el rango de la consulta. De este modo, solo necesitamos agregar o eliminar un solo elemento a la vez, que deberían ser operaciones bastante fáciles en nuestra estructura de datos.

Como cambiamos el orden de respuesta de las consultas, esto solo es posible cuando se permite responder las consultas en modo offline.

Implementación

En el algoritmo de Mo usamos dos funciones para agregar un índice y para eliminar un índice del rango que estamos manteniendo actualmente.


void remove(idx);  // TODO: eliminar el valor en idx de la estructura de datos
void add(idx);     // TODO: agregar el valor en idx a la estructura de datos
int get_answer();  // TODO: extraer la respuesta actual de la estructura de datos
 
int block_size;
 
struct Query {
    int l, r, idx;
    bool operator<(Query other) const
    {
        return make_pair(l / block_size, r) <
               make_pair(other.l / block_size, other.r);
    }
};
 
vector<int> mo_s_algorithm(vector<Query> queries) {
    vector<int> answers(queries.size());
    sort(queries.begin(), queries.end());
 
    // TODO: inicializar la estructura de datos
 
    int cur_l = 0;
    int cur_r = -1;
    // invariante: la estructura de datos siempre reflejará el rango [cur_l, cur_r]
    for (Query q : queries) {
        while (cur_l > q.l) {
            cur_l--;
            add(cur_l);
        }
        while (cur_r < q.r) {
            cur_r++;
            add(cur_r);
        }
        while (cur_l < q.l) {
            remove(cur_l);
            cur_l++;
        }
        while (cur_r > q.r) {
            remove(cur_r);
            cur_r--;
        }
        answers[q.idx] = get_answer();
    }
    return answers;
}

Según el problema podemos usar una estructura de datos distinta y modificar las funciones add/remove/get_answer en consecuencia. Por ejemplo, si se nos pide responder consultas de suma en un rango, usamos un entero simple como estructura de datos, que es $0$ al principio. La función add simplemente sumará el valor de la posición y luego actualizará la variable de respuesta. Por otro lado la función remove restará el valor de la posición y luego actualizará la variable de respuesta. Y get_answer solo devuelve el entero.

Para responder consultas de moda, podemos usar un árbol binario de búsqueda (p. ej. map<int, int>) para guardar cuántas veces aparece cada número en el rango actual, y un segundo árbol binario de búsqueda (p. ej. set<pair<int, int>>) para mantener los conteos de los números (p. ej. como pares conteo-número) en orden. El método add elimina el número actual del segundo BST, incrementa el conteo en el primero, e inserta el número de vuelta en el segundo. remove hace lo mismo, solo que decrementa el conteo. Y get_answer solo mira el segundo árbol y devuelve el mejor valor en $O(1)$ .

Complejidad

Ordenar todas las consultas tomará $O(Q \log Q)$ .

¿Y las demás operaciones? ¿Cuántas veces se llamará a add y remove?

Digamos que el tamaño del bloque es $S$ .

Si solo miramos todas las consultas cuyo índice izquierdo está en el mismo bloque, las consultas están ordenadas por el índice derecho. Por lo tanto llamaremos a add(cur_r) y remove(cur_r) solo $O(N)$ veces para todas estas consultas combinadas. Esto da $O(\frac{N}{S} N)$ llamadas para todos los bloques.

El valor de cur_l puede cambiar a lo sumo $O(S)$ entre dos consultas. Por lo tanto tenemos $O(S Q)$ llamadas adicionales de add(cur_l) y remove(cur_l).

Para $S \approx \sqrt{N}$ esto da $O((N + Q) \sqrt{N})$ operaciones en total. Así la complejidad es $O((N+Q)F\sqrt{N})$ donde $O(F)$ es la complejidad de las funciones add y remove.

Consejos para mejorar el tiempo de ejecución

Un tamaño de bloque de exactamente $\sqrt{N}$ no siempre ofrece el mejor tiempo de ejecución. Por ejemplo, si $\sqrt{N}=750$ entonces puede ocurrir que un tamaño de bloque de $700$ o $800$ se ejecute mejor. Más importante: no hay que calcular el tamaño del bloque en tiempo de ejecución; hay que dejarlo como const. La división por constantes está bien optimizada por los compiladores.
En los bloques impares ordenar el índice derecho en orden ascendente y en los bloques pares ordenarlo en orden descendente. Esto minimizará el movimiento del puntero derecho, ya que el ordenamiento normal moverá el puntero derecho desde el final de vuelta al principio al inicio de cada bloque. Con la versión mejorada este reinicio ya no es necesario.


bool cmp(pair<int, int> p, pair<int, int> q) {
    if (p.first / BLOCK_SIZE != q.first / BLOCK_SIZE)
        return p < q;
    return (p.first / BLOCK_SIZE & 1) ? (p.second < q.second) : (p.second > q.second);
}

Se puede leer sobre un enfoque de ordenamiento todavía más rápido aquí .

Descomposición por raíz cuadrada (Sqrt Decomposition)

Estructura de datos basada en descomposición por raíz cuadrada

Descripción

Implementación

Otros problemas

Algoritmo de Mo (Mo’s algorithm)

Implementación

Complejidad

Consejos para mejorar el tiempo de ejecución

Problemas de práctica