Intervall-Arithmetik auf Fix-Punkt-Zahlentyp

Am Beispiel des Komplements der Mandelbrotmenge z^2+c habe ich mal einen Zahlentyp (Fixpunkt) implementiert, der mittels Intervall-Arithmetik Rundungsfehler auffängt und so ein mathematisch verifiziertes Äußeres (bzw, genauer, eine Untermenge) berechnet.

Es war interessant zu sehen, dass der CPC nur etwa 3 Tage gebraucht hat, um ein 50x50-Bild zu berechnen: links mit c als komplexe Zahl, d.h. die Werte liegen auf einem Raster mit Abstand zueinander, und es ist unklar, was dazwischen geschieht, aber das Bild ähnelt schon entfernt den bekannten Darstellungen. Ein " "-Zeichen bedeutet, definitiv außerhalb, ein "X", aktuell unklar.

Benoit Mandelbrot hat ja seine erste grafische Darstellung Anfang der 1980er berechnet. Allerdings hat er bei IBM gearbeitet und wohl einen der damaligen Supercomputer verwenden können.

Das rechte Bild basiert auf c als komplexes Intervall, das die Fläche pflastert. Es ist in dieser Größe noch wenig aussagekräftig, dafür aber mathematisch akkurat ("ohne Lücken"😉. Eine C++-Version für größere Auflösung: https://github.com/marcm200/upperbound-rc. Dies ist aktuell aber mit dem Schneider nicht machbar.

Das angehängte Archiv hier enthält den Quellcode, eine kurze Beschreibung des Zahlentyps, der auf Tripel von Variablen basiert sowie eine Auflistung der Intervall- Arithmetik-Routinen, so dass man diese auch für andere Projekte verwenden kann. Die Routinen sind nur bedingt optimiert worden, da 3-4 Tage (geschätzt nach einem Lauf ohne Echtzeitsimulation) akzeptabel waren.
📎 _ps1_ps0.gif
📎 _retromset.zip

— marcm200
Das rechte Bild basiert auf c als komplexes Intervall, das die Fläche pflastert. Es ist in dieser Größe noch wenig aussagekräftig, dafür aber mathematisch akkurat ("ohne Lücken"😉. Eine C++-Version für größere Auflösung: https://github.com/marcm200/upperbound-rc. Dies ist aktuell aber mit dem Schneider nicht machbar.

Warum nicht? Woran liegt es? Speicher? Code nicht optimiert?
In MC sollte sich da doch noch viel rausholen lassen mMn.

Der Zeitbedarf wäre mir aktuell zu groß. Bei 200x200 könnte es theoretisch 16 mal so lange dauern, und da man mit höherer Auflösung näher an die Grenze zur Mandelbrotmenge kommt, müsste man auch eine höhere Maximaliterationszahl statt 15 nehmen, außerdem die Anzahl der Stellen des Zahlentyps (aktuell 8 ) auch graduell erhöhen. Bei 3 Tagen für das 50x50 Bild wären das bis zu 7 Wochen.

Optimierungspotential besteht eindeutig:

Auf Low-level könnte man versuchen, die vielen Variablen-Kopier-Aktionen zu minimieren oder evtl. direkt in Maschinensprache zu implementieren (quasi als memcopy, aber da fehlt mit die Erfahrung). Ich simuliere hier im wesentlichen C++-Referenzparameter und erstelle das BASIC-Programm durch C++-Makros.

Auf hoher Ebene bspw. mehr Fallunterscheidungen bei der Intervallmultiplikation, eine eigene Routine für die Quadrierung, besserer Rechnungsfluss (binomische Formel). Vielleicht auch Subdivision, aber da bin ich mir nicht sicher, ob das am Ende einen Zeitgewinn erbringt.

Außerdem wäre affine Arithmetik generell besser geeignet, da sie bereits bei kleinen Auflösungen schärfere Ergebnisse liefert. Das war mir aber in BASIC (noch) zu aufwendig, nachdem die C++-Version sehr kompliziert war.

Aber vielleicht hat ja bereits jemand so etwas programmiert?

Statt des Fix-Punkt-Zahlentyps könnte man auch den eingebauten float-Typ bei kleinen Auflösungen verwenden, wenn es die Möglichkeit gibt, den Rundungsmodus des CPC auf bspw. "Round up" oder "Round down" zu setzen.

(Ich habe gerade erst angefangen, die Literatur zu durchforsten): Weiß jemand, welchen Modus der CPC verwendet ("to nearest"?) und ob dieser geändert werden kann?

int() -> abrunden
round() -> nearest

Übrigens verstehe ich das Zeitproblem (3 Tage für so eine winzige Fläche??) gerade nicht.
Es gibt ja schon eine ganze Menge schneller Mandelbrot-Implementierungen, siehe z.B.:

https://www.octoate.de/2012/09/13/mavdelbrot-a-fast-mandelbrot-set-calculation-by-mav/

"The Amstrad CPC will need about 30 seconds to render the image."
Das ist noch nichtmal die schnellste Implementierung, im englischen CPCWiki-Forum wurde da noch mindestens eine schnellere vorgestellt und besprochen (müßte man mal raussuchen).

Ich habe mich allerdings nie in das Thema eingearbeitet. Vielleicht ist Deine Problemstellung eine ganz andere?

Mir geht es um mathematische Exaktheit, "Reliable Computing". Schnelle Berechnung eines Bildes, das aussieht wie eine Mandelbrotmenge, ist, wie Du sagst, recht einfach. Das Problem dort sind aber Rundungsfehler bei float-Berechnungen, wenn deren Ergebnis nicht mehr in float passen würde, d.h. "passend gemacht" wird.

Wenn ein Orbit den Kreis mit Radius 2 um den komplexen Nullpunkt verlässt, dann gehört der Startpunkt nicht zur Mandelbrotmenge ("escape time"😉. Nur, wenn man mit Rundung arbeiten muss, dann ist nicht sicher, dass dieses Verlassen auch wirklich das echte Verhalten ist. Der Punkt könnte gerade so durch Akkumulation von Rundungsfehlern aus dem Kreis "geflogen" sein.

Und dem entgegne ich mit Intervall-Arithmetik. Eine Berechnung, die nicht mehr  in den Datentyp passt, wird nach kleineren und größeren Werten so gerundet, dass  ein Intervall entsteht, das definitiv den echten Wert enthält. Und das bringt - durch die Verwendung von string$ als Basis-Datentyp, einen sehr hohen Zeitbedarf mit sich, da eine Multiplikation dann idR etwa 4 String$-Multiplikation bedarf. Dafür ist die Genauigkeit aber auch (fast) beliebig erhöhbar.

Könnte man den float-Rundungsmodus des CPC verändern, könnte man den float-Typ zumindest bis zu einer bestimmten Größe als Intervall-Grenztyp verwenden, ich müsste also nicht mehr mit string$ arbeiten.

Es ist hier also das alte Geben-und-Nehmen: Exaktheit wird erkauft durch Geschwindigkeitsverlust.

Ich mußte jetzt spontan an Mandelbrots Artikel "Wie lang ist die Küste Britanniens?" denken, und daß man mit Exaktheit alles bis ins Unendliche treiben kann, aber der Vergleich ist wohl nicht ganz richtig 🙂

— marcm200
Könnte man den float-Rundungsmodus des CPC verändern, könnte man den float-Typ zumindest bis zu einer bestimmten Größe als Intervall-Grenztyp verwenden, ich müsste also nicht mehr mit string$ arbeiten.

Wie genau muß er denn geändert werden?
round()
hilft Dir nicht? Bzw. wie genau sollte es sich verhalten?

— prodatron

— marcm200
Könnte man den float-Rundungsmodus des CPC verändern, könnte man den float-Typ zumindest bis zu einer bestimmten Größe als Intervall-Grenztyp verwenden, ich müsste also nicht mehr mit string$ arbeiten.

Wie genau muß er denn geändert werden?
round()
hilft Dir nicht? Bzw. wie genau sollte es sich verhalten?

round() liefert mir ja den abgerundeten Integer-Teil des float-Wertes. Der Rundungsmodus, den ich ändern möchte, bezieht sich aber auf das Ergebnis von float-Arithmetik.

Bspw. wenn man "double precision" in C++ nimmt, sind das, sagen wir, 17 dezimale Nachkommastellen. Wenn man nun 2 Mantissen multipliziert, erhält man oft als virtuelles Ergebnis eine Zahl mit deutlich mehr als 17 Nachkommastellen. Man muss diese also irgendwie in double "reinquetschen" (ignorieren wir mal den Exponenten).  Ähnliches gilt, wenn man einen sehr großen Wert und einen sehr kleinen addieren möchte, da dann meist der kleinere ignoriert werden muss (2^n + 2^-n für genügend große positive n)

Da spielt der Rundungsmodus rein. Man kann die überzähligen Ziffern einfach abschneiden (trunc), oder bspw. immer aufwärts runden, bis man bei max. 17 Nachkommastellen ist.

Da gibt es nun einen intelligenten Trick, durch wiederholtes Negieren (wenn dies exakt ist und nur auf dem Vorzeichenbit arbeitet), bei bspw. konstantem Aufwärtsrunden, Abwärtsrunden zu simulieren - und dann hätte man Intervall-Arithmetik auf double schnell implementiert. Denn das Umschalten des Rundungsmodus auf modernen Rechnern (für die Intervall-Untergrenze bräuchte man "round downard" und für die Obergrenze "round upward"😉 führt dazu, dass zuerst die FPU-Warteschlangen abgearbeitet werden müssen, der Prozessor kann also nicht die nächste Rechnung schon vorbereiten.

Falls es Dich interessiert, der zugrundeliegende Artikel ist: https://www.jstage.jst.go.jp/article/nolta/7/3/7_362/_article

round() könnte ich verwenden, das würde dann allerdings sehr breite Intervalle erzeugen, die dann keine definitive Aussage "liegt außerhalb der Mandelbrotmenge" mehr zulassen. IA ist um so besser, je enger die Intervalle sind.

— marcm200round() liefert mir ja den abgerundeten Integer-Teil des float-Wertes.

Wie gesagt, es rundet auf den "nearest", rundet also ab .5 auch auf:
round(2.4) = 2
round(2.5) = 3

— marcm200
Der Rundungsmodus, den ich ändern möchte, bezieht sich aber auf das Ergebnis von float-Arithmetik. [...] Man muss diese also irgendwie in double "reinquetschen" (ignorieren wir mal den Exponenten).

Das ist ne gute frage, wie er intern rundet, wenn das Ergebnis in eine niedrigere Auflösung gequetscht werden muss. Gibt es vielleicht ein Rechenbeispiel, mit dem man das mal testen kann, wenn man berücksichtigt, daß der CPC mit 32bit Mantissen arbeitet?

Irgendwie verstehe ich aber das grundsätzliche Problem nicht richtig:
Ein Computer ist "normal" sowieso nicht in der Lage, jede reele Zahl abzubilden, egal wie hochauflösend die Floats sind.
Beim CPC sind sie 40bits, bei modernen CPUs glaube ich 128bits oder so. D.h. es muß so oder so oft irgendwie was weggeschnitten werden.
Ob da jetzt beim letzten Bit richtig auf oder abgerundet wird, ändert ja nichts an dem Prinzip, daß eben trotzdem ein kleines bischen Genauigkeit wegfällt.
Ich bin allerdings auch kein Mathematiker.

— prodatron
Wie gesagt, es rundet auf den "nearest", rundet also ab .5 auch auf:

round(), int(), ja, die habe ich verwechselt. Habe sie seit fast Jahrzehnten nicht mehr verwendet, ebensowenig wie DEFINT, auch wieder-erlernt 🙂

— prodatron
Gibt es vielleicht ein Rechenbeispiel, mit dem man das mal testen kann, wenn man berücksichtigt, daß der CPC mit 32bit Mantissen arbeitet?

Ich habe aktuell nur ein C++-Programm (mein wiedererlerntes BASIC reicht dazu noch nicht aus).

void bsprd(void) {
    int32_t L=50;
    double a=1.0,b=1.0;
    for(int32_t k=0;k<L;k++) {
        a=a/2.0;
        b=b*2.0;
    }
    // d.h. a=2^-50 = (exakt lt wolfram-alpha.com) 0.00000000000000088817841970012523233890533447265625
    // b=2^50
    double c=a+b;

    printf("\na=%.40lg\nb=%.40lg\na+b=%.40lg\n",
        a,b,c);

    uint8_t *pa=(uint8_t*)&a;
    uint8_t *pb=(uint8_t*)&b;
    uint8_t *pc=(uint8_t*)&c;

    #define OUT8(TT,PP,NN) \
    {\
        printf("%16s= ",TT);\
        for(int32_t k=0;k<(NN);k++) {\
            int32_t w=(PP)[k];\
            printf("%3i ",w);\
        }\
        printf("\n");\
    }

    // Addition: c=b, d.h. a verschwindet
    OUT8("a",pa,8)
    OUT8("b",pb,8)
    OUT8("c=a+b",pc,8)

    // Multiplikation
    double d=a + 1.0;
    // d.h. 1+2^-50 = 1.00000000000000088817841970012523233890533447265625
    // lt wolfram-alpha.com

    double e=d*d; // d.h. a*a+2*a+1 = 2^-100+2^-49+1
    // exakt         : 1.0000000000000017763568394002512535387158899571179117285652827862...
    //                                                *
    // Ausgabe printf: 1.000000000000001776356839400250464677811
    printf("\nd=%.40lg\ne=d*d=%.40lg\n",d,e);
}

Addition:
2^50+2^-50 liefert als Ergebnis 2^50

Die 8 Bytes des double-Wertes als reiner Speicherinhalt. Man sieht, b und a+b sind bit-identisch. Damit wird a "vergessen".

               a=   0   0   0   0   0   0 208  60
               b=   0   0   0   0   0   0  16  67
           c=a+b=   0   0   0   0   0   0  16  67

Multiplikation:
d=(1+2^-50)^2

d=1.000000000000000888178419700125232338905
e=d*d=1.000000000000001776356839400250464677811
                                     *
exakt 1.0000000000000017763568394002512535387158899571179117285652827862...
(laut wolfram-alpha.com)

Hier wird allerdings auch die printf()-Konvertierungsroutine zur Ausgabe verwendet, die ebenfalls runden könnte. Aber eine bessere Ausgabe ist mir hier noch nicht eingefallen.

Irgendwie verstehe ich aber das grundsätzliche Problem nicht richtig:
Ein Computer ist "normal" sowieso nicht in der Lage, jede reele Zahl abzubilden, egal wie hochauflösend die Floats sind.

Ja, Computer können nur einige rationale Zahlen darstellen mit den float-Typen. Aber für normale Intervall-Arithmetik reicht das, die Grenzen dort sind immer repräsentierbare Zahlen. Die "echte" Zahl, bspw. wenn man sqrt(2) darstellen will, liegt dann mathematisch bewiesen immer innerhalb des Intervalls, auch nach bestimmten arithmetischen Operationen (das erinnert mich an den Matheunterricht in der Schule: Intervallschachtelung mit innerem Wert).

In Deinen beiden Beispielen wird jeweils die Mantisse zu lang, als daß es in 32 bit passen würde, weshalb beim CPC jeweils der "hintere" Teil verschwindet.
2^50 ist in dezimaler Darstellung 16stellig, eine 32bit Mantisse schafft halt nur maximal 10stellige Dezimalstellen.

Da aber gleich mehrere Bits hinten abgeschnitten werden, und in dem Bereich zwischendurch sowieso alles 0 ist, dürfte es doch egal sein, ob er beim "reinquetschen" auf- oder abrundet, hier wäre ein Abrunden eh korrekt.
Daher wäre es für Dich hier (in dem Beispiel jedenfalls) nicht relevant, das Runden zu beeinflussen?

— prodatron
In Deinen beiden Beispielen wird jeweils die Mantisse zu lang, als daß es in 32 bit passen würde, weshalb beim CPC jeweils der "hintere" Teil verschwindet.
2^50 ist in dezimaler Darstellung 16stellig, eine 32bit Mantisse schafft halt nur maximal 10stellige Dezimalstellen.

Man könnte hier auch den C++-Typ float verwenden (single precision), der nur 4 Bytes groß ist und eine kleinere Mantisse hat, wschl. ähnlicher zum CPC-Typ. Und dann L entsprechend anpassen.

— prodatron
... hier wäre ein Abrunden eh korrekt.

Wie definierst Du korrektes Runden? Mglst. wenig Fehler zum "echten" Wert?

— prodatron
Daher wäre es für Dich hier (in dem Beispiel jedenfalls) nicht relevant, das Runden zu beeinflussen?

Wenn ich definitiv weiß, dass abgerundet wird, dann muss ich hier den Rundungsmodus nicht beeinflussen und kann mit der Negierungsmethodik arbeiten. Leider haben moderne Rechner als Standad mWn "round-to-nearest-representable", und das kann in beide Richtungen gehen, und vielleicht macht der CPC das auch.


Mir gefällt die "Bewiesenheit" der Intervall-Arithmetik, dass die echte Zahl immer im Endintervall liegt. Damit lässt sich schon sehr viel anfangen (für den Preis von Performance-Verlust). Nimmt man einen Datentyp mit nur einer Dezimalnachkommastelle, und will die repräsentierbare Zahl 1.1 quadrieren, wird das exakte Ergebnis 1.21 zum Intervall [1.2 .. 1.3], und enthält das nicht darstellbare Ergebnis. Und dann gibt es ganze Rechenregeln für arithmetische Operationen, die schon mit Intervallen starten.

Experimenteller Code:

Hier habe ich den CPC-Typ float verwendet und nach Operationen die Mantisse des Ergebnisses um +-1 (Variable SIGADD, kann modifiziert werden) verändert.  Unter der Annahme, dass der CPC dem IEEE-Standard folgt, und ich diesen richtig interpretiere (Fehler +- 0.5 ulp bei Grundoperationen), müsste dann das korrekte Ergebnis innerhalb dieser Mantisse-veränderten Werte liegen.

Je nachdem, wofür ich die Float-Operation durchführe, nehme ich dann den kleineren oder größeren der beiden veränderten Werte.

Mein Ziel hier war es, zu sehen, ob man so einen substanziellen Zeitgewinn erreicht. Der Code ist nicht verifiziert, subnormale Float-Zahlen werden nicht berücksichtigt, mir hat es genügt, ein Ergebnisbild zu bekommen, das in etwa der Mandelbrotmenge entspricht, wenn in der Formel z^2+c der c-Wert ein komplexes Intervall ist (und einem Pixel entspricht; siehe Bild, ganz rechts).

Für ein 64x64-Bild (2. von links) brauchte der Emulator in Echtzeit nun knapp 5h anstelle von 3 Tagen  mit meinem STRING$-basierten Typ für 50x50 (ganz links). Das 256x256-Bild (3. von links) wurde nicht in Echtzeit berechnet (8h in cpcemu, fast)

Die Linien auf dem 2. und 3. Bild kommen wsch. daher, dass ich einige float-Ergebnisse als Fehler betrachte und den entsprechenden Pixel als UNCLEAR (grau) belasse. Bspw. verändere ich nur Mantissen, wenn das Ergebnis immer noch eine normalisierte Float-Zahl ist, ich also aus Gründen der Einfachheit des Codes den Exponenten nicht manipulieren muss (Zeilen 21500-21999). Fehler ist auch, wenn die Multiplikation 0 liefert, obwohl keiner der Faktoren 0 ist (underflow).

Man kann viel mit dem CPC machen, was Optimierungen angeht. Ich finde die  STRING$-Variante (Fixpunkt) aber deutlich einfacher vom Verständnis (Schulrechnen) und Programmieren als float/IEEE/ulp.
📎 _mset_cpc_ia.gif
📎 msetfm_basic.zip

Könnte man daraus berechnen wie viel schneller ein aktueller PC gegenüber dem CPC ist?

Ca 39 Millionen mal schneller.

Ich habe meine C++-Routine 1 Million mal laufen lassen für ein identisches Bild (64x64, 15 max it), damit die Orbits so ähnlich wie möglich sind. Dazu braucht diese 7:46 min (ggue 5:11 h fuer 1 Bild des CPC). Allerdings ist mein CPC-Code bei weitem nicht optimal und besteht aus sehr vielen IF-Tests, von denen mglw. manche nie wahr sein werden.

Setting:

[*]C++ Intervall-Datentyp basierend auf double precision und konsantem Aufwärtsrunden
[*]TDM-GCC-9 Compiler 64 bit
[*]AMD Ryzen 9 3900X 12-core
[*]Windows 10 64bit, Vordergrundprozess, single thread, normale Priorität

— marcm200
Ca 39 Millionen mal schneller.

Puh! Hätte eher mit so ca. 100.000 gerechnet. Aber was solls. Im Praxisbetrieb ist der CPC bei anständigem Ausbau und guten Programmen dem PC doch noch in einigem überlegen.  😉 🙂

— TFM

— marcm200
Ca 39 Millionen mal schneller.

Puh! Hätte eher mit so ca. 100.000 gerechnet. Aber was solls. Im Praxisbetrieb ist der CPC bei anständigem Ausbau und guten Programmen dem PC doch noch in einigem überlegen.  😉 🙂

Also mir fallen da nur zwei Sachen ein:
- Bootgeschwindigkeit
- Anfälligkeit für Viren/Trojaner

Will jetzt nicht vom Thema ablenken, möglicherweise Thema teilen?

Mit fällt da aus dem realen Leben viel ein:
- z.B. Brief schreiben (Ja, die Ämter wollen immer noch alles auf Papier). Wie mache ich das? CPC einschalten, Briefkopf aus der Akku-gepufferten RAM disc laden, Brief ergänzen, ausdrucken und fertig! Dazu nutze ich z.B. Protext im ROM (oder FutureTex). So schnell ist kein PC.

- Nachschlagen wichtiger Daten (die ich oft benötige). Ja, auch hier zeigen sich die Stärken von Software im ROM und Daten auf der RAM Disc (ja, natürlich Batterie-gepuffert). Alternativ zur RAM-Disk kann man auch z.B. die SD-Karte des M4 nutzen.

- es gibt noch viel viel viel mehr...  :jubelaola: