Algoritmos de búsqueda en el modelo de comparación

Capítulo 5 — Análisis de Algoritmos con Julia

Autor/a

Eric S. Téllez

1 Algoritmos de búsqueda en el modelo de comparación

En este capítulo analizaremos algoritmos de búsqueda en arreglos ordenados basados en funciones de comparación, con el objetivo de localizar elementos y posiciones específicas. Veremos algoritmos de peor caso, cuyo costo depende únicamente del tamaño \(n\) del arreglo, y algoritmos adaptables u oportunistas, cuyo costo depende de la posición \(p\) donde se encuentra la respuesta; estos últimos son muy veloces cuando la respuesta está cerca de donde se comienza a buscar, lo cual sucede con frecuencia en problemas reales.

Objetivos de aprendizaje

Al terminar este capítulo, el lector deberá ser capaz de:

Plantear el problema de búsqueda como la localización de una posición de inserción y derivar su cota inferior en el modelo de comparación.
Explicar la búsqueda binaria a través de su invariante, demostrar que es correcta y que es óptima en el peor caso.
Distinguir entre el costo en términos de \(n\) (tamaño del arreglo) y en términos de \(p\) (posición de la respuesta), y argumentar cuándo conviene cada tipo de algoritmo.
Implementar y analizar los algoritmos de búsqueda no acotada \(B_0\), \(B_1\) y \(B_2\) de Bentley y Yao (1976), y comprender la construcción general \(B_k\).
Medir experimentalmente el número de comparaciones y el tiempo de los algoritmos, y relacionar los resultados con el análisis.

¿Por qué dedicar un capítulo a buscar en un arreglo ordenado? En el Capítulo 4 vimos que ordenar cuesta \(\Theta(n \log n)\) comparaciones; ese costo se paga una sola vez y se amortiza si después hacemos muchas búsquedas. En Recuperación de Información, por ejemplo, un índice invertido guarda para cada término una lista de posteo ordenada con los identificadores de los documentos donde aparece; resolver una consulta implica buscar en estas listas miles o millones de veces, y la mayoría de las búsquedas comienzan muy cerca de donde terminó la anterior. Los algoritmos de este capítulo son los bloques de construcción de los algoritmos de intersección del Capítulo 6.

1.1 Problema

Sea \(A[1..n] = a_1, \cdots, a_n\) un arreglo ordenado con \(n \geq 1\) y un operador \(<\) (menor que); por simplicidad, también usaremos \(\leq\) (menor o igual que). Supondremos que no hay elementos duplicados en \(A\); nótese que esto no implica una pérdida de generalidad.

Si se permiten duplicados, la definición basada en \(x \leq A[p]\) localiza la primera ocurrencia de \(x\) (o el lugar donde debería estar). Para localizar la última ocurrencia basta con cambiar la condición a \(A[p] < x\) y ajustar el resultado; en Julia estas dos variantes son searchsortedfirst y searchsortedlast. Si además podemos preprocesar el arreglo, e.g., para conocer las zonas con duplicados, se pueden mejorar aún más los tiempos.

La tarea será: dado el valor \(x\) a ser localizado en \(A\), el problema consiste en determinar la posición de inserción \(p\), es decir, la posición donde \(x\) debería colocarse para que el arreglo siga ordenado. Formalmente, \(p\) es la única posición que cumple alguna de las siguientes condiciones:

si \(p = 1\) entonces \(x \leq A[p]\).
si \(2 \leq p \leq n\) entonces \(A[p-1] < x \leq A[p]\).
si \(p=n+1\) entonces \(A[n] < x\).

Dicho de otra forma, \(p\) es la posición del primer elemento que es mayor o igual que \(x\); si no existe tal elemento, \(p = n+1\). Note que esta formulación resuelve dos problemas a la vez: si \(x\) está en el arreglo, entonces \(A[p] = x\) y \(p\) es su posición; si no está, \(p\) nos dice dónde insertarlo, cuál es su sucesor \(A[p]\) y cuál su predecesor \(A[p-1]\). Por ejemplo, para \(A = [10, 20, 30, 40, 50, 60, 70]\) tenemos:

Tabla 1: Posiciones de inserción para distintos valores de \(x\) en un arreglo de \(n=7\) elementos.

\(x\)	\(p\)	Justificación
0	1	\(x \leq A[1] = 10\)
69	7	\(A[6] = 60 < 69 \leq A[7] = 70\)
70	7	\(A[6] = 60 < 70 \leq A[7] = 70\), y de hecho \(A[7] = x\)
71	8	\(A[7] = 70 < 71\), i.e., \(p = n+1\)

Usaremos los cuatro valores de Tabla 1 como prueba rápida de todos los algoritmos del capítulo, ya que cubren los casos frontera: antes del primer elemento, entre dos elementos, un elemento existente, y después del último elemento.

1.1.1 Costo del problema

Para \(A[1..n]\) y el valor \(x\) cuya posición de inserción se desea localizar, el resultado puede ser cualquiera de las \(n+1\) posiciones posibles, i.e., cada una de ellas es una instancia distinta del problema que el algoritmo debe distinguir de las demás.

En el modelo de comparación, la única manera de obtener información sobre \(x\) es compararlo con algún \(A[i]\), y cada comparación tiene dos resultados posibles: \(x \leq A[i]\) o \(x > A[i]\). Podemos entonces representar cualquier algoritmo como un árbol de decisión binario: cada nodo interno es una comparación, cada rama es uno de los dos resultados y cada hoja es una respuesta. Un árbol con \(n+1\) hojas tiene altura al menos \(\lceil \log_2 (n+1) \rceil\), por lo que cualquier algoritmo requiere al menos esa cantidad de comparaciones en su peor caso. La Figura 1 muestra el árbol de decisión óptimo para \(n = 7\): con \(3\) comparaciones se distinguen las \(8\) respuestas.

Figura 1: Árbol de decisión para localizar la posición de inserción en un arreglo de \(n=7\) elementos. Cada nodo interno compara \(x\) con un elemento del arreglo; las hojas son las \(n+1=8\) posiciones de inserción posibles. Cualquier algoritmo de comparación debe tener al menos \(n+1\) hojas, y por lo tanto altura \(\geq \lceil \log_2(n+1) \rceil = 3\).

Este argumento es el mismo que usamos para el problema de ordenamiento en el Capítulo 4, donde las hojas eran las \(n!\) permutaciones; la diferencia es que aquí solo hay \(n+1\) respuestas posibles, y por eso el problema es mucho más barato: \(\Theta(\log n)\) en lugar de \(\Theta(n \log n)\). Note que esta cota supone que \(n\) es conocido; más adelante veremos qué sucede cuando no lo es.

1.1.2 Búsqueda secuencial

Un algoritmo naïve compara \(x\) contra \(A[1], A[2], \ldots\) hasta encontrar el primer elemento mayor o igual que \(x\), o hasta agotar el arreglo.

Listado 1: Búsqueda secuencial

"""
    seqsearch(A, x, sp=1)

Búsqueda exhaustiva con inicio en `sp`; supone que todo elemento antes de `sp` es menor que `x`.
"""
function seqsearch(A, x, sp=1)
    n = length(A)
1    while sp <= n && x > A[sp]
        sp += 1
    end

2    sp
end

let S=[10, 20, 30, 40, 50, 60, 70]
    (seqsearch(S, 0), seqsearch(S, 69), seqsearch(S, 70), seqsearch(S, 71))
end

1: Avanzamos mientras no nos salgamos del arreglo y el elemento actual sea menor que \(x\).
2: Al salir del ciclo, o bien \(sp = n+1\), o bien \(A[sp]\) es el primer elemento con \(x \leq A[sp]\); en ambos casos \(sp\) es la posición de inserción.

