Fakultät


Schnelle Faktorielle - große Fakultätsberechnungen mit einem ATmega1284
10.6.2019

Anhand der relativ einfachen, aber sehr stack- und rechenintensiven Ackermann-Funktion wurde bereits gezeigt, dass strukturiertes AVR-Assembler-Programmieren mit s'AVR nicht nur Spaß machen kann, sondern auch noch eine sehr gute Performanz ergibt.

Nun, rekursive Funktionen sind eher etwas für Mathematiker, Informatiker und Programmierer, denn man kann sich solche Funktionen (insbesondere wie im Falle Ackermann) nicht immer einfach vorstellen.

Deshalb soll nun ein weiteres, deutlich komplexeres Beispiel schließlich ohne Rekursion folgen, das zum Bestimmen sehr großer Fakultäten einen ATmega1284 als Rechenknecht einsetzt, der mit s’AVR besonders beeindruckt.

Die Fakultät (teilweise auch Faktorielle genannt, im Englischen Factorial) einer Zahl z ergibt sich aus der Multiplikation aller natürlichen Zahlen von 1 bis z, also 1 * 2 * 3 * ... * z, die mit steigendem z sehr schnell große Werte annimmt. Verkürzt schreibt man für dieses Produkt einfach "z!".

Je nach verfügbarem (Taschen-) Rechner kann man per Fakultät-Taste/Funktion folgende Falkultäten gerade noch ganzzahlig als Dezimalzahl (also ohne 10er-Exponent) berechnen und darstellen lassen:

11! = 39916800‬ (8 Digits),
12! = 479001600 (9 Digits),
13! = 6227020800‬ (10 Digits, Kc-S56, SR-50, TI-30, TI-83 und sehr viele mehr),
14! = 87178291200 (11 Digits, Excel, Google direkt im Suchfenster)
17! = 355687428096000 (15 Digits, wissenschaftlicher Galaxy-Rechner)
18! = 6402373705728000‬ (16 Digits, wissenschaftlicher iPhone-Rechner) oder gar
29! = 8841761993739701954543616000000 (31 Digits, wissenschaftlicher Windows-Rechner)

Der wissenschaftliche Rechner namens SpeedCrunch (verfügbar für Windows, Linux und macOS) mit teilweise sehr ausgefallenen Funktionen und mit sehr vielen Stellen Genauigkeit zeigt per Gamma-Funktion Γ(z+1) = z! maximal Γ(58) = 57! mit 77 Stellen und 13 Endnullen ganzzahlig an.

Der Online-Rechner JavaScripter schreckt schließlich auch vor 12000! mit 43742 Stellen (und 2998 Endnullen) nicht zurück. Bei 12001! ist aber auch hier Schluss.

Die ganzzahlige Darstellung ist unter anderem dann hilfreich, wenn man wirklich alle Dezimalstellen sucht oder auch nur die Anzahl der Endnullen einer großen Fakultät bestimmen will - letzteres eine kleine Denksportaufgabe aus meiner Schulzeit (die man ggf. in kurzer Zeit sogar für 1000! leicht durch Kopfrechnen lösen kann).

Auch rekursiv

Nun soll nicht verschwiegen werden, dass man eine Fakultät auch rekursiv berechnen kann, nämlich anhand der Definition z! = z * (z-1)! mit 1! = 0! = 1, was man sich (im Unterschied zur Ackermann-Funktion) sogar sehr gut vorstellen kann.

Auch programmiertechnisch wird das sehr einfach, z. B. in Java oder C-Sprache wie folgt (man beachte, dass "fakul" sich selbst aufruft = rekursiv):

    static long fakul(long z)
    {
       if (z>1) return (z * fakul(z-1));
       else return(1);
    }

Es gibt davon viele Varianten und natürlich auch C-typische Verstümmelungen ...

Die Fakultät in Pascal/Delphi als rekursive Funktion formuliert, könnte so aussehen:

    function fakul(z: int64): int64;
    begin
       if z = 1 then fakul := 1
       else fakul := z * fakul(z-1);
    end;

Zum rekursiven Berechnen einer Fakultät fängt man also mit dem höchsten Wert z an und multipliziert dann immer weiter mit der nächst kleineren Ganzzahl bis herunter zu 1.

