Votre Haskell est (non) factoriel uniquement et bon pour

En ce qui concerne les langages préférés, je dis généralement que, toutes choses étant égales par ailleurs, je préfère C ++ pour les concasseurs de nombres et Haskell pour tout le reste. Il est utile de vérifier périodiquement si cette division est justifiée, et ce n'est que récemment qu'une question oiseuse et très simple s'est posée: comment la somme de tous les diviseurs d'un nombre se comportera-t-elle avec la croissance de ce nombre, disons, pour le premier milliard de nombres. Il est facile d'intimider cette tâche (il est dommage de l'appeler le broyeur de numéros résultant), donc cela semble être une excellente option pour un tel contrôle.

De plus, je n'ai toujours pas la possibilité de prédire avec précision les performances du code Haskell, il est donc utile d'essayer des approches sciemment mauvaises pour voir comment les performances se dégraderont.

De plus, vous pouvez facilement montrer un algorithme plus efficace que la recherche frontale de diviseurs pour chaque nombre de $$$1$ avant $$$n$ .

Algorithme

Commençons donc par l'algorithme.

Comment trouver la somme de tous les diviseurs d'un nombre $$$n$ ? Vous pouvez parcourir tout $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ et pour chaque $$$k_1$ vérifier le reste de la division $$$n$ sur $$$k_1$ . Si le reste est $$$0$ , puis ajoutez à la batterie $$$k_1 + k_2$ où $$$k_2 = \ frac {n} {k_1}$ si $$$k_1 \ neq k_2$ , et juste $$$k_1$ sinon.

Cet algorithme peut-il être appliqué $$$n$ fois, pour chaque numéro de $$$1$ avant $$$n$ ? Vous pouvez bien sûr. Quelle sera la difficulté? Facile à voir cette commande $$$O (n ^ \ frac {3} {2})$ divisions - pour chaque nombre, nous faisons exactement la racine des divisions, et nous avons des nombres $$$n$ . Pouvons-nous faire mieux? Il s'avère que oui.

L'un des problèmes de cette méthode est que nous gaspillons trop d'efforts. Trop de divisions ne nous mènent pas au succès, donnant un reste non nul. Il est naturel d’essayer d’être un peu plus paresseux et d’aborder la tâche de l’autre côté: générons simplement toutes sortes de candidats pour les diviseurs et voyons quels nombres ils satisfont?

Alors, maintenant, il nous faut un seul coup pour chaque numéro de $$$1$ avant $$$n$ calculer la somme de tous ses diviseurs. Pour ce faire, parcourez tous $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ , et pour chacun de ces $$$k_1$ passons par tout $$$k_2 \ in \ {k_1 \ dots \ lfloor \ frac {n} {k} \ rfloor \}$ . Pour chaque paire $$$(k_1, k_2)$ ajouter à la cellule avec l'index $$$k_1 \ cdot k_2$ valeur $$$k_1 + k_2$ si $$$k_1 \ neq k_2$ et $$$k_1$ sinon.

Cet algorithme fait exactement $$$n ^ \ frac {1} {2}$ les divisions, et chaque multiplication (qui est moins chère que la division) nous mène au succès: à chaque itération nous augmentons quelque chose. C'est beaucoup plus efficace que l'approche frontale.

De plus, ayant cette même approche frontale, vous pouvez comparer les deux implémentations et vous assurer qu'elles donnent les mêmes résultats pour des nombres assez petits, ce qui devrait ajouter un peu de confiance.

Première implémentation

Et en passant, c'est directement presque un pseudo-code de l'implémentation initiale dans 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 module est simple, et je ne l'apporte pas.

De plus, nous montrons ici le montant uniquement pour les $$$n$ pour faciliter la comparaison avec d'autres implémentations. Malgré le fait que Haskell soit une langue paresseuse, dans ce cas, tous les montants seront calculés (bien que la justification complète de cela dépasse le cadre de cet article), donc cela ne fonctionne pas que nous ne comptons rien par inadvertance.

À quelle vitesse ça marche? Sur mon i7 3930k, en un seul flux, 100 000 éléments sont élaborés en 0,4 s. Dans ce cas, 0,15 s est consacré aux calculs et 0,25 s au GC. Et nous occupons environ 8 mégaoctets de mémoire, bien que, comme la taille de l'int est de 8 octets, nous devrions idéalement avoir 800 kilo-octets.

Bon (pas vraiment). Comment ces nombres vont-ils croître avec l'augmentation, euh, des nombres? Pour 1'000'000 éléments, il fonctionne depuis environ 7,5 secondes, consacre trois secondes à l'informatique et 4,5 secondes au GC, et occupe 80 mégaoctets (10 fois plus que nécessaire). Et même si nous prétendons être des développeurs de logiciels Java principaux pendant une seconde et commençons à régler GC, nous ne changerons pas de manière significative l'image. Dommage. Il semble que nous n'attendrons jamais un milliard de chiffres, et nous n'entrerons pas non plus dans la mémoire: il n'y a que 64 gigaoctets de RAM sur ma machine, et cela prendra environ 80 si la tendance se poursuit.

