Este seu Haskell é (não) apenas fatorial e bom para

Quando se trata de idiomas favoritos, costumo dizer que, sendo todas as coisas iguais, prefiro C ++ para trituradores de números e Haskell para todo o resto. É útil verificar periodicamente se essa divisão é justificada, e apenas recentemente surgiu uma questão ociosa e muito simples: como a soma de todos os divisores de um número se comportará com o crescimento desse número, digamos, para o primeiro bilhão de números. É fácil intimidar essa tarefa (é uma pena chamá-la de multiplicador de números resultante), por isso parece uma ótima opção para essa verificação.

Além disso, ainda não tenho a capacidade de prever com precisão o desempenho do código Haskell, por isso é útil tentar abordagens conscientemente ruins para ver como o desempenho diminui.

Além disso, você pode exibir facilmente um algoritmo mais eficiente do que a busca frontal por divisores para cada número de $$$1$ antes $$$n$ .

Algoritmo

Então, vamos começar com o algoritmo.

Como encontrar a soma de todos os divisores de um número $$$n$ ? Você pode passar por todos $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ e para todos esses $$$k_1$ verifique o restante da divisão $$$n$ em $$$k_1$ . Se o restante for $$$0$ e adicione à bateria $$$k_1 + k_2$ onde $$$k_2 = \ frac {n} {k_1}$ se $$$k_1 \ neq k_2$ e apenas $$$k_1$ caso contrário.

Esse algoritmo pode ser aplicado $$$n$ vezes, para cada número de $$$1$ antes $$$n$ ? Você pode, é claro. Qual será a dificuldade? Fácil de ver esse pedido $$$O (n ^ \ frac {3} {2})$ divisões - para cada número, fazemos exatamente a raiz das divisões e temos números $$$n$ . Podemos fazer melhor? Acontece que sim.

Um dos problemas com esse método é que estamos desperdiçando muito esforço. Muitas divisões não nos levam ao sucesso, dando um restante diferente de zero. É natural tentar ser um pouco mais preguiçoso e abordar a tarefa do outro lado: vamos gerar todos os tipos de candidatos a divisores e ver quais números eles satisfazem?

Então, agora vamos precisar de um golpe para cada número de $$$1$ antes $$$n$ calcular a soma de todos os seus divisores. Para fazer isso, passe por todos $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ e para cada um desses $$$k_1$ vamos passar por tudo $$$k_2 \ in \ {k_1 \ dots \ lfloor \ frac {n} {k} \ rfloor \}$ . Para cada par $$$(k_1, k_2)$ adicionar à célula com o índice $$$k_1 \ cdot k_2$ valor $$$k_1 + k_2$ se $$$k_1 \ neq k_2$ e $$$k_1$ caso contrário.

Esse algoritmo faz exatamente $$$n ^ \ frac {1} {2}$ divisões e cada multiplicação (que é mais barata que a divisão) nos leva ao sucesso: a cada iteração aumentamos algo. Isso é muito mais eficaz que a abordagem frontal.

Além disso, com essa mesma abordagem frontal, você pode comparar as duas implementações e garantir que elas apresentem os mesmos resultados para números razoavelmente pequenos, o que deve adicionar um pouco de confiança.

Primeira implementação

E, a propósito, isso é diretamente quase um pseudocódigo da implementação inicial em 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 é simples e não o trago.

Além disso, aqui mostramos o valor apenas para os mais $$$n$ para facilitar a comparação com outras implementações. Apesar de Haskell ser uma linguagem preguiçosa, nesse caso, todos os valores serão calculados (embora a justificativa completa esteja além do escopo deste artigo), portanto, não resulta que não contaremos nada inadvertidamente.

Quão rápido ele funciona? No meu i7 3930k, em um fluxo, 100.000 elementos são trabalhados em 0,4 s. Nesse caso, 0,15 s são gastos em cálculos e 0,25 s em GC. E ocupamos cerca de 8 megabytes de memória, embora, como o tamanho do int seja 8 bytes, o ideal seja termos 800 kilobytes.

Bom (não realmente). Como esses números crescerão com números cada vez maiores? Para 1'000'000 elementos, ele trabalha há cerca de 7,5 segundos, gastando três segundos em computação e 4,5 segundos em GC, além de ocupar 80 megabytes (10 vezes mais que o necessário). E mesmo se fingirmos ser desenvolvedores sênior de software Java por um segundo e começarmos a ajustar o GC, não mudaremos significativamente o cenário. Que pena. Parece que nunca esperaremos um bilhão de números e também não entraremos na memória: há apenas 64 gigabytes de RAM na minha máquina e levará cerca de 80 se a tendência continuar.

Parece hora de fazer

Opção C ++