(1, 7, 7, 8)

El parámetro sp (starting position) permite indicar una pista de dónde comenzar; lo usaremos en todos los algoritmos del capítulo, ya que modela el caso común donde se realizan búsquedas sucesivas y sabemos que la respuesta no está antes de cierta posición.

Si la respuesta es la posición \(p\), el algoritmo realiza exactamente \(p\) comparaciones de elementos (la \(p\)-ésima es la que falla), o \(n\) si \(p = n+1\); por lo tanto, su costo de peor caso es \(\Theta(n)\). Es un algoritmo adaptable en un sentido trivial: es muy rápido cuando \(p\) es muy pequeña.

En la práctica, la búsqueda secuencial es la más rápida para arreglos muy pequeños (unas decenas de elementos), ya que sus accesos a memoria son contiguos y predecibles, y el procesador puede ejecutar varias comparaciones por ciclo. Por esta razón, muchas implementaciones de búsqueda binaria cambian a búsqueda secuencial cuando el rango es pequeño. En listas no ordenadas, además, las heurísticas de auto-organización como mover al frente el elemento encontrado mejoran notablemente la búsqueda secuencial (Bentley y McGeoch 1985).

Sin embargo, dado que el arreglo está ordenado, se puede mejorar mucho el tiempo de búsqueda: una sola comparación \(x \leq A[i]\) descarta de golpe todas las posiciones a un lado de \(i\).

1.1.3 Búsqueda binaria

La idea de la búsqueda binaria es mantener un rango de posiciones candidatas donde sabemos que está la respuesta, y en cada paso comparar \(x\) con el elemento central del rango para descartar la mitad de los candidatos. El costo de búsqueda para cualquier instancia es \(O(\log n)\).

Listado 2: Búsqueda binaria

"""
    binarysearch(A, x, sp=1, ep=length(A))

Encuentra la posición de inserción de `x` en `A` suponiendo que la respuesta está en el rango `sp:ep+1`,
i.e., que todo elemento antes de `sp` es menor que `x` y que `x <= A[ep+1]` si existe.
"""
function binarysearch(A, x, sp=1, ep=length(A))
1    lo, hi = sp, ep + 1
2    while lo < hi
3        mid = (lo + hi) ÷ 2
4        if x <= A[mid]
            hi = mid
        else
            lo = mid + 1
        end
    end

5    lo
end

let S=[10, 20, 30, 40, 50, 60, 70]
    (binarysearch(S, 0), binarysearch(S, 69), binarysearch(S, 70), binarysearch(S, 71))
end

1: Invariante: la posición de inserción está en el rango lo:hi. Inicialmente son las \(n+1\) posiciones posibles (note que hi puede valer \(n+1\)).
2: Iteramos mientras haya más de un candidato.
3: Punto central del rango. Como lo < hi, se cumple lo <= mid < hi <= n+1, por lo que A[mid] siempre es un acceso válido.
4: Si \(x \leq A[mid]\) la respuesta está en lo:mid (la posición mid sigue siendo candidata); de lo contrario, \(A[mid] < x\) y la respuesta está en mid+1:hi. En ambos casos el invariante se preserva y el rango se reduce.
5: Al terminar, lo == hi y es el único candidato que queda.

(1, 7, 7, 8)

Para comprender la dinámica del algoritmo, es útil observar la traza de una ejecución. La siguiente función es idéntica a Listado 2, pero imprime el rango de candidatos en cada iteración.

function binarysearch_trace(A, x, sp=1, ep=length(A))
    lo, hi = sp, ep + 1
    while lo < hi
        mid = (lo + hi) ÷ 2
        if x <= A[mid]
            decision = "x <= A[$mid] = $(A[mid]) -> hi = $mid"
        else
            decision = "x > A[$mid] = $(A[mid]) -> lo = $(mid+1)"
        end
        println("candidatos $lo:$hi ($(hi - lo + 1) posiciones), $decision")
        if x <= A[mid]
            hi = mid
        else
            lo = mid + 1
        end
    end
    println("respuesta: $lo")
    lo
end

let S=[10, 20, 30, 40, 50, 60, 70]
    binarysearch_trace(S, 69)
end;

candidatos 1:8 (8 posiciones), x > A[4] = 40 -> lo = 5
candidatos 5:8 (4 posiciones), x > A[6] = 60 -> lo = 7
candidatos 7:8 (2 posiciones), x <= A[7] = 70 -> hi = 7
respuesta: 7

Compare la traza con la Figura 1: la búsqueda binaria es el árbol de decisión óptimo, recorrido de la raíz a una hoja.

1.1.3.1 Corrección y costo

La corrección de la búsqueda binaria se desprende del invariante del punto 1: es verdadero al inicio, cada iteración lo preserva (punto 4), y el ciclo termina porque el rango se reduce estrictamente en cada paso: mid < hi garantiza que hi = mid lo reduce, y lo = mid + 1 > lo también. Cuando termina, el invariante con un solo candidato nos da la respuesta.

Es muy fácil escribir mal una búsqueda binaria; los errores típicos son ciclos infinitos (cuando el rango no se reduce) y respuestas desplazadas en uno. Razonar explícitamente sobre el invariante es la mejor defensa. Un error famoso es el desbordamiento de (lo + hi) ÷ 2 cuando los índices son enteros de 32 bits y \(n\) es cercano a \(2^{31}\); la forma lo + (hi - lo) ÷ 2 lo evita. En Julia los índices son Int64, por lo que en la práctica no es un problema, pero conviene conocerlo.

Para el costo, sea \(r = hi - lo + 1\) el número de candidatos. La comparación con el punto central deja \(\lfloor r / 2 \rfloor\) candidatos de un lado y \(\lceil r / 2 \rceil\) del otro; en el peor caso quedan \(\lceil r/2 \rceil\). Partiendo de \(r = n+1\), después de \(k\) comparaciones quedan a lo más \(\lceil (n+1)/2^k \rceil\) candidatos, y el algoritmo termina cuando queda uno; esto sucede en cuanto \(2^k \geq n+1\). Por lo tanto, el número de comparaciones en el peor caso es exactamente

\[C_\text{bin}(n) = \lceil \log_2 (n+1) \rceil,\]

que coincide con la cota inferior que derivamos para el problema: la búsqueda binaria es óptima en el peor caso. Dicho de otra manera, ningún algoritmo de comparación puede garantizar menos comparaciones para todas las instancias. Es común escribir simplemente \(C_\text{bin}(n) = \log_2 n + O(1)\) o, de manera más gruesa, \(O(\log n)\). También es importante hacer notar que la búsqueda binaria es muy eficiente en memoria, ya que solo necesita dos índices adicionales.

Mejor caso, peor caso y caso promedio

Una característica interesante de Listado 2 es que todas las instancias cuestan casi lo mismo: \(\lfloor \log_2 (n+1) \rfloor\) o \(\lceil \log_2 (n+1) \rceil\) comparaciones, ya que la profundidad de todas las hojas del árbol de decisión difiere a lo más en uno. Esto significa que el peor caso, el mejor caso y el caso promedio coinciden en el término dominante. La variante que compara con tres resultados (\(<\), \(=\), \(>\)) y termina en cuanto encuentra \(x\) mejora el mejor caso a una sola comparación, pero no mejora el peor caso ni cambia significativamente el promedio, y hace dos comparaciones por iteración. La búsqueda binaria es un algoritmo de peor caso por excelencia: no se adapta a la instancia.

Julia ya cuenta con esta operación en su biblioteca estándar: searchsortedfirst(A, x) calcula exactamente la posición de inserción que definimos; la utilizaremos como oráculo para verificar nuestras implementaciones.

let S=[10, 20, 30, 40, 50, 60, 70]
    (searchsortedfirst(S, 0), searchsortedfirst(S, 69),
     searchsortedfirst(S, 70), searchsortedfirst(S, 71))
end

(1, 7, 7, 8)

1.1.3.2 Búsqueda \(m\)-aria y la jerarquía de memoria

Si fuera posible probar varios puntos a la vez, i.e., dividir el rango en \(m\) segmentos con una sola operación, el número de operaciones estaría acotado por \(\lceil \log_{m} (n+1) \rceil\). En el modelo de comparación esto no ayuda: cada operación requiere \(m-1\) comparaciones y el total \((m-1)\log_m n\) es mayor que \(\log_2 n\) para toda \(m > 2\) (ver ejercicios). Sin embargo, tiene mucho sentido para estructuras de datos que trabajan en diferentes niveles de memoria. En hardware moderno, la diferencia de velocidad entre el caché, la memoria principal y el disco es de varios órdenes de magnitud; el costo dominante no es la comparación, sino traer un bloque de un nivel lento al nivel rápido. Si un bloque contiene \(m-1\) claves que particionan el rango en \(m\) divisiones, el número de accesos a bloques se reduce a \(\log_m n\); esta es la idea detrás de los árboles B y de los índices de las bases de datos.

Otra manera de reducir las comparaciones es salirse del modelo: la búsqueda por interpolación usa los valores de \(A[lo]\) y \(A[hi]\) para estimar dónde debería estar \(x\), como cuando buscamos una palabra en un diccionario. Si los datos están distribuidos uniformemente, el costo esperado es \(O(\log \log n)\) (Perl et al. 1978; Knuth 1998, sec. 6.2.1), pero su peor caso es \(O(n)\) y en la práctica es sensible a la distribución de los datos.

1.1.3.3 Ejercicios

Modifique Listado 2 para que regrese también un valor booleano que indique si \(x\) está en el arreglo. ¿Cuántas comparaciones adicionales necesita?
Escriba binarysearchlast(A, x) que regrese la posición del último elemento menor o igual que \(x\) (o \(0\) si no existe); compare con searchsortedlast. Use ambas funciones para contar en \(O(\log n)\) el número de ocurrencias de \(x\) en un arreglo ordenado con duplicados.
Una variante común de búsqueda binaria mantiene el rango sp:ep sobre los elementos (no sobre las posiciones de inserción), itera mientras sp < ep y hace una comparación final x <= A[sp] para decidir entre sp y sp+1. Escríbala, verifique que es correcta y muestre que en el peor caso hace \(\lceil \log_2 n \rceil + 1\) comparaciones; ¿para qué valores de \(n\) hace una comparación más que Listado 2? ¿A qué se debe?
La búsqueda ternaria divide el rango en tres partes con dos comparaciones por iteración. Derive su número de comparaciones en el peor caso y muestre que es mayor que el de la búsqueda binaria. Generalice a \(m\) partes.
En el modelo de comparación, ¿cuál es el costo de peor caso de determinar si \(x\) está en un arreglo no ordenado? Argumente con un adversario que responde a las comparaciones de la peor manera posible.
Implemente la búsqueda por interpolación. Compare su número de comparaciones contra la búsqueda binaria en arreglos de enteros aleatorios uniformes y en arreglos con distribución muy sesgada, e.g., \(A[i] = 2^i\) o \(A[i] = i^3\). Construya una instancia donde haga \(\Omega(n)\) comparaciones.
Mida el tiempo de seqsearch y binarysearch para arreglos pequeños (\(n \in \{4, 8, 16, 32, 64, 128\}\)), promediando sobre muchas búsquedas aleatorias. ¿A partir de qué \(n\) la búsqueda binaria es más rápida? Modifique binarysearch para que cambie a búsqueda secuencial cuando el rango sea menor a ese umbral y mida la mejora.

1.2 Contando comparaciones

Para analizar experimentalmente los algoritmos de este capítulo nos interesa contar el número de comparaciones de cada búsqueda, y no solo medir tiempos. En nuestras implementaciones cada acceso A[i] va seguido de exactamente una comparación con \(x\), así que basta con contar los accesos. Una manera de hacerlo sin modificar los algoritmos es envolver el arreglo en un tipo que implemente la interfaz de AbstractVector y cuente cada lectura.

Listado 3: Arreglo que cuenta accesos; cada acceso corresponde a una comparación en nuestros algoritmos

1mutable struct CountingArray{T} <: AbstractVector{T}
    data::Vector{T}
    count::Int
end

CountingArray(data) = CountingArray(data, 0)
2Base.size(C::CountingArray) = size(C.data)
Base.IndexStyle(::Type{<:CountingArray}) = IndexLinear()

3function Base.getindex(C::CountingArray, i::Int)
    C.count += 1
    C.data[i]
end

"""
    comparaciones(f, A, x, args...)

Ejecuta la búsqueda `f(A, x, args...)` y regresa la tupla (posición, número de comparaciones)
"""
function comparaciones(f, A::CountingArray, x, args...)
    A.count = 0
    p = f(A, x, args...)
    (p, A.count)
end

let S = CountingArray([10, 20, 30, 40, 50, 60, 70])
    (comparaciones(seqsearch, S, 69), comparaciones(binarysearch, S, 69))
end

1: El tipo se declara como subtipo de AbstractVector, por lo que funciona en cualquier lugar donde se espera un arreglo; es mutable para poder actualizar el contador.
2: La interfaz mínima de AbstractArray en Julia requiere size y getindex; con esto obtenemos gratis length, iteración, impresión, etc.
3: Cada lectura incrementa el contador y delega al arreglo interno.

((7, 7), (7, 3))

Este es un ejemplo del uso de la programación genérica de Julia que vimos en los Capítulos 1 y 3: los algoritmos no cambiaron, pero ahora podemos instrumentarlos. Los costos observados, \(7\) comparaciones para seqsearch y \(3 = \lceil \log_2 8 \rceil\) para binarysearch, coinciden con el análisis. Verifiquemos ahora el peor caso de la búsqueda binaria para varios tamaños, probando las \(n+1\) instancias posibles de cada tamaño:

for n in (1, 7, 8, 100, 1000, 10^6)
    A = CountingArray(collect(1:n))
    peor = maximum(comparaciones(binarysearch, A, x)[2] for x in 1:n+1)
    println("n = $n:\tpeor caso = $peor\tceil(log2(n+1)) = $(ceil(Int, log2(n+1)))")
end

n = 1:  peor caso = 1   ceil(log2(n+1)) = 1
n = 7:  peor caso = 3   ceil(log2(n+1)) = 3
n = 8:  peor caso = 4   ceil(log2(n+1)) = 4
n = 100:    peor caso = 7   ceil(log2(n+1)) = 7
n = 1000:   peor caso = 10  ceil(log2(n+1)) = 10
n = 1000000:    peor caso = 20  ceil(log2(n+1)) = 20

1.3 Búsqueda no acotada

Cuando el tamaño del arreglo es demasiado grande, o la relación \(p / n\) es significativamente pequeña, la búsqueda acotada no es la mejor opción. Aun cuando en la práctica el límite superior \(n\) podría estar determinado, y por lo tanto se pueden resolver búsquedas en \(O(\log n)\), es posible obtener una cota relativa a \(p\), independiente de \(n\), por lo que los casos de interés se verán beneficiados. Dos situaciones motivan este tipo de búsqueda:

El límite \(n\) es desconocido o no existe. Por ejemplo, buscar en un flujo de datos ordenado, o encontrar el menor entero \(i\) tal que \(f(i) \geq x\) para una función \(f\) monótona que podemos evaluar pero cuyo dominio no está acotado. En este escenario la búsqueda binaria ni siquiera puede comenzar.
El límite \(n\) es conocido, pero la respuesta suele estar cerca del inicio. Cuando se realizan búsquedas sucesivas de valores crecientes \(x_1 < x_2 < \cdots\) sobre el mismo arreglo, cada búsqueda puede comenzar en la posición donde terminó la anterior (la pista sp), y la respuesta suele estar a pocas posiciones de distancia. Esto es exactamente lo que ocurre al intersectar listas de posteo (ver Capítulo 6). Aquí \(p\) debe entenderse como la posición relativa a la pista.

A este tipo de búsqueda también se le conoce como búsqueda exponencial, galloping search o finger search: la pista sp es un “dedo” que señala dónde comenzar, y el costo depende de la distancia entre el dedo y la respuesta, no del tamaño del arreglo.

En ambos casos queremos un algoritmo cuyo costo esté acotado por una función de \(p\). Una estrategia simple y poderosa es la siguiente:

Determinar un buen rango que contenga la respuesta.
Aplicar búsqueda binaria en ese rango para obtener la respuesta.

Bentley y Yao (1976) describen en detalle una familia de algoritmos casi óptimos para la búsqueda no acotada siguiendo la estrategia anteriormente mencionada; en particular, poniendo un énfasis importante en la determinación del rango. Lo consiguen mediante la definición de los algoritmos \(B_0, B_1, B_2, \ldots\) que se describen en el resto de la sección. Para cada algoritmo \(B_k\) denotaremos con \(F_k\) la secuencia de posiciones que se prueban para determinar el rango, y con \(C_k(p)\) su costo en comparaciones cuando la respuesta es \(p\); para simplificar la notación supondremos que la pista es \(sp = 1\), de modo que \(p\) es tanto la posición absoluta como la relativa.

1.3.1 Algoritmo \(B_0\) (búsqueda unaria)

Es el algoritmo más simple, y ya lo vimos con anterioridad (Listado 1): realiza una búsqueda exhaustiva de la posición de inserción, haciendo pruebas para toda posición \(x \leq A[1], x \leq A[2], \cdots, x \leq A[p]\), por lo que su costo es \(p\). Entonces:

\(F_0 = 1, 2, 3, \cdots\);
\(C_0(p) = p\); no requiere búsqueda binaria.

Note que \(B_0\) es un algoritmo de búsqueda no acotada perfectamente válido: no necesita conocer \(n\). Su problema es que el costo crece linealmente con \(p\).

1.3.2 Algoritmo \(B_1\) (búsqueda doblada: doubling search/galloping)

La idea de \(B_1\) es probar posiciones cada vez más alejadas, duplicando la distancia en cada paso: \(1, 2, 4, 8, \ldots\). Esto es, se compara \(x\) con \(A[2^0], A[2^1], A[2^2], \ldots\) hasta encontrar el primer exponente \(j\) tal que \(x \leq A[2^{j}]\) (o hasta salirse del arreglo). En ese momento sabemos que \(A[2^{j-1}] < x \leq A[2^{j}]\), es decir, que la respuesta está en el rango \(2^{j-1}+1 : 2^{j}\), y aplicamos búsqueda binaria en él.

\(F_1 = 2^0, 2^1, 2^2, \cdots\), i.e., las potencias de dos;
\(C_1(p) = \underbrace{\lceil \log_2 p \rceil + 1}_\text{rango} + \underbrace{C_\text{bin}(2^{\lceil \log_2 p \rceil - 1})}_\text{búsqueda binaria} = 2 \lceil \log_2 p \rceil + O(1) = 2\log_2 p + O(1).\)

La explicación viene a continuación. El primer exponente \(j\) con \(2^j \geq p\) es \(j = \lceil \log_2 p \rceil\), y para encontrarlo se realizan \(j+1\) comparaciones (los exponentes \(0, 1, \ldots, j\)). El rango resultante tiene \(2^{j} - 2^{j-1} = 2^{j-1}\) posiciones, por lo que la búsqueda binaria cuesta \(\lceil \log_2 (2^{j-1} + 1) \rceil = j\) comparaciones. En total, \(2j + 1 = 2\lceil \log_2 p \rceil + 1\) comparaciones: aproximadamente el doble de lo que costaría una búsqueda binaria si supiéramos que \(n = p\), pero sin necesidad de conocer \(n\).

Hay una bonita relación entre \(B_1\) y la codificación de enteros. El código \(\gamma\) de Elias representa un entero \(p\) con \(2\lfloor \log_2 p \rfloor + 1\) bits: una parte unaria que indica \(\lfloor \log_2 p \rfloor\) y una parte binaria con los bits restantes. La secuencia de respuestas sí/no que produce \(B_1\) es esencialmente este código, y la de \(B_2\) corresponde al código \(\delta\) de Elias. Bentley y Yao (1976) notan esta conexión, que es muy relevante en Recuperación de Información, donde estos códigos se usan para comprimir las listas de posteo.

El siguiente código implementa \(B_1\); note que las posiciones probadas son relativas a la pista sp.

Listado 4: Algoritmo \(B_1\), búsqueda doblada o galloping

function doublingsearch(A, x, sp=1)
    n = length(A)
    p = 0
    i = 1

1    while sp+i <= n && A[sp+i] < x
        p = i
2        i += i
    end

3    binarysearch(A, x, sp + p, min(n, sp+i))
end

let S=[10, 20, 30, 40, 50, 60, 70]
    (doublingsearch(S, 0), doublingsearch(S, 69), doublingsearch(S, 70), doublingsearch(S, 71))
end

1: Determinación del rango: avanzamos mientras el elemento probado sea menor que \(x\) y no nos salgamos del arreglo. Al salir, p es el último desplazamiento probado con A[sp+p] < x (o cero si no hubo ninguno) e i es el primero con x <= A[sp+i] (o se salió del arreglo).
2: El desplazamiento se duplica en cada iteración: \(1, 2, 4, 8, \ldots\)
3: Aplicar un algoritmo de búsqueda eficiente en el rango que contiene la respuesta; el min evita salirnos del arreglo cuando la respuesta está cerca del final.

(1, 7, 7, 8)

Veamos el número de comparaciones para un arreglo de \(n = 100\) elementos con respuestas en distintas posiciones; compare con las \(\lceil \log_2 101 \rceil = 7\) comparaciones que siempre realiza la búsqueda binaria.

let A = CountingArray(collect(1:100))
    for x in (1, 5, 10, 37, 90)
        b0 = comparaciones(seqsearch, A, x)[2]
        bin = comparaciones(binarysearch, A, x)[2]
        b1 = comparaciones(doublingsearch, A, x)[2]
        println("p = $x:\tB0 = $b0\tbinaria = $bin\tB1 = $b1")
    end
end

p = 1:  B0 = 1  binaria = 7 B1 = 3
p = 5:  B0 = 5  binaria = 7 B1 = 5
p = 10: B0 = 10 binaria = 6 B1 = 9
p = 37: B0 = 37 binaria = 7 B1 = 12
p = 90: B0 = 90 binaria = 7 B1 = 12