Im Unterschied zur Ackermann-Funktion steigt die Rekursionstiefe stetig an, um dann bei Erreichen von z=1 rasch wieder stetig abgebaut zu werden.

Bei dieser rekursiven Methode wird man sehr schnell (nämlich für z>20) den zulässigen Wertebereich von "long" bzw. "int64" überschreiten. Das Maximum mit der rekursiven Methode wäre bei den genannten Programmbeispielen also 20! = 2432902008176640000. Das sind immerhin 19 Dezimalstellen (und 4 Endnullen).

Die klassische Methode

Möchte man noch größere Fakultäten exakt berechnen, z.B. 100! oder gar 1000!, wird man mit der oben gezeigten rekursiven Methode nicht weiter kommen und muss für die Darstellung aller Dezimalstellen zur klassischen Multiplikation "von Hand" greifen ("schriftliches Multiplizieren"), wie man es in der Grundschule gelernt hat. D. h. man erstellt pro Zehnerpotenz des jeweiligen z (z. B. 4-stellig, angefangen mit 1) ein eigenes (mit größer werdendem z sehr langes) Teilprodukt und addiert am Schluss alle Teilprodukte zum Zwischenergebnis, das man dann mit dem nächst größeren z multipliziert, u.s.w.

Aber schnell wird man erkennen, dass das Handrechnen bei so großen z-Werten recht mühsam ist, denn allein zum schnellen Niederschreiben des bekannten Ergebnisses von 1000! mit 2568 Stellen wäre man etwa 20 Minuten beschäftigt (bei kleiner Handschrift könnte eine DIN-A-4-Seite damit ziemlich voll werden). Und selbst wenn man Tag und Nacht rechnen würde, wird ein Monat für das schriftliche Berechnen von 1000! kaum reichen.

Man nimmt hierfür also besser einen schnellen Rechenknecht!

Ein moderner PC würde (je nach verwendeter Programmiersprache) das Ergebnis von 1000! in deutlich weniger als einer Zehntelsekunde berechnet haben. Mit einem 8-bit-Mikrocontroller dauert es aber schon etwas länger - vorausgesetzt, er hat auch genügend Arbeitsspeicher.

Und genau das wollen wir untersuchen, z.B. anhand eines ATmega1284, jenem 8-bit-AVR mit einem sehr großen Arbeitsspeicher[6] (nämlich 16 kByte), den es auch noch in handhabbaren Gehäusen gibt.

Diese "händische" Methode hatte ich auf dem PC bereits in den 90er Jahren mit QBasic und QuickC unter DOS versucht und gleich danach mit Pascal (Delphi 5) unter Windows, womit die Rechenzeiten drastisch kürzer wurden.

Letzteres ist heute noch meine Referenz, aber die einfache QuickC-Version für DOS soll mit einigen Erweiterungen zunächst als Basis für GCC und den ATmega1284 dienen, der bei meinen Untersuchungen mit 18 MHz getaktet wird.

Und so schaut das Quellprogramm in GCC aus (Kopie aus Atmel Studio), das fast beliebig große Fakuläten berechnet (bis der Speicher ausgeht):

[Hinweis: Auf Smartphones werden teilweise leider Inlineframes nicht richtig dargestellt. Deshalb erscheint dort der Quellcode in voller Länge und nicht in einem kleineren Scroll-Fenster.]
 

Wie im GCC-Quellprogramm am Ende erwähnt, benötigt der ATmega1284 zur Berechnung von 1000! in der ursprünglichen Version stolze 2:16 Minuten (mehr dazu siehe Ergänzung vom 30.6.2019).

Zur Berechnung anderer Fakultäten (z. B. 1227!, das Maximum für GCC auf dem ATmega1284) müssen im Quellprogramm bei Bedarf zwei Änderungen gemacht werden, nämlich unter Main der Wert für z:

und für optimale Speicherausnützung und kürzeste Laufzeit ganz am Anfang:

Falls ExpLen unnötig groß gewählt wird, hat das nur den Nachteil längerer Laufzeiten oder im schlimmsten Fall stürzt das Programm ab (der Compiler streikt bei ExpLen>3269).

