Este es tu Haskell (no) solo factorial y bueno para

Cuando se trata de lenguajes favoritos, generalmente digo que, en igualdad de condiciones, prefiero C ++ para los trituradores de números y Haskell para todo lo demás. Es útil verificar periódicamente si esta división está justificada, y recientemente surgió una pregunta ociosa y muy simple: cómo se comportará la suma de todos los divisores de un número con el crecimiento de este mismo número, por ejemplo, para los primeros mil millones de números. Es fácil intimidar esta tarea (es una lástima llamarlo el generador de números resultante), por lo que parece una gran opción para tal verificación.

Además, todavía no tengo la capacidad de predecir con precisión el rendimiento del código Haskell, por lo que es útil intentar enfoques deliberadamente malos para ver cómo se degradará el rendimiento.

Bueno, además, puede mostrar fácilmente un algoritmo más eficiente que una búsqueda frontal de divisores para cada número de $$$1$ antes $$$n$ .

Algoritmo

Entonces, comencemos con el algoritmo.

Cómo encontrar la suma de todos los divisores de un número $$$n$ ? Puedes pasar por todo $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ y por cada tal $$$k_1$ verifica el resto de la división $$$n$ en $$$k_1$ . Si el resto es $$$0$ , luego agregue a la batería $$$k_1 + k_2$ donde $$$k_2 = \ frac {n} {k_1}$ si $$$k_1 \ neq k_2$ y solo $$$k_1$ de lo contrario

¿Se puede aplicar este algoritmo? $$$n$ veces, para cada número de $$$1$ antes $$$n$ ? Puedes, por supuesto. ¿Cuál será la dificultad? Fácil de ver ese orden $$$O (n ^ \ frac {3} {2})$ divisiones: para cada número hacemos exactamente la raíz de las divisiones, y tenemos números $$$n$ . ¿Podemos hacerlo mejor? Resulta que sí.

Uno de los problemas con este método es que estamos desperdiciando demasiado esfuerzo. Demasiadas divisiones no nos llevan al éxito, dando un resto distinto de cero. Es natural tratar de ser un poco más perezoso y abordar la tarea desde el otro lado: ¿vamos a generar todo tipo de candidatos para divisores y ver qué números satisfacen?

Entonces, ahora necesitamos una caída repentina para cada número de $$$1$ antes $$$n$ calcular la suma de todos sus divisores. Para hacer esto, revisa todo $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ y para cada uno $$$k_1$ repasemos todo $$$k_2 \ in \ {k_1 \ dots \ lfloor \ frac {n} {k} \ rfloor \}$ . Para cada par $$$(k_1, k_2)$ agregar a la celda con el índice $$$k_1 \ cdot k_2$ valor $$$k_1 + k_2$ si $$$k_1 \ neq k_2$ y $$$k_1$ de lo contrario

Este algoritmo hace exactamente $$$n ^ \ frac {1} {2}$ divisiones, y cada multiplicación (que es más barata que la división) nos lleva al éxito: en cada iteración aumentamos algo. Esto es mucho más efectivo que el enfoque frontal.

Además, con este mismo enfoque frontal, puede comparar ambas implementaciones y asegurarse de que den los mismos resultados para números bastante pequeños, lo que debería agregar un poco de confianza.

Primera implementación

Y, por cierto, esto es directamente casi un pseudocódigo de la implementación inicial en Haskell:

module Divisors.Multi(divisorSums) where import Data.IntMap.Strict as IM divisorSums :: Int -> Int divisorSums n = IM.fromListWith (+) premap IM.! n where premap = [ (k1 * k2, if k1 /= k2 then k1 + k2 else k1) | k1 <- [ 1 .. floor $ sqrt $ fromIntegral n ] , k2 <- [ k1 .. n `quot` k1 ] ] 

Main módulo es simple y no lo traigo.