Es importante saber cuándo usar un algoritmo u otro, y por lo tanto determinar cuándo \(2\log_2{p} + O(1) < \log_2 n + O(1).\) Para simplificar este análisis ignoraremos las constantes: \[\begin{align} 2\log_2{p} & < \log_2 n, \\ \log_2{p^2} & < \log_2 n, \\ p^2 & < n, \\ p & < \sqrt{n}; \\ \end{align}\] esto indica que si \(p\) es menor a \(\sqrt{n}\) entonces hay una ventaja al usar \(B_1\); lo cual nos dice que para posiciones cercanas al inicio (o a la pista) el uso de \(B_1\) puede llevar a búsquedas más veloces. En el ejemplo anterior, con \(n = 100\), el cruce ocurre un poco antes de \(\sqrt{100} = 10\): en \(p = 5\) \(B_1\) todavía gana y en \(p = 10\) ya perdió; el desplazamiento respecto a \(\sqrt{n}\) se debe a las constantes que ignoramos. Note que en la práctica es necesario tener en cuenta la memoria; interesantemente, para valores de \(p\) pequeños es posible que esto beneficie aún más al algoritmo, ya que los accesos están concentrados en una zona pequeña del arreglo que probablemente se mantenga en caché.

Es cierto que estos algoritmos son oportunistas, pero hay aplicaciones donde esto realmente sucede. En el peor caso, el costo será apenas dos veces el óptimo.

1.3.3 Algoritmo \(B_2\) (búsqueda doblemente doblada, doubling-doubling search)

Para entender \(B_2\) conviene reinterpretar \(B_1\): la fase de determinación del rango es una búsqueda secuencial (\(B_0\)) sobre los exponentes \(j = 0, 1, 2, \ldots\) para encontrar \(m = \lceil \log_2 p \rceil\), el exponente de la potencia de dos que acota a \(p\); una vez conocido \(m\), la búsqueda binaria dentro del rango cuesta \(m\) comparaciones. Como \(B_0\) es lineal, encontrar \(m\) cuesta \(m + 1\) comparaciones y de ahí el factor \(2\) en \(C_1(p)\).

La pregunta natural es: ¿por qué buscar el exponente \(m\) con \(B_0\) y no con \(B_1\)? Eso es exactamente \(B_2\): para determinar \(m\), se prueban los exponentes \(1, 2, 4, 8, \ldots\), es decir, las posiciones \(2^{2^0}, 2^{2^1}, 2^{2^2}, \ldots = 2, 4, 16, 256, 65536, \ldots\), hasta encontrar la primera con \(x \leq A[2^{2^i}]\); con esto sabemos que \(m\) está en el rango de exponentes \(2^{i-1}+1 : 2^i\) y lo localizamos con una búsqueda binaria sobre los exponentes, comparando \(x\) con \(A[2^e]\) para el exponente central \(e\) del rango. Finalmente, con \(m\) determinado, hacemos la búsqueda binaria usual en el rango de posiciones \(2^{m-1}+1 : 2^m\).

\(F_2 = 2^{2^0}, 2^{2^1}, 2^{2^2}, \cdots\) para la primera fase, seguida de una búsqueda binaria sobre exponentes y una sobre posiciones;
\(C_2(p) = \log_2 p + 2\log_2 \log_2 p + O(1)\), que resulta de sumar los costos de las tres fases:

\[C_2(p) = \underbrace{\lceil \log_2 \log_2 p \rceil + 1}_\text{exponente del exponente} + \underbrace{\lceil \log_2 \log_2 p \rceil}_\text{búsqueda binaria en exponentes} + \underbrace{\lceil \log_2 p \rceil}_\text{búsqueda binaria en posiciones} + O(1).\]

Note cómo el término de mayor peso, \(\log_2 p\), es ahora el mismo que el de una búsqueda binaria que conociera \(n = p\); el precio de no conocer \(n\) se redujo de \(\log_2 p\) a \(2 \log_2 \log_2 p\), que en términos prácticos crece tan lento que se puede ver como una constante pequeña.

Listado 5: Algoritmo \(B_2\), búsqueda doblemente doblada

function doublingdoublingsearch(A, x, sp=1)
    n = length(A)
1    pos(e) = e >= 62 ? typemax(Int) : sp + (1 << e)
2    reached(e) = pos(e) > n || x <= A[pos(e)]

3    lo, hi = -1, 1
    while !reached(hi)
        lo = hi
        hi += hi
    end

4    lo += 1
    while lo < hi
        mid = (lo + hi) ÷ 2
        if reached(mid)
            hi = mid
        else
            lo = mid + 1
        end
    end

5    first = hi == 0 ? sp : pos(hi - 1)
    binarysearch(A, x, first, min(n, pos(hi)))
end

let S=[10, 20, 30, 40, 50, 60, 70]
    (doublingdoublingsearch(S, 0), doublingdoublingsearch(S, 69),
     doublingdoublingsearch(S, 70), doublingdoublingsearch(S, 71))
end

1: pos(e) es la posición absoluta que corresponde al exponente \(e\), i.e., \(sp + 2^e\). Para exponentes enormes satura en lugar de desbordar el entero; recuerde que \(2^{63}\) no cabe en un Int64.
2: reached(e) responde si la respuesta está en o antes de la posición pos(e); una posición fuera del arreglo cuenta como alcanzada.
3: Fase 1 (\(B_0\) sobre los exponentes de los exponentes): se prueban los exponentes hi \(= 1, 2, 4, 8, \ldots\), i.e., las posiciones \(sp + 2, sp + 4, sp + 16, sp + 256, \ldots\), hasta alcanzar la respuesta. Al salir, el exponente buscado \(m\) está en lo+1:hi; el valor inicial \(-1\) de lo indica que no hay información sobre el exponente \(0\).
4: Fase 2: búsqueda binaria sobre el rango de exponentes lo:hi para encontrar el menor exponente \(m\) que alcanza la respuesta. Esta fase es idéntica a Listado 2, pero en lugar de preguntar por A[mid] pregunta por A[pos(mid)].
5: Fase 3: ya sabemos que \(A[pos(m-1)] < x \leq A[pos(m)]\) (o que \(pos(m)\) está fuera del arreglo); búsqueda binaria en ese rango. El caso hi == 0 ocurre cuando la respuesta está en \(sp\) o \(sp+1\).

(1, 7, 7, 8)

Comparemos ahora los tres algoritmos para un arreglo de un millón de elementos (la búsqueda binaria realiza \(20\) o \(21\) comparaciones en todas las instancias), con respuestas a distintas distancias de la pista:

let A = CountingArray(collect(1:2^20))
    println("p\tbinaria\tB1\tB2\t2*ceil(log2 p)+1")
    for x in (2, 16, 256, 1024, 4096, 65536, 2^20)
        bin = comparaciones(binarysearch, A, x)[2]
        b1 = comparaciones(doublingsearch, A, x)[2]
        b2 = comparaciones(doublingdoublingsearch, A, x)[2]
        println("$x\t$bin\t$b1\t$b2\t$(2ceil(Int, log2(x)) + 1)")
    end
end

p   binaria B1  B2  2*ceil(log2 p)+1
2   21  3   4   3
16  20  8   7   9
256 20  16  13  17
1024    20  20  17  21
4096    20  24  19  25
65536   20  32  23  33
1048576 20  39  26  41

Observe que (i) para \(p\) pequeña los algoritmos no acotados ganan por un amplio margen, (ii) \(B_1\) empata con la búsqueda binaria justo en \(p = \sqrt{n} = 1024\), como predijo el análisis, y (iii) \(B_2\) sigue siendo competitivo mucho más allá de \(\sqrt{n}\), y en el peor caso (\(p = n\)) paga un sobrecosto moderado respecto a la búsqueda binaria. La última columna muestra la cota \(C_1(p) = 2\lceil \log_2 p \rceil + 1\) que derivamos para \(B_1\); para \(B_2\) puede verificar que los valores son consistentes con \(\log_2 p + 2\log_2 \log_2 p + O(1)\). Note que las pruebas que caen fuera del arreglo no cuestan comparaciones, ya que se resuelven comparando índices; por eso, cuando la respuesta está cerca del final, los algoritmos no acotados cuestan un poco menos de lo que indican las fórmulas.