Das Hauptprogramm läuft in einer Endlosschleife. Die Laufzeit ist jene Zeit vom Erlöschen der LED bis zum nächsten Aufleuchten der LED.

Da bei meinen GCC-Untersuchungen eine einfache Platine mit ausschließlich dem ATmega1284 verwendet wurde und kein Display angeschlossen war, wurden zur Kontrolle nur die ersten 4 Digits des Ergebnisses per

und

ausgegeben und an den Pins nachgemessen (Low bzw. High per Voltmeter, Scope oder Logic-Analyzer). Die Werte "40" und "23" sind für 1000! ein sicheres Indiz, dass die Berechnung aller 2568 Stellen stimmt. Bei anderen Werten (insbesondere "00") ist etwas schief gelaufen.

Natürlich kann man je nach Lust & Laune und entsprechendem Hardware-Aufwand auch das per ATmega1284 berechnete Ergebnis auf irgend einem Medium in seiner voller Länge einschließlich der vielen Endnullen ausgeben lassen. Per Laufschrift auf einer längeren 7-Segment-LED-Zeile (mit 500 msec Schrittzeit) wäre vielleicht eine Idee, um die eingangs erwähnten ca. 20 Minuten für eine Abschrift zu überprüfen ...

Tweaks

Die Abfragen von zd[j] und fakul[k] auf die Werte 0 und 1 bringen bereits eine große Steigerung bezüglich Rechengeschwindigkeit, da u. a. viele Multiplikationen entfallen.

Überraschenderweise hat aber der größere Aufwand, die indizierte FOR-Schleife

zur Berechnung der Teilprodukte zu strecken ("Unrolling", "Unwinding") überhaupt keinen Laufzeitgewinn gebracht, so dass diese Methode im obigen GCC-Programm schließlich doch nicht angewendet wurde, wodurch das Programm natürlich deutlich kürzer und etwas übersichtlicher wird.

Aber: Der GCC-Compiler hat genau dieselbe Idee und wendet bei dieser FOR-Schleife zum Optimieren das Unrolling standardmäßig an - und deshalb gibt es keinen Unterschied zu meinem eigenen Versuch (den ich mir hätte schenken können, wenn ich die LSS-Datei des Compilers gleich angeschaut hätte)!

Strukturiert in AVR-Assembler programmiert

Jetzt geht es natürlich noch darum, die Ergebnisse von GCC mit s’AVR zu vergleichen.

Um wie beim GCC-Programm ebenfalls einen großen Laufzeitgewinn zu erzielen, muss man die Methode des "Entrollens" beim s'AVR-Programm also selbst umsetzen.

Das eigentliche s’AVR-Quellprogramm kann man entweder vom GCC-Quellprogramm ableiten, aber wegen der ähnlichen Syntax noch einfacher auch von einem Basic- oder Pascal-Programm (so man hat).

Nachfolgend schließlich das per Copy&Paste 4-fach gestreckte (und wie immer strikt strukturiert programmierte) s'AVR-Quellprogramm, das durch das Strecken zwar deutlich länger ist im Vergleich zu einer indizierten FOR-Schleife, aber dafür erstaunliche zehnmal schneller als das GCC-Programm unter sonst gleichen Bedingungen:

Softwaretricks

Es ist aber nicht nur das Unrolling der genannten FOR-Schleife, sondern auch das trickreiche Anordnen und Berechnen der Teilprodukte, u. a. mit einem speziell hierfür entworfenen Macro MOD_DIV10, das die Einerziffern und den Übertrag in die nächste Dezimalstelle extrem schnell (in nur 10 Maschinenzyklen) per MOD 10 und DIV 10 einer 8-Bit-Zahl (0..255) bestimmt, das im GCC-Programm den beiden Zeilen

entspricht.

Hier das gesamte MOD_DIV10-Macro nebst ausführlicher Beschreibung (vielleicht kann das jemand für eigene AVR-Programme gebrauchen):

Die Konstante 205 ist die größtmögliche 8-bit-Zahl (nämlich für maximale Genauigkeit), die sich durch Multiplikation von 0,1 mit einer 2er-Potenz ergibt, sprich 0,1 * 2^11 = 204,8, aufgerundet auf 205.

