Dies ist Ihr Haskell (nicht) nur faktoriell und gut für

Wenn es um Lieblingssprachen geht, sage ich normalerweise, dass ich bei gleichen Bedingungen C ++ für Zahlenbrecher und Haskell für alles andere bevorzuge. Es ist nützlich, regelmäßig zu überprüfen, ob diese Unterteilung gerechtfertigt ist, und erst kürzlich stellte sich eine müßige und sehr einfache Frage: Wie wird sich die Summe aller Teiler einer Zahl mit dem Wachstum dieser Zahl beispielsweise für die ersten Milliarden Zahlen verhalten? Es ist leicht, diese Aufgabe einzuschüchtern (es ist eine Schande, sie als den resultierenden Zahlenschleifer zu bezeichnen), daher scheint sie eine großartige Option für eine solche Prüfung zu sein.

Außerdem kann ich die Leistung des Haskell-Codes immer noch nicht genau vorhersagen. Daher ist es hilfreich, wissentlich schlechte Ansätze zu versuchen, um festzustellen, wie sich die Leistung verschlechtert.

Außerdem können Sie auf einfache Weise einen effizienteren Algorithmus als die frontale Suche nach Teilern für jede Zahl von vorführen $$$1$ vorher $$$n$ .

Algorithmus

Beginnen wir also mit dem Algorithmus.

So finden Sie die Summe aller Teiler einer Zahl $$$n$ ? Sie können alles durchgehen $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ und für jeden solchen $$$k_1$ Überprüfen Sie den Rest der Abteilung $$$n$ auf $$$k_1$ . Wenn der Rest ist $$$0$ , dann zur Batterie hinzufügen $$$k_1 + k_2$ wo $$$k_2 = \ frac {n} {k_1}$ wenn $$$k_1 \ neq k_2$ und nur $$$k_1$ sonst.

Kann dieser Algorithmus angewendet werden? $$$n$ mal für jede nummer von $$$1$ vorher $$$n$ ? Das kannst du natürlich. Was wird die Schwierigkeit sein? Diese Reihenfolge ist leicht zu erkennen $$$O (n ^ \ frac {3} {2})$ Divisionen - für jede Zahl machen wir genau die Wurzel der Divisionen, und wir haben Zahlen $$$n$ . Können wir es besser machen? Es stellt sich heraus, dass ja.

Eines der Probleme bei dieser Methode ist, dass wir zu viel Aufwand verschwenden. Zu viele Abteilungen führen uns nicht zum Erfolg und geben einen Rest ungleich Null. Es ist natürlich zu versuchen, etwas fauler zu sein und sich der Aufgabe von der anderen Seite zu nähern: Lassen Sie uns einfach alle Arten von Kandidaten für Teiler generieren und sehen, welche Zahlen sie erfüllen?

Also, jetzt brauchen wir einen Sturzflug für jede Nummer von $$$1$ vorher $$$n$ Berechnen Sie die Summe aller Teiler. Gehen Sie dazu alles durch $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ und für jeden solchen $$$k_1$ Lass uns alles durchgehen $$$k_2 \ in \ {k_1 \ dots \ lfloor \ frac {n} {k} \ rfloor \}$ . Für jedes Paar $$$(k_1, k_2)$ zur Zelle mit dem Index hinzufügen $$$k_1 \ cdot k_2$ Wert $$$k_1 + k_2$ wenn $$$k_1 \ neq k_2$ und $$$k_1$ sonst.

Dieser Algorithmus macht genau $$$n ^ \ frac {1} {2}$ Divisionen und jede Multiplikation (die billiger als Division ist) führt uns zum Erfolg: Bei jeder Iteration erhöhen wir etwas. Dies ist viel effektiver als der frontale Ansatz.

Darüber hinaus können Sie mit demselben Frontalansatz beide Implementierungen vergleichen und sicherstellen, dass sie für relativ kleine Zahlen dieselben Ergebnisse liefern, was ein wenig Vertrauen schaffen sollte.

Erste Implementierung

Übrigens ist dies direkt fast ein Pseudocode der ersten Implementierung in 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 Modul ist einfach und ich bringe es nicht mit.

Außerdem zeigen wir hier den Betrag nur für die meisten $$$n$ zum leichteren Vergleich mit anderen Implementierungen. Trotz der Tatsache, dass Haskell eine faule Sprache ist, werden in diesem Fall alle Beträge berechnet (obwohl die vollständige Rechtfertigung dafür den Rahmen dieses Artikels sprengt), sodass es nicht funktioniert, dass wir nichts versehentlich zählen.

Wie schnell geht es? Auf meinem i7 3930k werden in einem Stream 100.000 Elemente in 0,4 s ausgearbeitet. In diesem Fall werden 0,15 s für Berechnungen und 0,25 s für GC aufgewendet. Und wir belegen ungefähr 8 Megabyte Speicher, obwohl wir, da die Größe des int 8 Bytes beträgt, idealerweise 800 Kilobyte haben sollten.

Gut (nicht wirklich). Wie werden diese Zahlen mit steigenden Zahlen wachsen? Für 1'000'000 Elemente hat es ungefähr 7,5 Sekunden gearbeitet, drei Sekunden für das Rechnen und 4,5 Sekunden für das GC aufgewendet und 80 Megabyte (10-mal mehr als nötig) belegt. Und selbst wenn wir uns für eine Sekunde als Senior Java Software Developers ausgeben und mit der Optimierung von GC beginnen, werden wir das Bild nicht wesentlich ändern. Schade. Es scheint, dass wir niemals eine Milliarde Nummern warten werden und auch nicht in den Speicher gelangen werden: Auf meinem Computer befinden sich nur 64 Gigabyte RAM, und es werden ungefähr 80 benötigt, wenn sich der Trend fortsetzt.