Simplificando las expresiones, la relación que nos describe cuándo es mejor usar \(B_2\) que la búsqueda binaria es como sigue:

\[\begin{align} \log_2{p} + 2\log_2{\log_2{p}} &< \log_2 n\\ \log_2{p} + \log_2{\log^2_2{p}} &< \log_2 n\\ \log_2{(p \log^2_2{p})} &< \log_2 {n}\\ p \log^2_2{p} & < n.\\ \end{align}\]

Es decir, \(B_2\) es preferible mientras \(p\) sea menor que, aproximadamente, \(n / \log_2^2 n\): un rango mucho más amplio que el \(p < \sqrt{n}\) de \(B_1\). Por ejemplo, si \(p = \sqrt{n}\) entonces \(p \log^2_2 p = \frac{1}{4}\sqrt{n} \log^2_2 n\), que es claramente menor que \(n\) para valores relativamente grandes de \(n\) (a partir de \(n \approx 2^{16}\)), por lo que \(B_2\) funciona mejor que \(B_1\) para \(p\) relativamente grandes. Para \(p\) muy pequeña, en cambio, la diferencia entre ambos es de una o dos comparaciones en uno u otro sentido, ya que \(B_2\) paga el sobrecosto de sus fases adicionales sin obtener mucho beneficio; en los experimentos anteriores la ventaja de \(B_2\) se vuelve sistemática a partir de \(p\) de unas pocas decenas.

1.3.4 Algoritmo \(B_k\) y el límite de la estrategia

La construcción anterior se puede repetir: \(B_k\) determina el exponente \(m = \lceil \log_2 p \rceil\) usando \(B_{k-1}\) sobre la secuencia de potencias de dos, y después realiza una búsqueda binaria de \(m\) comparaciones en el rango \(2^{m-1}+1 : 2^{m}\). Esto da la recurrencia

\[C_k(p) = C_{k-1}(\lceil \log_2 p \rceil) + \lceil \log_2 p \rceil + O(1), \qquad C_0(p) = p,\]

cuya solución es

\[C_k(p) = \log_2 p + \log_2^{(2)} p + \cdots + \log_2^{(k-1)} p + 2 \log_2^{(k)} p + O(k),\]

donde \(\log_2^{(j)}\) denota el logaritmo iterado \(j\) veces, i.e., \(\log_2^{(1)} p = \log_2 p\) y \(\log_2^{(j)} p = \log_2 \log_2^{(j-1)} p\). Puede verificar que para \(k=1\) y \(k=2\) se recuperan \(C_1\) y \(C_2\). Cada nivel adicional reemplaza el factor \(2\) del último término por un término adicional mucho más pequeño.

El logaritmo iterado \(\log^\star_2 n\) se define como el número de veces que hay que aplicar \(\log_2\) a \(n\) para obtener un valor menor o igual que \(1\), i.e., la \(k\) más pequeña tal que \(\log_2^{(k)} n \leq 1\). Crece extraordinariamente lento: \(\log^\star_2 2^{65536} = 5\), y \(2^{65536}\) es un número con casi veinte mil dígitos decimales. Para cualquier \(n\) de interés práctico, \(\log^\star_2 n \leq 5\).

Bentley y Yao (1976) muestran que al tomar \(k = \log^\star_2 p\) se obtiene un algoritmo con costo \(\log_2 p + \log_2^{(2)} p + \log_2^{(3)} p + \cdots + O(\log^\star_2 p)\), y demuestran una cota inferior para la búsqueda no acotada que muestra que ningún algoritmo puede mejorar esencialmente esta expresión; de ahí el título de su artículo, an almost optimal algorithm. En términos prácticos, \(B_1\) y \(B_2\) son los algoritmos relevantes: la constante escondida en \(O(k)\) crece con \(k\), mientras que la ganancia \(\log^{(k)}_2 p\) se vuelve despreciable muy rápido.

1.3.5 Comparación experimental

La Figura 2 resume el comportamiento de los algoritmos en un arreglo de \(n = 2^{20}\) elementos, midiendo el número de comparaciones en función de la posición \(p\) de la respuesta (en escala logarítmica). La búsqueda binaria es una línea horizontal: no se adapta a la instancia. \(B_1\) y \(B_2\) crecen con \(\log_2 p\) con pendientes \(2\) y \(1\) respectivamente; \(B_1\) cruza la búsqueda binaria en \(\sqrt{n}\) y \(B_2\) se mantiene cerca de ella hasta el final.

using Plots, LaTeXStrings

n = 2^20
A = CountingArray(collect(1:n))
P = unique(round.(Int, 2 .^ (0:0.25:20)))
cbin = [comparaciones(binarysearch, A, p)[2] for p in P]
cb1 = [comparaciones(doublingsearch, A, p)[2] for p in P]
cb2 = [comparaciones(doublingdoublingsearch, A, p)[2] for p in P]

plot(P, cbin, xscale=:log2, label="búsqueda binaria", lw=2,
     xlabel=L"p", ylabel="comparaciones", legend=:topleft)
plot!(P, cb1, label=L"B_1", lw=2)
plot!(P, cb2, label=L"B_2", lw=2)
vline!([sqrt(n)], label=L"\sqrt{n}", linestyle=:dash, color=:gray)

Figura 2: Número de comparaciones en función de la posición de la respuesta \(p\), para un arreglo de \(n = 2^{20}\) elementos. La línea punteada vertical marca \(\sqrt{n}\), donde \(B_1\) y la búsqueda binaria empatan.

Recuerde que el número de comparaciones es una medida del modelo, no del tiempo real. Los algoritmos no acotados acceden a posiciones cercanas a la pista, lo que favorece al caché, pero realizan más operaciones aritméticas por comparación; la Actividad de este capítulo le pedirá medir ambas cosas.

1.3.5.1 Ejercicios

Derive paso a paso el costo \(C_2(p)\) a partir de la descripción en tres fases de Listado 5, cuidando los términos de redondeo. Compare su expresión con los números obtenidos experimentalmente.
¿Cuándo \(B_1\) es mejor que \(B_2\)? Plantee la desigualdad entre \(C_1(p)\) y \(C_2(p)\), incluyendo constantes razonables para los términos \(O(1)\), y verifique experimentalmente con CountingArray el valor de \(p\) donde se cruzan.
Haga un pseudocódigo para \(B_k\) en general. Sugerencia: escriba una función recursiva unbounded(k, reached) que reciba el predicado reached(i) que indica si la respuesta está en o antes de la posición \(i\); \(B_0\) es la búsqueda secuencial sobre \(i\), y \(B_k\) llama a \(B_{k-1}\) con el predicado e -> reached(2^e) para encontrar el exponente.
Resuelva la recurrencia de \(C_k\) y verifique la fórmula cerrada para \(k = 3\). ¿Para qué valores de \(p\) (si existen) \(B_3\) hace menos comparaciones que \(B_2\) en su implementación?
En Listado 4 la búsqueda binaria recibe el rango sp+p : sp+i, aunque ya sabemos que A[sp+p] < x. Ajuste el rango para no repetir información, mida la diferencia en comparaciones y explique por qué no cambia el análisis asintótico.
Suponga que debe buscar \(m\) valores ordenados \(x_1 < x_2 < \cdots < x_m\) en el mismo arreglo \(A\) de tamaño \(n\). Compare el costo total de (a) \(m\) búsquedas binarias independientes y (b) \(m\) búsquedas con \(B_1\) donde cada una comienza donde terminó la anterior. Muestre que en (b) el costo total es \(O(m \log (n/m))\) cuando las respuestas están uniformemente espaciadas; ¿qué sucede cuando todas las respuestas están concentradas en una zona pequeña del arreglo? Este problema reaparece en el Capítulo 6.
Use doublingsearch para calcular \(\lfloor \sqrt{x} \rfloor\) para un entero \(x \geq 0\) sin conocer una cota superior, definiendo un arreglo virtual \(A[i] = i^2\) como el VirtualRange que se usó en las pruebas. ¿Cuántas comparaciones realiza en función de \(x\)?