Damit ergibt sich 205/2^11 = 0,100097656, d. h. 1 ‰ bzw. 1/4 Bit Genauigkeit bei einem 8-bit-Wert. Dieses Macro bestimmt MOD10/DIV10 demnach exakt für alle 8-bit-Zahlen.

Interessant ist in diesem Zusammenhang vielleicht noch die Feststellung, dass die stellenweise Summe der vier Teilprodukte[1] bei 1000! nie größer als 29 ist.

Dieses (universell verwendbare) Macro belegt nur 9 Worte Programmspeicher pro Aufruf bzw. wegen 9 Aufrufen insgesamt 81 Worte für das gesamte s'AVR-Programm.

Im s'AVR-Hauptprogramm ist AKKU im oberen AVR-Registerbereich definiert, CARRY darf auch im unteren AVR-Registerbereich sein.

Programmspeicherbedarf

Das GCC-Programm benötigt für die reine Fakultätberechnung 728 Bytes, zzgl. 76 Bytes für das Unterprogramm (divmodhi4) für die MOD10/DIV10-Berechnung. Das ist zwar etwas sparsamer als ein vielfach verwendetes Macro, aber bei häufigen Aufrufen deutlich langsamer.

Beim s'AVR-Programm sind es für die reine Fakultätberechnung 330 Bytes inklusive der neun MOD_DIV10-Macros für die Fakultätberechnung.

Insgesamt ist das GCC-Programm (nur die eigentliche Fakultätberechnung!) also mehr als doppelt so groß wie das s'AVR-Programm.

Beim s'AVR-Programm kann man die zu berechnende Fakultät in einer Definitionsdatei zwischen 100!, 500!, 1000! und 1230! auswählen und die diversen Programmvariablen werden vom Assembler automatisch angepasst.

Außerdem werden beim s'AVR-Programm die ersten 8 Digits des Ergebnisses, die Anzahl der Ergebnis-Digits und die Zahl der Endnullen (beides per Programm aus dem Ergebnis ausgezählt) auf einem TM1640-basierenden[2] 16-Digit-LED-Display dezimal dargestellt, hier für 1000!:

Factorial1000

Auch hier läuft das Hauptprogramm nach dem ersten Berechnen und Anzeigen in einer Endlosschleife zum Messen der Laufzeit, wiederum ab Erlöschen der LED bis zum nächsten Aufleuchten der LED.

Um den Arbeitsspeicher maximal zu nützen (nämlich 99,8%), kann man mit dem s'AVR-Programm bis zu 1230! mit 3269 Stellen in ca. 21 Sekunden berechnen.

Beim GCC-Programm ist aufgrund eines größeren Stack-Bedarfs nur maximal 1227! mit 3260 Stellen (und 304 Endnullen) möglich (Laufzeit 3:19 Minuten).

Fazit s’AVR vs. GCC

Im Vergleich zum GCC-Programm zur Berechnung von sehr großen Fakultäten (z. B. 1000!) auf einem ATmega1284 ist das vorgestellte s'AVR-Programm zirka zehnmal schneller und benötigt dabei nur etwa 41% des Programmspeichers. Bei beiden Programmen wird dabei derselbe Berechnungsalgorithmus angewendet.

Der C-Compiler verwendet (grundsätzlich) ungeschickterweise für die Konstante 0[3] das AVR-Register R1, das aber (neben R0) bei vielen AVR-CPUs auch für die Multiplikation verwendet wird. Das Retten bzw. Neubeschreiben von R1 mit dem Wert 0 kostet unnötig Zeit.

Bei vergleichsweise ähnlich guter strukturierter Programmierung wie bei GCC kann man per s’AVR aber alle Spezialitäten der AVR-CPU berücksichtigen und damit einige Fallstricke umgehen.

So habe ich zwar für schnelle Zugriffe und optimale Ausnützung der oberen AVR-Register ebenfalls einige Konstante in den weniger flexiblen unteren AVR-Registern abgelegt, aber R0 und R1 natürlich für die Multiplikation freigehalten, die im vorgestellten Programm sehr häufig[4] aufgerufen wird.

Auch das ausgefeilte MOD_DIV10-Macro (das ebenfalls eine Multiplikation enthält) trägt zum Laufzeitgewinn bei.