Il semble temps de faire

Option C ++

Essayons de nous faire une idée de ce qu'il est logique de rechercher, et pour cela, nous allons écrire le code sur les avantages.

Eh bien, puisque nous avons déjà un algorithme débogué, alors tout est 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, un million d'éléments sont traités en 0,024 s, occupant les 8 mégaoctets de mémoire alloués. Milliard - 155 secondes, 8 gigaoctets de mémoire, comme prévu. Ouch. Haskell n'est pas bon. Haskell doit être expulsé. Seuls les factoriels et les prépromorphismes dessus et écrivent! Ou pas?

Deuxième option

Évidemment, exécuter toutes les données générées via IntMap , c'est-à-dire, en fait, une carte relativement ordinaire - pour le dire légèrement, n'est pas la décision la plus sage (oui, c'est la même option évidemment moche qui a été mentionnée au début). Pourquoi n'utilisons-nous pas un tableau comme dans le code C ++?

Essayons:

 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 

Ici, nous utilisons immédiatement la version sans boîte du tableau, car Int assez simple, et nous n'avons pas besoin de paresse. La version en boîte ne différerait que par le type arr , donc nous ne perdons pas non plus en idiome. En outre, la liaison pour bound est effectuée séparément ici, mais pas parce que le compilateur est stupide et ne fait pas de LICM, mais parce qu'alors vous pouvez spécifier explicitement son type et éviter un avertissement du compilateur concernant la défaillance de l'argument floor .

0,045 s pour un million d'éléments (seulement deux fois pire que les avantages!). 8 mégaoctets de mémoire, zéro millisecondes en GC (!). À des tailles plus grandes, la tendance persiste - environ deux fois plus lente que C ++, et la même quantité de mémoire. Excellent résultat! Mais pouvons-nous faire mieux?

Il s'avère que oui. accumArray vérifie les indices, ce que nous n'avons pas besoin de faire dans ce cas - les indices sont corrects dans la construction. Essayons de remplacer l'appel à accumArray par 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 

Comme vous pouvez le voir, les changements sont minimes, à l'exception de la nécessité d'être indexé à partir de zéro (ce qui, à mon avis, est un bogue dans l'API de la bibliothèque, mais c'est une autre question). Quelle est la performance?

Un million d'éléments - 0,021 s (wow, dans la marge d'erreur, mais plus vite que les pros!). Naturellement, les mêmes 8 mégaoctets de mémoire, les mêmes 0 ms dans le GC.

Milliards d'éléments - 152 s (il semble que c'est vraiment plus rapide que les avantages!). Un peu moins de 8 gigaoctets. 0 ms en GC. Le code est toujours idiomatique. Je pense que nous pouvons dire que c'est une victoire.

En conclusion

Tout d'abord, j'ai été surpris de constater que le remplacement de accumArray par une version unsafe entraînerait une telle augmentation. Il serait plus raisonnable de s'attendre à 10 à 20% (après tout, dans les avantages, remplacer l' operator[] par at() ne donne pas une baisse significative des performances), mais pas de moitié!

Deuxièmement, à mon avis, c'est vraiment cool qu'un code propre assez idiomatique sans dépassement mutable atteigne ce niveau de performance.

Troisièmement, bien entendu, de nouvelles optimisations sont possibles, à tous les niveaux. Je suis sûr, par exemple, que vous pouvez extraire un peu plus du code sur les avantages. Cependant, à mon avis, dans tous ces repères, l'équilibre entre l'effort dépensé (et la quantité de code) et l'échappement qui en résulte est important. Sinon, tout converge finalement vers le défi de LLVM JIT ou quelque chose comme ça. En outre, il existe certainement des algorithmes plus efficaces pour résoudre ce problème, mais le résultat d'une courte réflexion présentée ici fonctionnera également pour cette petite aventure du dimanche.

Quatrièmement, mon préféré: les systèmes de type doivent être développés. unsafe n'est pas nécessaire ici, en tant que programmeur, je peux prouver que k_1 * k_2 <= n pour tous les k_1, k_2 rencontrés dans la boucle. Dans un monde idéal de langages typés de manière dépendante, je construirais cette preuve statiquement et la transmettrais à la fonction correspondante, ce qui supprimerait le besoin de vérifications lors de l'exécution. Mais, hélas, à Haskell, il n'y a pas de zavtipov à part entière, et dans les langues où se trouvent les zavtipy (et que je connais), il n'y a pas de array et d'analogues.

Cinquièmement, je ne connais pas suffisamment d'autres langages de programmation pour pouvoir prétendre à des quasi-benchmarks dans ces langages, mais un de mes amis a écrit un analogue en python. Presque exactement cent fois plus lentement et pire de mémoire. Et l'algorithme lui-même est extrêmement simple, donc si quelqu'un compétent écrit un analogue dans Go, Rust, Julia, D, Java, Malbolge dans les commentaires, ou partage une comparaison, par exemple, avec C ++ - du code sur leur machine - ce sera probablement génial .

PS: Désolé pour l'en-tête légèrement clickbait. Je ne pourrais rien trouver de mieux.