Además, aquí mostramos la cantidad solo para la mayoría $$$n$ para facilitar la comparación con otras implementaciones. A pesar del hecho de que Haskell es un lenguaje perezoso, en este caso se calcularán todas las cantidades (aunque la justificación completa para esto está más allá del alcance de este artículo), por lo que no resulta que no contaremos nada inadvertidamente.

¿Qué tan rápido funciona? En mi i7 3930k, en una secuencia, 100,000 elementos se resuelven en 0.4 s. En este caso, se gastan 0,15 s en cálculos y 0,25 s en GC. Y ocupamos aproximadamente 8 megabytes de memoria, aunque, dado que el tamaño del int es de 8 bytes, idealmente deberíamos tener 800 kilobytes.

Bueno (no realmente) ¿Cómo crecerán estos números con números crecientes? Para 1'000'000 elementos, ha estado trabajando durante aproximadamente 7.5 segundos, gastando tres segundos en computación y 4.5 segundos gastando en GC, y también ocupando 80 megabytes (10 veces más de lo necesario). E incluso si pretendemos ser Desarrolladores de Software Java Senior por un segundo y comenzamos a ajustar GC, no cambiaremos significativamente la imagen. Que mal. Parece que nunca esperaremos mil millones de números, y tampoco entraremos en la memoria: solo hay 64 gigabytes de RAM en mi máquina, y tomará alrededor de 80 si la tendencia continúa.

Parece tiempo de hacer

Opción C ++

Tratemos de tener una idea de lo que tiene sentido luchar, y para esto escribiremos el código en las ventajas.

Bueno, como ya tenemos un algoritmo depurado, entonces todo es simple:

 #include <vector> #include <string> #include <cmath> #include <iostream> int main(int argc, char **argv) { if (argc != 2) { std::cerr << "Usage: " << argv[0] << " maxN" << std::endl; return 1; } int64_t n = std::stoi(argv[1]); std::vector<int64_t> arr; arr.resize(n + 1); for (int64_t k1 = 1; k1 <= static_cast<int64_t>(std::sqrt(n)); ++k1) { for (int64_t k2 = k1; k2 <= n / k1; ++k2) { auto val = k1 != k2 ? k1 + k2 : k1; arr[k1 * k2] += val; } } std::cout << arr.back() << std::endl; } 

-O3 -march=native , clang 8, se procesan un millón de elementos en 0.024 s, ocupando los 8 megabytes de memoria asignados. Mil millones: 155 segundos, 8 gigabytes de memoria, como se esperaba. Ouch Haskell no es bueno. Haskell necesita ser expulsado. ¡Solo factoriales y prepromorfismos en él y escribir! O no?

Segunda opción

Obviamente, ejecutar todos los datos generados a través de IntMap , que es, de hecho, un mapa relativamente ordinario, por decirlo suavemente, no es la decisión más sabia (sí, esta es la misma opción obviamente mala que se mencionó al principio). ¿Por qué no usamos una matriz como en el código C ++?

Probemos

 module Divisors.Multi(divisorSums) where import qualified Data.Array.IArray as A import qualified Data.Array.Unboxed as A divisorSums :: Int -> Int divisorSums n = arr A.! n where arr = A.accumArray (+) 0 (1, n) premap :: A.UArray Int Int premap = [ (k1 * k2, if k1 /= k2 then k1 + k2 else k1) | k1 <- [ 1 .. floor bound ] , k2 <- [ k1 .. n `quot` k1 ] ] bound = sqrt $ fromIntegral n :: Double 

Aquí usamos inmediatamente la versión sin caja de la matriz, ya que Int bastante simple y no necesitamos pereza. La versión en caja solo diferiría en el tipo arr , por lo que no perdemos en el idioma también. Además, el enlace para el bound se realiza por separado aquí, pero no porque el compilador sea estúpido y no haga LICM, sino porque entonces puede especificar explícitamente su tipo y evitar una advertencia del compilador sobre el incumplimiento del argumento de floor .