Es scheint Zeit zu tun

C ++ Option

Lassen Sie uns versuchen, eine Vorstellung davon zu bekommen, was es sinnvoll ist, danach zu streben, und dafür werden wir den Code auf die Pluspunkte schreiben.

Nun, da wir bereits einen debuggten Algorithmus haben, ist alles einfach:

 #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, eine Million Elemente werden in 0,024 s verarbeitet und belegen die zugewiesenen 8 Megabyte Speicher. Milliarden - 155 Sekunden, 8 Gigabyte Speicher, wie erwartet. Autsch. Haskell ist nicht gut. Haskell muss rausgeworfen werden. Nur Fakultäten und Präpromorphismen drauf und schreiben! Oder nicht?

Zweite Option

Es ist offensichtlich nicht die klügste Entscheidung, alle generierten Daten über IntMap , das heißt in der Tat eine relativ gewöhnliche Karte - gelinde gesagt, dies ist dieselbe offensichtlich miese Option, die zu Beginn erwähnt wurde). Warum verwenden wir kein Array wie im C ++ - Code?

Versuchen wir mal:

 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 

Hier verwenden wir sofort die Unboxed-Version des Arrays, da Int recht einfach ist und wir keine Faulheit darin brauchen. Die Box-Version würde sich nur im arr Typ unterscheiden, sodass wir auch nicht an Redewendungen verlieren. Außerdem wird die Bindung für bound hier separat vorgenommen, jedoch nicht, weil der Compiler dumm ist und kein LICM ausführt, sondern weil Sie dann seinen Typ explizit angeben und eine Warnung des Compilers über die Standardeinstellung des floor Arguments vermeiden können.

0,045 s für eine Million Elemente (nur zweimal schlechter als die Pluspunkte!). 8 Megabyte Speicher, null Millisekunden in GC (!). Bei größeren Größen bleibt der Trend bestehen - etwa doppelt so langsam wie bei C ++ und bei gleichem Speicherplatz. Tolles Ergebnis! Aber können wir es besser machen?

Es stellt sich heraus, dass ja. accumArray überprüft die Indizes, was in diesem Fall nicht erforderlich ist - die Indizes sind korrekt aufgebaut. Versuchen wir, den Aufruf von accumArray durch unsafeAccumArray zu 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 

Wie Sie sehen können, sind die Änderungen minimal, mit Ausnahme der Notwendigkeit, von Grund auf neu indiziert zu werden (was meiner Meinung nach ein Fehler in der Bibliotheks-API ist, aber dies ist eine andere Frage). Was ist die Leistung?

Eine Million Elemente - 0,021 s (wow, innerhalb der Fehlergrenze, aber schneller als die Profis!). Natürlich die gleichen 8 Megabyte Speicher, die gleichen 0 ms im GC.

Milliarden Elemente - 152 s (es sieht so aus, als wäre es wirklich schneller als die Pluspunkte!). Etwas weniger als 8 Gigabyte. 0 ms in GC. Der Code ist immer noch idiomatisch. Ich denke, wir können sagen, dass dies ein Sieg ist.

Abschließend

Erstens war ich überrascht, dass das Ersetzen von accumArray durch eine unsafe Version zu einer solchen Steigerung führen wird. Es wäre vernünftiger, 10 bis 20 Prozent zu erwarten (schließlich führt das Ersetzen von operator[] durch at() nicht zu einer signifikanten Leistungsminderung), aber nicht um die Hälfte!

Zweitens ist es meiner Meinung nach wirklich cool, dass ein ziemlich idiomatischer sauberer Code ohne veränderliches Herausragen dieses Leistungsniveau erreicht.

Drittens sind natürlich weitere Optimierungen auf allen Ebenen möglich. Ich bin mir zum Beispiel sicher, dass Sie etwas mehr aus dem Code auf den Pluspunkten herausholen können. Meiner Meinung nach ist jedoch bei all diesen Benchmarks das Gleichgewicht zwischen dem Aufwand (und der Menge an Code) und dem daraus resultierenden Auspuff wichtig. Ansonsten konvergiert letztendlich alles mit der Herausforderung von LLVM JIT oder ähnlichem. Darüber hinaus gibt es sicherlich effizientere Algorithmen zur Lösung dieses Problems, aber das Ergebnis eines hier vorgestellten kurzen Gedankens wird auch für dieses kleine Sonntagsabenteuer funktionieren.

Viertens mein Favorit: Typsysteme müssen entwickelt werden. unsafe wird hier nicht benötigt, als Programmierer kann ich beweisen, dass k_1 * k_2 <= n für alle k_1, k_2 in der Schleife angetroffen wird. In einer idealen Welt von abhängig typisierten Sprachen würde ich diesen Beweis statisch konstruieren und an die entsprechende Funktion übergeben, wodurch die Notwendigkeit von Überprüfungen zur Laufzeit entfällt. Aber leider gibt es in Haskell kein vollwertiges Zavtipov, und in Sprachen, in denen Zavtipy ist (und die ich kenne), gibt es kein array und keine Analoga.

Fünftens kenne ich andere Programmiersprachen nicht genug, um mich für Benchmarks in diesen Sprachen zu qualifizieren, aber einer meiner Freunde hat ein Analogon in Python geschrieben. Fast genau hundertmal langsamer und schlimmer aus dem Gedächtnis. Und der Algorithmus selbst ist extrem einfach. Wenn also jemand mit Kenntnissen ein Analogon in Go, Rust, Julia, D, Java, Malbolge in die Kommentare schreibt oder einen Vergleich beispielsweise mit C ++ - Code auf seinem Computer teilt, ist er wahrscheinlich großartig .

PS: Entschuldigung für den leicht Clickbait-Header. Ich könnte mir nichts Besseres einfallen lassen.