Vamos tentar ter uma idéia do que faz sentido nos esforçarmos e, para isso, escreveremos o código nas vantagens.

Bem, como já temos um algoritmo depurado, tudo é simples:

 #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 , cl 8, um milhão de elementos é processado em 0,024 s, ocupando os 8 megabytes de memória alocados. Bilhões - 155 segundos, 8 gigabytes de memória, conforme o esperado. Ai. Haskell não é bom. Haskell precisa ser jogado fora. Apenas fatoriais e prepromorfismos nele e escreva! Ou não?

Segunda opção

Obviamente, executar todos os dados gerados através do IntMap , que é, de fato, um mapa relativamente comum - para dizer o mínimo, não é a decisão mais sensata (sim, esta é a mesma opção obviamente péssima mencionada no início). Por que não usamos uma matriz como no código C ++?

Vamos tentar:

 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 

Aqui, imediatamente usamos a versão sem caixa da matriz, pois Int bastante simples e não precisamos de preguiça. A versão em caixa diferiria apenas no tipo arr , por isso não perdemos o idioma também. Além disso, a ligação para bound é feita separadamente aqui, mas não porque o compilador seja estúpido e não faz o LICM, mas porque então você pode especificar explicitamente seu tipo e evitar um aviso do compilador sobre o padrão do argumento de floor .

0,045 s para um milhão de elementos (apenas duas vezes pior que as vantagens!). 8 megabytes de memória, zero milissegundos em GC (!). Em tamanhos maiores, a tendência persiste - cerca de duas vezes mais lenta que o C ++ e a mesma quantidade de memória. Ótimo resultado! Mas podemos fazer melhor?

Acontece que sim. accumArray verifica os índices, o que não precisamos fazer neste caso - os índices estão corretos na construção. Vamos tentar substituir a chamada para accumArray por 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 você pode ver, as alterações são mínimas, exceto pela necessidade de serem indexadas do zero (o que, na minha opinião, é um bug na API da biblioteca, mas essa é outra questão). Qual é o desempenho?

Um milhão de elementos - 0,021 s (uau, dentro da margem de erro, mas mais rápido que os profissionais!). Naturalmente, os mesmos 8 megabytes de memória, os mesmos 0 ms no GC.

Bilhões de elementos - 152 s (parece que é realmente mais rápido que os benefícios!). Um pouco menos de 8 gigabytes. 0 ms no GC. O código ainda é idiomático. Eu acho que podemos dizer que isso é uma vitória.

Em conclusão

Em primeiro lugar, fiquei surpreso que a substituição do accumArray por uma versão unsafe trará esse aumento. Seria mais razoável esperar de 10 a 20% (afinal, nos trunfos, substituir o operator[] por at() não causa uma diminuição significativa no desempenho), mas não pela metade!

Em segundo lugar, na minha opinião, é muito legal que um código limpo bastante idiomático, sem uma alteração mutável, atinja esse nível de desempenho.

Em terceiro lugar, é claro, outras otimizações são possíveis e em todos os níveis. Estou certo, por exemplo, de que você pode extrair um pouco mais do código das vantagens. No entanto, na minha opinião, em todos esses benchmarks, o equilíbrio entre o esforço despendido (e a quantidade de código) e o escape resultante é importante. Caso contrário, tudo acabará convergindo para o desafio do LLVM JIT ou algo assim. Além disso, certamente existem algoritmos mais eficientes para resolver esse problema, mas o resultado de uma breve reflexão apresentada aqui também funcionará para esta pequena aventura de domingo.

Quarto, meu favorito: os sistemas de tipos precisam ser desenvolvidos. unsafe não é necessário aqui, como programador posso provar que k_1 * k_2 <= n para todos os k_1, k_2 encontrados no loop. Em um mundo ideal de linguagens de tipo dependente, eu construíva essa prova estaticamente e a passava para a função correspondente, o que eliminaria a necessidade de verificações em tempo de execução. Mas, infelizmente, em Haskell não há zavtipov de pleno direito, e em idiomas em que há zavtipy (e que eu sei), não há array e análogos.

Em quinto lugar, não conheço outras linguagens de programação suficientes para se qualificar para quase benchmarks nessas linguagens, mas um dos meus amigos escreveu um análogo em python. Quase exatamente cem vezes mais devagar e pior da memória. E o algoritmo em si é extremamente simples; portanto, se alguém com conhecimento escrever um analógico em Go, Rust, Julia, D, Java, Malbolge nos comentários, ou compartilhar uma comparação, por exemplo, com o código C ++ em sua máquina - provavelmente será ótimo .

PS: Desculpe pelo cabeçalho ligeiramente isca de clique. Eu não consegui pensar em nada melhor.