0.045 s para un millón de elementos (¡solo dos veces peor que las ventajas!). 8 megabytes de memoria, cero milisegundos en GC (!). En tamaños más grandes, la tendencia persiste: aproximadamente dos veces más lenta que C ++ y la misma cantidad de memoria. Gran resultado! ¿Pero podemos hacerlo mejor?

Resulta que sí. accumArray verifica los índices, lo cual no es necesario que hagamos en este caso: los índices son correctos en la construcción. Intentemos reemplazar la llamada a accumArray con unsafeAccumArray :

 module Divisors.Multi(divisorSums) where import qualified Data.Array.Base as A import qualified Data.Array.IArray as A import qualified Data.Array.Unboxed as A divisorSums :: Int -> Int divisorSums n = arr A.! (n - 1) where arr = A.unsafeAccumArray (+) 0 (0, n - 1) premap :: A.UArray Int Int premap = [ (k1 * k2 - 1, if k1 /= k2 then k1 + k2 else k1) | k1 <- [ 1 .. floor bound ] , k2 <- [ k1 .. n `quot` k1 ] ] bound = sqrt $ fromIntegral n :: Double 

Como puede ver, los cambios son mínimos, excepto por la necesidad de ser indexados desde cero (que, en mi opinión, es un error en la API de la biblioteca, pero esta es otra pregunta). ¿Cuál es el rendimiento?

Un millón de elementos: 0.021 s (¡guau, dentro del margen de error, pero más rápido que los profesionales!). Naturalmente, los mismos 8 megabytes de memoria, los mismos 0 ms en el GC.

Mil millones de elementos: 152 s (¡parece que es realmente más rápido que las ventajas!). Ligeramente menos de 8 gigabytes. 0 ms en GC. El código sigue siendo idiomático. Creo que podemos decir que esto es una victoria.

En conclusión

En primer lugar, me sorprendió que el reemplazo de accumArray con una versión unsafe aumentaría tanto. Sería más razonable esperar 10-20 por ciento (después de todo, en las ventajas, reemplazar el operator[] con at() no da una disminución significativa en el rendimiento), ¡pero no a la mitad!

En segundo lugar, en mi opinión, es realmente genial que un código limpio bastante idiomático sin un mutable sobresaliente alcance este nivel de rendimiento.

En tercer lugar, por supuesto, son posibles otras optimizaciones, y en todos los niveles. Estoy seguro, por ejemplo, de que puede extraer un poco más del código de las ventajas. Sin embargo, en mi opinión, en todos los puntos de referencia, el equilibrio entre el esfuerzo realizado (y la cantidad de código) y el escape resultante es importante. De lo contrario, todo finalmente converge al desafío de LLVM JIT o algo así. Además, existen ciertamente algoritmos más eficientes para resolver este problema, pero el resultado de una breve reflexión presentada aquí también funcionará para esta pequeña aventura dominical.

Cuarto, mi favorito: los sistemas de tipos necesitan ser desarrollados. unsafe no es necesario aquí, como programador puedo demostrar que k_1 * k_2 <= n para todos los k_1, k_2 encontrados en el bucle. En un mundo ideal de lenguajes de tipo dependiente, construiría esta prueba estáticamente y la pasaría a la función correspondiente, lo que eliminaría la necesidad de verificaciones en tiempo de ejecución. Pero, por desgracia, en Haskell no hay zavtipov completo, y en los idiomas donde hay zavtipy (y que sé), no hay una array y análogos.

Quinto, no conozco otros lenguajes de programación lo suficiente como para calificar para puntos de referencia cercanos en estos idiomas, pero uno de mis amigos escribió un análogo en Python. Casi exactamente cien veces más lento, y peor de memoria. Y el algoritmo en sí mismo es extremadamente simple, por lo que si alguien conocedor escribe un análogo en Go, Rust, Julia, D, Java, Malbolge en los comentarios, o comparte una comparación, por ejemplo, con el código C ++ en su máquina, probablemente será genial .

PD: Perdón por el encabezado ligeramente clickbait. No se me ocurrió nada mejor.