1.4 Más allá de los arreglos

Los arreglos ordenados permiten búsquedas muy eficientes, pero insertar un nuevo elemento cuesta \(O(n)\), ya que hay que desplazar los elementos posteriores. Cuando el conjunto cambia con frecuencia, se requieren estructuras de datos dinámicas que mantengan el orden: los árboles binarios de búsqueda balanceados garantizan búsqueda, inserción y eliminación en \(O(\log n)\) en el peor caso, y la skip list que vimos en el Capítulo 4 obtiene las mismas cotas con alta probabilidad usando aleatorización (Pugh 1990).

Las ideas de este capítulo se trasladan directamente a estas estructuras. Una búsqueda en un árbol balanceado es un árbol de decisión como el de Figura 1, con la diferencia de que su forma la determina la estructura y no el algoritmo. La búsqueda no acotada corresponde a la búsqueda con dedo (finger search): si guardamos una referencia (el dedo) al último nodo visitado, podemos localizar un elemento a distancia \(d\) del dedo en \(O(\log d)\) en lugar de \(O(\log n)\); las estructuras que lo soportan se conocen como finger trees, y la skip list lo permite con una pequeña modificación.

1.4.0.1 Ejercicios

¿Qué es un árbol binario de búsqueda? ¿Cuál es el costo de búsqueda en un árbol? ¿Qué se debe hacer para asegurar los costos, i.e., por qué un árbol binario de búsqueda sin más garantías puede degradarse a una lista?
¿Qué es un finger tree? Describa cómo cambia el algoritmo de búsqueda cuando se parte del dedo en lugar de la raíz.
¿Cuál es el costo de búsqueda de la skip list? Relacione el argumento con los niveles de la estructura y la probabilidad \(p\) con que se construyen.
¿Cómo se puede hacer la skip list adaptativa, i.e., con búsqueda con dedo? ¿Qué otra forma de adaptarse a las consultas podría aplicar, por ejemplo, si algunos elementos se consultan con mucha mayor frecuencia que otros?
Investigue, implemente y pruebe un árbol binario de búsqueda y compare su número de comparaciones contra la búsqueda binaria en un arreglo con los mismos elementos, insertándolos en orden aleatorio y en orden creciente.

1.5 Lecturas

El tratamiento clásico de la búsqueda en arreglos ordenados es el capítulo 6 de Knuth (1998), Searching; en particular, la sección 6.1 cubre la búsqueda secuencial y la sección 6.2.1 la búsqueda binaria, sus variantes y la búsqueda por interpolación, con un análisis cuidadoso del número de comparaciones. En Cormen et al. (2022) la búsqueda binaria aparece en los ejercicios del capítulo 2 y los árboles binarios de búsqueda en el capítulo 12.

El artículo de Bentley y Yao (1976) es corto y muy legible; además de definir la familia \(B_k\), demuestra la cota inferior para la búsqueda no acotada y discute la relación con la representación de enteros. La lectura de Bentley y McGeoch (1985) complementa el panorama con heurísticas de auto-organización para la búsqueda secuencial, mostrando cómo adaptarse a la distribución de las consultas en lugar de a la posición de la respuesta. La búsqueda por interpolación y su costo \(O(\log \log n)\) esperado se analizan en Perl et al. (1978). Finalmente, Pugh (1990) presenta la skip list, cuya búsqueda es una generalización probabilística de las ideas de este capítulo.

1.6 Material audio-visual

En el siguiente video se adentrará en diferentes estrategias de búsqueda, notoriamente aquellas que llamaremos oportunistas o adaptables (adaptive). Estas técnicas nos permitirán tomar provecho de instancias sencillas de problemas e incrementar el desempeño en ese tipo de instancias.

Tenga en cuenta que, honrando la literatura, usaremos de forma indiscriminada listas ordenadas como sinónimo de arreglos ordenados.

1.7 Actividades

1.7.1 Actividad 1: Búsquedas sucesivas con dedo

En aplicaciones como la intersección de listas de posteo, las búsquedas no llegan aisladas sino en secuencias ordenadas \(x_1 < x_2 < \cdots < x_m\) sobre el mismo arreglo \(A[1..n]\). Esto abre dos estrategias: resolver cada consulta con una búsqueda binaria independiente, o resolverlas en orden usando como pista sp la posición donde terminó la consulta anterior (búsqueda con dedo), con lo que el costo de cada búsqueda depende de la distancia \(p\) a la respuesta anterior y no de \(n\). El objetivo de esta actividad es determinar experimentalmente cuándo conviene cada estrategia, y relacionarlo con los análisis del capítulo.

Todos los datos son sintéticos y se generan con el código de Listado 6, con semillas fijas, para que todos los reportes trabajen sobre las mismas instancias; no utilice archivos externos ni cambie las semillas ni los parámetros predeterminados. Las consultas siguen tres distribuciones, elegidas para que la distancia \(p\) entre respuestas consecutivas se comporte de manera distinta:

Uniforme: \(m\) valores distintos tomados uniformemente del universo \(1:U\); las distancias \(p\) se concentran alrededor de \(n/m\).
Ráfagas: \(k = 8\) centros aleatorios y \(m/k\) consultas en una ventana pequeña alrededor de cada uno; las distancias son casi siempre minúsculas, con unos pocos saltos enormes.
Cola pesada: \(m\) huecos \(h_i = u_i^{-1/\alpha}\) con \(u_i\) uniforme en \((0,1]\) y \(\alpha = 1.5\), reescalados al universo; la mayoría de los huecos son pequeños pero existen huecos muy grandes, lo que imita a las listas de posteo reales.

Listado 6: Generadores de datos sintéticos para la actividad; todas las funciones usan semillas fijas

using Random

const SEMILLA = 2026

"""
    genera_arreglo(n; g=8, semilla=SEMILLA)

Arreglo ordenado de `n` enteros distintos, acumulando huecos uniformes en `1:g`.
"""
function genera_arreglo(n; g=8, semilla=SEMILLA)
    rng = Xoshiro(semilla)
    cumsum(rand(rng, 1:g, n))
end

"""
    consultas_uniformes(U, m; semilla)

`m` valores distintos tomados uniformemente del universo `1:U`, ordenados.
"""
function consultas_uniformes(U, m; semilla=SEMILLA + 1)
    rng = Xoshiro(semilla)
    sort!(randperm(rng, U)[1:m])
end

"""
    consultas_rafagas(U, m; k=8, semilla)

`m` consultas repartidas en `k` ráfagas; cada ráfaga toma `m ÷ k` valores distintos
de una ventana de tamaño `4(m ÷ k)` que comienza en un punto aleatorio del universo.
"""
function consultas_rafagas(U, m; k=8, semilla=SEMILLA + 2)
    rng = Xoshiro(semilla)
    mk = m ÷ k
    w = 4mk
    X = Int[]
    for c in rand(rng, 1:U-w, k)
        append!(X, c .+ randperm(rng, w)[1:mk])
    end
    sort!(unique!(X))
end

"""
    consultas_cola_pesada(U, m; α=1.5, semilla)

`m` huecos con distribución de cola pesada `u^(-1/α)`, reescalados al universo `1:U`.
"""
function consultas_cola_pesada(U, m; α=1.5, semilla=SEMILLA + 3)
    rng = Xoshiro(semilla)
    h = rand(rng, m) .^ (-1 / α)
    c = cumsum(h)
    unique!(ceil.(Int, U .* c ./ c[end]))