Und schließlich erlaubt das im s’AVR-Programm angewandte Verschachteln von Zwischenergebnis und den vier Teilprodukten in einem einzigen großen Array die elegante Adressierung aller Arrays gleichzeitig mit nur einem Pointer-Register und per Displacement, das bei AVR leider nur einen kleinen Bereich von 0 bis 63 zulässt und deshalb bei fünf separaten Arrays für das sehr lange Zwischenergebnis und die ebenso langen vier Teilprodukte nicht angewendet werden könnte.

Bestimmt gibt es beim vorgestellten GCC-Programm noch ein paar Details zu verbessern, um kürzere Laufzeiten zu erzielen. Aber eine Steigerung um den Faktor 10 scheint mir nicht machbar.

Überraschend ist vielleicht noch die Tatsache, dass der gestreckte s’AVR-Quellcode für die reine Fakultät-Berechnung nur knapp über dreimal so viele Programmzeilen benötigt wie der ungestreckte GCC-Quellcode. So betrachtet, wäre das s’AVR-Quellprogramm sogar noch etwas kompakter als das GCC-Quellprogramm.

Ziemliche Beschleunigung durch ein kleines Detail (30.6.2019)

Nun schien der Laufzeitunterschied zwischen GCC und s’AVR doch etwas sehr groß, so dass ich das ursprünglich für den PC geschriebene GCC-Programm nochmals genauer angeschaut habe.

Und siehe da, eine kleine Änderung bei der Definition zweier Variablen ergibt einen sehr großen Schritt nach vorn, auch wenn beim Ergebnis immer noch ein deutlicher Abstand zu s’AVR übrig bleibt (und nun der Laufzeitvergleich ähnlich wie bei der Ackermann-Funktion ausfällt):

// int a, b;             // Zwischenwerte bei der Fakultät-Berechnung
uint8_t a, b;            // Zwischenwerte bei der Fakultät-Berechnung

Die großen Arrays fakul[ExpLen+3] und tp[4][ExpLen+3] wurden von Anfang an als 8-Bit-Variable definiert, denn sonst wäre die Berechnung von 1000! aufgrund des begrenzten Arbeitsspeichers des ATmega1284 überhaupt nicht möglich gewesen.

Überraschend, dass der GCC-Compiler bei den Variablen a und b bei der ursprünglichen Definition "int a, b;" keine Optimierung gefunden hat.

Mit dieser Änderung wird vom GCC-Programm deutlich mehr Programmspeicher belegt (nämlich nun 946 Bytes für die reine Fakultätberechnung), da statt dem Unterprogramm divmodhi4 bis auf ein einziges Mal nun (wie bei s’AVR) immer eine offensichtlich Macro-basierende Befehlsfolge eingefügt wird.

Schließlich beträgt die GCC-Laufzeit für 1000! jetzt ca. 25 sec und ist bei 2,9-fachem Programmspeicherbedarf somit "nur" noch um den Faktor 1,8 langsamer als das s’AVR-Programm.

Algorithmus verbessern? (10.6.2019)

Den Berechnungsalgorithmus könnte man für alle Programme noch insoweit verbessern, dass man die Teilprodukte nicht immer über die volle Länge berechnet, sondern - ähnlich wie beim Rechnen mit Papier und Bleistift - mit einem dynamischen Pointer beim Zwischenergebnis nur bis zur höchsten Ziffer <> 0.

Da aber bei den Teilprodukten bereits grundsätzlich auf die Ziffern 0 und 1 abgefragt wird, bringt eine solche Vorgehensweise vielleicht nicht allzuviel Laufzeitgewinn.

Kein großer Gewinn (16.6.2019)

Inzwischen habe ich mit meinem Delphi-Programm ein Verfahren untersucht, das das Berechnen der Teilprodukte abkürzt, sobald eine führende Null beim Zwischenergebnis erreicht wird, mit dem Ergebnis, dass es (bedingt durch zusätzliche Abfragen) nur dann einen signifikanten Laufzeitgewinn gibt, wenn die Arrays deutlich größer sind als für die jeweilige Fakultät unbedingt nötig wäre.

Deshalb werde ich mir diesen extra Aufwand bei den AVR-Programmen sparen, da die Programme dadurch umfangreicher werden und nicht mehr so übersichtlich sind.