end

"""
    distancias(A, X)

Distancias `p` entre las posiciones de inserción de consultas consecutivas de `X` en `A`.
"""
function distancias(A, X)
    pos = searchsortedfirst.(Ref(A), X)
    diff(pos)
end

cuantil(P, q) = sort(P)[max(1, round(Int, q * length(P)))]

Note que las funciones reciben el universo \(U\) y no el arreglo, y que las ráfagas y la cola pesada pueden producir unas pocas consultas menos que \(m\) al eliminar duplicados. La siguiente salida sirve como control: antes de medir nada, reproduzca estos valores con sus propios datos para verificar que generó las mismas instancias. Si trabaja en Python, replique la lógica de los generadores; las secuencias aleatorias no serán idénticas pero las distribuciones de \(p\) deben ser muy similares a las de este cuadro.

A = genera_arreglo(2^20)
U = A[end]
println("n = $(length(A)), U = $U")
generadores = (("uniforme", consultas_uniformes), ("ráfagas", consultas_rafagas),
               ("cola pesada", consultas_cola_pesada))
for (nombre, gen) in generadores, m in (2^8, 2^14, 2^20)
    X = gen(U, m)
    P = distancias(A, X)
    etiqueta = rpad(nombre, 12) * rpad("m = 2^$(Int(log2(m)))", 10)
    etiqueta *= rpad("|X| = $(length(X))", 16)
    q = (cuantil(P, 0.5), cuantil(P, 0.9), cuantil(P, 0.99), maximum(P))
    println(etiqueta, "p: mediana = $(q[1]), q90 = $(q[2]), q99 = $(q[3]), máx = $(q[4])")
end

n = 1048576, U = 4721224
uniforme    m = 2^8   |X| = 256       p: mediana = 2715, q90 = 9422, q99 = 17499, máx = 20902
uniforme    m = 2^14  |X| = 16384     p: mediana = 44, q90 = 150, q99 = 295, máx = 633
uniforme    m = 2^20  |X| = 1048576   p: mediana = 1, q90 = 2, q99 = 5, máx = 15
ráfagas     m = 2^8   |X| = 256       p: mediana = 1, q90 = 3, q99 = 131250, máx = 176250
ráfagas     m = 2^14  |X| = 16384     p: mediana = 1, q90 = 2, q99 = 4, máx = 174156
ráfagas     m = 2^20  |X| = 1028146   p: mediana = 1, q90 = 2, q99 = 4, máx = 40102
cola pesada m = 2^8   |X| = 256       p: mediana = 2698, q90 = 7353, q99 = 18519, máx = 48708
cola pesada m = 2^14  |X| = 16384     p: mediana = 33, q90 = 96, q99 = 416, máx = 32400
cola pesada m = 2^20  |X| = 1048576   p: mediana = 1, q90 = 2, q99 = 7, máx = 5976

Con estos datos, realice lo siguiente:

Implementación. Implemente la búsqueda binaria acotada y las búsquedas no acotadas \(B_0\), \(B_1\) y \(B_2\), todas con el parámetro de pista sp. Verifique cada una contra searchsortedfirst sobre arreglos aleatorios con y sin duplicados y con distintas pistas válidas. Escriba una función buscar_todas(f, A, X) que resuelva en orden todas las consultas de X con el algoritmo f, comenzando cada búsqueda donde terminó la anterior; para la búsqueda binaria independiente la pista será siempre \(1\). Puede usar CountingArray (Listado 3) para contar comparaciones.
Datos. Genere el arreglo A = genera_arreglo(2^20) y, para cada \(m \in \{2^4, 2^6, 2^8, \ldots, 2^{20}\}\), los tres conjuntos de consultas con los generadores de Listado 6. Si su equipo no permite trabajar con \(n = 2^{20}\), reduzca \(n\) a la potencia de dos más grande que pueda y repórtelo; en ese caso sus valores de control serán distintos.
Caracterización de las instancias. Antes de medir, calcule para cada configuración las distancias \(p\) entre respuestas consecutivas con distancias y muestre su distribución, e.g., con histogramas en escala logarítmica o con los cuantiles \(0.5\), \(0.9\) y \(0.99\). Compare estas distancias con los umbrales \(\sqrt{n}\) y \(n / \log_2^2 n\) del capítulo.
Medición. Para cada algoritmo, tamaño \(m\) y distribución, mida el número total de comparaciones y el tiempo acumulado para resolver las \(m\) consultas; repita las mediciones de tiempo (e.g., con @benchmark) para obtener valores fiables. Reporte comparaciones y tiempo por separado, en gráficas con \(m\) en el eje horizontal (escala logarítmica) y una curva por algoritmo, una gráfica por distribución.
Puntos de cruce. A partir de sus gráficas, determine para cada distribución los valores de \(m/n\) donde: (a) \(B_1\) deja de superar a la búsqueda binaria independiente, (b) \(B_2\) deja de superarla, y (c) \(B_0\) se vuelve competitivo. Contraste con las predicciones del capítulo (\(p < \sqrt{n}\) para \(B_1\) y \(p < n/\log_2^2 n\) para \(B_2\)), tomando en cuenta la distribución de \(p\) que caracterizó en el punto 3.
Discusión. Explique las diferencias entre lo que indican las comparaciones y lo que indica el tiempo. Considere que los algoritmos no acotados realizan más operaciones aritméticas por comparación pero acceden a posiciones cercanas entre sí, mientras que las búsquedas binarias independientes saltan por todo el arreglo. ¿Qué algoritmo recomendaría para intersectar dos listas ordenadas de tamaños \(m\) y \(n\), y cómo depende la recomendación de \(m/n\)? Guarde sus implementaciones: las reutilizará en el Capítulo 6.

1.7.2 Entregable

Su trabajo se entregará en PDF y con el notebook fuente, se entregará en el sistema en archivos separados. El código se debe documentar, así como mantener una estructura del documento que permita a un lector interesado entender el problema, sus experimentos y metodología, así como sus conclusiones. Tenga en cuenta que los notebooks pueden alternar celdas de texto y código.

No olvide estructurar su reporte, en particular el reporte debe cubrir los siguientes puntos:

Título del reporte, su nombre.
Introducción.
Desarrollo de las actividades.
Conclusiones
Lista de referencias. Nota, una lista de referencias que no fueron utilizadas en el cuerpo del texto será interpretada como una lista vacía.

Nota: no olvides revisar la rúbrica del curso.

Referencias

Bentley, Jon L., y Catherine C. McGeoch. 1985. “Amortized analyses of self-organizing sequential search heuristics”. Commun. ACM (New York, NY, USA) 28 (4): 404–11. https://doi.org/10.1145/3341.3349.

Bentley, Jon Louis, y Andrew Chi-Chih Yao. 1976. “An almost optimal algorithm for unbounded searching”. Information processing letters 5 (SLAC-PUB-1679).

Cormen, Thomas H, Charles E Leiserson, Ronald L Rivest, y Clifford Stein. 2022. Introduction to algorithms. MIT press.

Knuth, Donald. 1998. The Art Of Computer Programming, vol. 3 (2nd ed): Sorting And Searching. Vol. 3. Addison Wesley Longman Publishing Co. Inc.

Perl, Yehoshua, Alon Itai, y Haim Avni. 1978. “Interpolation search—a log log N search”. Communications of the ACM 21 (7): 550–53. https://doi.org/10.1145/359545.359557.

Pugh, William. 1990. “Skip lists: a probabilistic alternative to balanced trees”. Commun. ACM (New York, NY, USA) 33 (6): 668–76. https://doi.org/10.1145/78973.78977.