Nur knapp schneller mit Delphi und i5-CPU

Das erwähnte, aber hier nicht im Detail beschriebene Pascal/Delphi-Programm gibt nicht nur das Ergebnis von bis zu 2000! (bei Bedarf und Neukompilierung sogar bis 9999! [5]) in voller Länge auf dem Bildschirm aus, sondern bestimmt auch die Anzahl der Ergebnisstellen und der Endnullen, ebenso die Laufzeit.

Fakultät

Dieses Delphi-Programm benötigt mit demselben Berechnungsalgorithmus für 1000! bei optimierten Array-Größen (oder einem dynamischem Pointer, siehe Anmerkung vom 16.6.2019) und einer i5-CPU @3,3 GHz nur 64 msec (sonst ca. 110 msec).

Umgerechnet auf die 18 MHz des ATmega1284 dürfte dieser für 1000! demnach 11,7 sec benötigen.

So gesehen sind die mit dem s'AVR-Programm erzielten 14 Sekunden schon ziemlich nahe dran, das GCC-Programm ist mit 2:18 Minuten (per int a, b;) bzw. 25 Sekunden (per uint8_t a, b;) Laufzeit für dieselbe Berechnung aber deutlich abgeschlagen.

In diesem Fall der "händischen" Berechnung einer großen Fakultät gewinnt eben die 64-bit-CPU knapp gegen den kleinen 8-bit-Mikrocontroller, der sich aber bei geschickter Programmierung trotzdem ganz wacker schlägt.

Bei Interesse kann das Delphi-Programm für Windows (WIN98SE bis WIN10) gerne als EXE- oder ZIP-Datei zur Verfügung gestellt werden.

Die Profis

Wenn man richtig viel Rechen-Power nebst großem Arbeitsspeicher und geeignete Mathe-Libraries zur Verfügung hat, lassen sich (extrem) große Fakultäten mit ziemlich ausgefeilten Methoden auch noch schneller berechnen. Diese sind aber für einen kleinen 8-bit-Mikrocontroller nicht mehr geeignet.

nach oben


[1]  Bei 2000! ist die stellenweise Summe der vier Teilprodukte nie größer als 39, bei 100! nie größer als 19 und bei 10! sogar nie größer als 8.

Lediglich beim Multiplizieren von Ziffern <>0 und <>1 und anschließendem Addieren des Carrys gibt es etwas größere Ergebnisse, allerdings (mindestens) bis 2000! (und das erst ab 59!) nie größer als 89 (von 9*9 + 8), so dass ein 8-bit-Register sowohl zum Multiplizieren mit Carry, zum Aufaddieren der Teilprodukte als auch für MOD_DIV10 allemal ausreicht.

[2] Zur Ansteuerung wird eine modifizierte TM1638-Library verwendet.

[3] Bei MIPS und anderen CPUs gibt es ein festverdrahtetes 0-Register (bei MIPS ist es R0).

[4] Wenn ich mich nicht verzählt habe, gibt es bei der Berechnung von 1000! trotz der Abfragen von zd[j] und fakul[k] auf Werte < 2 immer noch 2177871 reine Multiplikationen. Durch die Berechnung von MOD10 und DIV10 bei den Summenbildungen kommt ein Vielfaches davon dazu, so dass es insgesamt mehr als 10 Millionen Multiplikationen sein dürften.

[5] Der wissenschaftliche Windows-Rechner kann immerhin bis zu 3248! = 1,9736342530860425312047080034031e+9997 berechnen, allerdings nur in Exponentialdarstellung.

Beim wissenschaftlichen iPhone-Rechner ist bereits bei 103! = 9,90290071649e163 das Maximum erreicht und bei wissenschaftlichen Rechnern mit 8-stelliger Mantisse ist es meist nur 69! = 1,7112245e98.

[6] Genau genommen gibt es noch einen ATmega256RFR2 im VQFN-64-Gehäuse, der sogar 32 kByte SRAM (und einen 2,4-GHz-Transceiver) integriert hat. Ein solches Gehäuse ist für Hobyisten kaum brauchbar, wenn es nicht bereits auf einer Platine verlötet ist.