[C++] kann ich mich auf die Float-Arithmetik verlassen?

Programmiersprachen, APIs, Bibliotheken, Open Source Engines, Debugging, Quellcode Fehler und alles was mit praktischer Programmierung zu tun hat.
Antworten
Benutzeravatar
Schrompf
Moderator
Beiträge: 4838
Registriert: 25.02.2009, 23:44
Benutzertext: Lernt nur selten dazu
Echter Name: Thomas Ziegenhagen
Wohnort: Dresden
Kontaktdaten:

[C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Schrompf »

Moin,

mein Projekt soll prozedural generierte Spielwelten haben. Da das aber ziemliche Datenmassen sind, wollte ich die eigentlich nicht wie Minecraft und Konsorten komplett auf dem Server generieren und über's Netzwerk an alle Clients schieben. Stattdessen habe ich mir überlegt, dass auf jedem Client der selbe Landschaftsgenerator mit dem selben Seed und den selben Parametern läuft und ich quasi nur das Changeset dazu übertrage.

Contentgenerierung passiert komplett auf der CPU. Ich werde also zumindest nicht in die selbe Falbe wie Zudo treten :) Aber ich vermute, es gibt noch andere Fallen. Dazu meine Frage:

Kann ich mich darauf verlassen, dass identische float-Operanden und die selben Operationen auch identische Ergebnisse ergeben?

Muss ja nicht bitweise identisch sein, aber wenn z.b. ein Compiler beschließt, die Operationen ein bisschen umzustapeln und dadurch manche Operationen in einer anderen Präzision ausgeführt werden, bin ich doch auch schon am Arsch. Selbst kleinste Abweichungen könnten dann zufällig andere Pfade im Generator nehmen und komplett unterschiedliche Spielwelten ergeben. Zielhardware ist vorerst x86_64, das erfasst neben Windows/Linux/OSX ja inzwischen auch die beiden großen Konsolen. Scheiß auf Mobile.
Früher mal Dreamworlds. Früher mal Open Asset Import Library. Heutzutage nur noch so rumwursteln.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Krishty »

Du musst
  1. mit fp:precise kompilieren (sonst passiert das Umstapeln der Operationen durch den Compiler)
  2. vor jedem Beginn der Berechnung und nach jeder Unterbrechung durch Funktionsaufrufe mit globaler Wirkung deine FPU- (bzw. SSE- und AVX-)Flags wiederherstellen (ganz einfach aus dem Grund, dass z.B. ein Verschieben deines Fensters die Shell laden könnte, die wiederum alle Shell Extensions lädt, und eine davon stellt die FPU-Flags um)
  3. vorsichtig sein mit Intrinsics, da oft nicht vollständig IEEE-kompatibel und damit nicht portabel
Dann müsste das laufen wie Butter.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Schrompf
Moderator
Beiträge: 4838
Registriert: 25.02.2009, 23:44
Benutzertext: Lernt nur selten dazu
Echter Name: Thomas Ziegenhagen
Wohnort: Dresden
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Schrompf »

Danke. Ich habe beim Weitergoogeln den hier gefunden: https://randomascii.wordpress.com/2013/ ... terminism/

Und der hat ne Menge weiterer Kopfkratzer parat. Er meint aber auch, laut IEEE-Dings-Bumms könnte ich mich zumindest auf Plus, Minus, Mal, Durch, Wurzel verlassen. Damit müsste ich mich auch auf Float-To-Int-Konvertierungen verlassen können. Und damit kann ich mir irgendwo ne hinreichend große Sinus-Tabelle hinlegen und habe die letzte Operation erschlagen, die ich (vorerst) brauche.
Früher mal Dreamworlds. Früher mal Open Asset Import Library. Heutzutage nur noch so rumwursteln.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Krishty »

The result of a + b is stored in a temporary destination of unspecified precision.
Debug versus release versus levels of optimization
Macht unter x64 kein Compiler mehr. Die speichern alle temporären doubles als double und floats als floats, schon allein weil die Konvertierung Zeit kostet (was ja unter x87 genau andersrum war).

Womit du tatsächlich Probleme kriegen wirst sind sin(), cos(), usw, da hat er recht. [Er hatte sowieso recht, nur haben sich die Compiler seit dem Artikel zum Glück stark verbessert.] Wenn dir float ausreicht, würde ich die double-Versionen nutzen und zurückrunden. Oder halt deine eigene schreiben ;) aber das ist PITA mit Genauigkeitsmessungen und für double irre viel Aufwand.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Spiele Programmierer
Establishment
Beiträge: 426
Registriert: 23.01.2013, 15:55

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Spiele Programmierer »

Hast du schonmal von dem Spiel Factorio gehört? Deren gesamte Multiplayer-Implementierung funktioniert, indem die Simulation absolut deterministisch abläuft und nur die Aktionen der Spieler übertragen werden. Also es ist definitiv möglich.

Soweit ich weiß, hat das über die Zeit allerdings zu einer sehr großen Menge an Problemen geführt. Zum Debugging vergleichen bei Factorio Server und Client Hash-Werte der Spielewelt. Ich kann gerade keinen technischen Artikel mehr dazu finden, aber wenn ich mich recht erinnere, waren sie auch gezwungen sie sehr viele Funktionen (zB. Sinus und Kosinus) selbst zu Implementieren.

Des Weiteren würde ich dir empfehlen, den Compiler sogar auf /fp:strict zu stellen.
/fp:precise aktiviert einige nicht IEEE Standard konforme Optimierungen:
https://msdn.microsoft.com/en-us/library/e7s85ffb.aspx hat geschrieben:The following floating-point behavior is enabled with /fp:precise:
Contractions—that is, using a composite operation that has just one rounding at the end to replace multiple operations.
Auf x86 werden scheinbar ggf. die alten x87 Befehle ohne Zwischenrundung(!) in Ausdrücken benützt.
Auf x64 kann wohl zB. a * b + c zu Fused Multiply Add zusammengefasst werden, falls die CPU das unterstützt und es im Compiler aktiviert wurde.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Krishty »

mad() kann auch explizit (als Intrinsic) problematisch werden, weil z.B. 3D-Skalarprodukte normalisierter Vektoren größer als 1 ausfallen können …
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Schrompf
Moderator
Beiträge: 4838
Registriert: 25.02.2009, 23:44
Benutzertext: Lernt nur selten dazu
Echter Name: Thomas Ziegenhagen
Wohnort: Dresden
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Schrompf »

Uiuiui. Danke nochmal für euren Input. Auch das Internet meint dazu: ja, es ist möglich, deterministische Algorithmen zu schreiben, aber es gibt zwei Millionen Stolperstellen. Mein Plan lautet jetzt also, am Ende die globalen Strukturen doch in Integer-Mathematik zu produzieren und beim Verfeinern dann in einem lokalen Koordinatensystem mit floats zu arbeiten. Mal gucken, ob das praktikabel ist. Irgendne Art von Verifikation wollte ich aber auch einbauen - notfalls pumpe ich zur Überprüfung in bester Minecraft-Manier das komplette Segment über's Netzwerk. Dank der inzwischen nicht mal mehr neuen Kompression passen 2048^3 Voxel in vielleicht 10MB. Das machen heutige VServer stressfrei, ist ja kein Dauerbetrieb bei den Spielern.
Früher mal Dreamworlds. Früher mal Open Asset Import Library. Heutzutage nur noch so rumwursteln.
DerAlbi
Establishment
Beiträge: 269
Registriert: 20.05.2011, 05:37

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von DerAlbi »

Fixkomma FTW
Halte ich wirklich für eine sinnvolle Alternative. Gerade, wenn es um Welten geht, wo der Wertebereich der Koordinaten deterministisch auf eine Größeneinheit normiert ist.
Float ist toll, wenn man Werte gleichverteilt auf einer logarithmischen Achse hat..
Aber joar.. man muss dafür halt ne Lib schreiben, damits bequem nutzbar ist.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Krishty »

DerAlbi hat geschrieben:Fixkomma FTW
Halte ich wirklich für eine sinnvolle Alternative. Gerade, wenn es um Welten geht, wo der Wertebereich der Koordinaten deterministisch auf eine Größeneinheit normiert ist.
Float ist toll, wenn man Werte gleichverteilt auf einer logarithmischen Achse hat..
Aber joar.. man muss dafür halt ne Lib schreiben, damits bequem nutzbar ist.
Dass Weltkoordinaten fast nie logarithmisch sind, stimmt. Schrompf wird aber durch nichts davon abgehalten, die Welt in Abschnitte von 16³ Voxel zu kacheln und auf denen jeweils eine float-Simulation laufen zu lassen. float ist halt für die meisten Anwendungsfälle einfacher nutzbar (sin() und Konsorten sind fertig in jeder Standardbibliothek), und die Gleitkommaleistung der CPUs ist deutlich höher als die Integerleistung.

Ich setze bei mir auch Festkommazahlen für Kollisionen ein, und die funktionieren quasi perfekt. Dummerweise hat aber die Entwicklung vier Monate gebraucht statt einer Woche.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Schrompf
Moderator
Beiträge: 4838
Registriert: 25.02.2009, 23:44
Benutzertext: Lernt nur selten dazu
Echter Name: Thomas Ziegenhagen
Wohnort: Dresden
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Schrompf »

Jupp, genau dieser Zusatzaufwand lässt mich gruseln. Ich hab jetzt ne kleine Tochter, ich komme pro Tag bestenfalls ne Stunde zum Coden. Da werde ich mir keine Festkomma-Arithmetik mit dazugehörigen Unittest hinzaubern. Es soll funktionieren, und zwar möglichst aufwandarm. Ne kleine Tabelle mit bitweise definierten Floats als Sinus-Ersatz ist da ne deutlich entspanntere Lösung als alles mit dem Rumgeshifte für Festkommazahlen neu zu implementieren.
Früher mal Dreamworlds. Früher mal Open Asset Import Library. Heutzutage nur noch so rumwursteln.
Benutzeravatar
B.G.Michi
Establishment
Beiträge: 163
Registriert: 07.03.2006, 20:38
Alter Benutzername: B.G.Michi
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von B.G.Michi »

Na du brauchst das Rad ja nicht neu erfinden. Wenn du wirklich Festkommaarithmetik verwenden möchtest, wette ich du findest mit Google und Stackoverflow in 10 min eine (halbwegs) brauchbare Library die alle benötigten Funktionen liefert und mit etwas Glück mit passender Lizenz. Festkommalibraries gibts bestimmt wie Sand am Meer...
DerAlbi
Establishment
Beiträge: 269
Registriert: 20.05.2011, 05:37

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von DerAlbi »

Ich sitz grade dran ^_^
Ist mein zweiter Versuch sowas zu implementieren und mitlerweile siehts recht gut aus. Die Begrenzungen, die ich das letze mal hatte, fliegen dank einer art lazy-evaluation aus.
Insgesamt sieht der Funktionsumfang bzw die Operator-Optimierung im Moment so aus:

Code: Alles auswählen

//Assign korrektheit
Fixpoint<24> = Fixpoint<10> * Fixpoint<20>; // wird zu
Fixpoint<24> = Fixpoint<30> >> 6; 

Fixpoint<16> = Fixpoint<12> + Fixpoint<20>; // wird zu
Fixpoint<16> = Fixpoint<12> << 4 + Fixpoint<20> >> 4; 

Fixpoint<16> = Fixpoint<12> + Fixpoint<12>;  //wird zu
Fixpoint<16> = Fixpoint<12> << 4; 

Fixpoint<16> = Fixpoint<8> + Fixpoint<12>;  //wird zu
Fixpoint<16> = (Fixpoint<8> << 4 + Fixpoint<12>) << 4; 

Fixpoint<12> = Fixpoint<16> + Fixpoint<20>;  //wird zu
Fixpoint<12> = Fixpoint<16> + Fixpoint<20> >> 4) >> 4; 


//Left-Hand-Side Dominanz
Fixpoint<22> = Fixpoint<10> *  ( Fixpoint<20_lhs> + Fixpoint<15> << 5 ) //wird zu 
Fixpoint<22> = Fixpoint<10> *  Fixpoint<20_lhs> //wird zu 
Fixpoint<22> = Fixpoint<30> >> 8

Fixpoint<22> = Fixpoint<10> *  ( Fixpoint<15_lhs> + Fixpoint<20> >> 5) //wird zu 
Fixpoint<22> = Fixpoint<10> *  Fixpoint<15_lhs>  // Wird zu
Fixpoint<22> = Fixpoint<30> >> 8


Multiply-Add - Korrektheit
Fixpoint<24> = Fixpoint<10>*Fixpoint<20> + Fixpoint<20>    //wird zu
Fixpoint<24> = Fixpoint<30> >> 6 + Fixpoint<20> << 4 

Fixpoint<24> = Fixpoint<10>*Fixpoint<20> + Fixpoint<15>*Fixpoint<15>   //wird zu
Fixpoint<24> = ( Fixpoint<30> + Fixpoint<30> ) >> 6  

Fixpoint<24> = Fixpoint<10>*Fixpoint<20> + Fixpoint<15>*Fixpoint<15>   + Fixpoint<10>*Fixpoint<20>  //wird zu
Fixpoint<24> = ( Fixpoint<30> + Fixpoint<30> + Fixpoint<30> ) >> 6  

Fixpoint<24> = Fixpoint<10>*Fixpoint<10> + Fixpoint<15>*Fixpoint<15>   + Fixpoint<10>*Fixpoint<20>  //wird zu
Fixpoint<24> = ( Fixpoint<20_lhs> + Fixpoint<30> >> 10 ) + Fixpoint<30>    //Left-Hand-Side-Dominanz
Fixpoint<24> = Fixpoint<20_lhs> << 4 + Fixpoint<30> >> 6

Fixpoint<24> = Fixpoint<20> * Fixpoint<10> + Fixpoint<10>*Fixpoint<10>   + Fixpoint<10>*Fixpoint<20>  //wird zu
Fixpoint<24> = ( Fixpoint<30_lhs> + Fixpoint<20> << 10) + Fixpoint<30>    //Left-Hand-Side-Dominanz
Fixpoint<24> = ( Fixpoint<30_lhs> + Fixpoint<30>  ) >> 6


//IntegerKonstanten werden direkt umgesetzt

Fixpoint<24> = Fixpoint<10> + int;  //oder
Fixpoint<24> = int + Fixpoint<10>;  // wird zu
Fixpoint<24> = Fixpoint<10> + Fixpoint<0>(int); 
Fixpoint<24> = Fixpoint<10> << 14 + Fixpoint<0>(int) << 24; 

Fixpoint<24> = Fixpoint<10> * int;  //oder
Fixpoint<24> = int * Fixpoint<10>;  // wird zu
Fixpoint<24> = Fixpoint<0>(int) * Fixpoint<10>; 
Fixpoint<24> = Fixpoint<10> << 14;

Fixpoint<24> = Fixpoint<20>*int + Fixpoint<10>*int + int*Fixpoint<20>; //wird zu
Fixpoint<24> = (Fixpoint<20_lhs> + Fixpoint<10> << 10 ) + Fixpoint<20>; 
Fixpoint<24> = (Fixpoint<20_lhs>  + Fixpoint<20> ) << 4;


//Float-Konstanten werden zu Fixpoint<> konvertiert und verarbeitet
Fixpoint<24> = Fixpoint<20>*float + Fixpoint<10>*double + float*Fixpoint<20>; //wird zu
Fixpoint<24> = Fixpoint<20>*Fixpoint<20>(float) + Fixpoint<10>*Fixpoint<10>(double) + Fixpoint<20>*Fixpoint<20>(float); 
Fixpoint<24> = (Fixpoint<40_lhs> + Fixpoint<20> << 20 ) + Fixpoint<40>;
Fixpoint<24> = (Fixpoint<40_lhs> + Fixpoint<40>) >> 16
Soweit sieht das ganz gut aus. Die Division wird noch hart. Natürlich ist der Basistyp frei wählbar (char, short, int) +/- signed/unsigned und die Multiplikationen werden in doppelt so großen Datentypen evaluiert. Falls man Typen mixt, gilt die gleiche Dominanz-Regel, wie bei der Präzision.
Benutzeravatar
Jonathan
Establishment
Beiträge: 2352
Registriert: 04.08.2004, 20:06
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Jonathan »

Schrompf hat geschrieben:Jupp, genau dieser Zusatzaufwand lässt mich gruseln.
Wirklich Angst machen würde mir bei der Geschichte, dass sich irgendwann einmal irgendetwas am Compiler ändert und auf einen Schlag das Spiel kaputt ist. Schlimmer noch, vielleicht ist es nur sehr subtil kaputt und man merkt es erst überhaupt nicht. Und dann exakt herauszufinden, wie damals was genau definiert war, und was jetzt genau irgendwie anders berechnet sein könnte, wäre ein ziemlicher Debugging-Alptraum.
Ich habe genügend Spiele die auf modernen Betriebssystemen nicht mehr laufen. Da reizt mich irgendwie die Idee, zumindest meinen Kram in 10 Jahren neu kompilieren zu können und er läuft dann immer noch.
Lieber dumm fragen, als dumm bleiben!
https://jonathank.de/games/
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Krishty »

Jonathan hat geschrieben:Ich habe genügend Spiele die auf modernen Betriebssystemen nicht mehr laufen.
Das liegt aber nicht am Compiler, sondern am Abschaffen der FPU. Das war auch lange absehbar.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Schrompf
Moderator
Beiträge: 4838
Registriert: 25.02.2009, 23:44
Benutzertext: Lernt nur selten dazu
Echter Name: Thomas Ziegenhagen
Wohnort: Dresden
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Schrompf »

Habe jetzt doch eine FixPoint-Klasse geschrieben. Template-Konfigurierbarkeit und den ganzen Quatsch habe ich mir erspart, das Ding ist erstmal hardcoded auf 40.24 in einem uint64_t. Dafür hat die Klasse aber funktionierende und überlaufvermeidende Multiplikation und Division. Alles, was ich an fertigem Code da draußen gefunden habe, hat die Multiplikation z.b. komplett ausgelassen oder ganz plump mit (A * B) >> SHIFT implementiert, was Dir bei nem uint64_t natürlich sofort überläuft. Von der Division brauchen wir gar nicht erst zu reden.

Wer mit Gruselcode (weil in krudem Deutsch) klar kommt, darf gern mal einen Blick drauf werfen. Die UnitTests sind im zweiten Block weiter unten. Ich bin dankbar für nahezu jedes fachliche Gemecker.

Code: Alles auswählen

/// @file FestPunkt.h
/// Repräsentiert eine Festkomma-Zahl mit passender Arithmetik

#pragma once

struct FestPunkt
{
  uint64_t w; // Festkommawert mit AnzKommaBits Nachkomma. Ohne Vorzeichen, weil damit das Bitgewusel einfacher zu kontrollieren ist.

  static constexpr size_t AnzGesamtBits = 8 * sizeof( decltype( w));
  static constexpr size_t AnzKommaBits = 24;
  static constexpr size_t AnzVorKommaBits = AnzGesamtBits - AnzKommaBits;
  static constexpr uint64_t NachKommaMaske = BitHelfer::MachEins<uint64_t>(AnzKommaBits);
  static constexpr float FloatZuFP = float( 1u << AnzKommaBits);
  static constexpr double DoubleZuFP = double( 1u << AnzKommaBits);
  static constexpr float FPZuFloat = 1.0f / FloatZuFP;
  static constexpr double FPZuDouble = 1.0 / DoubleZuFP;

  /// Standard-Konstruktor mit 0
  constexpr FestPunkt() : w( 0) { }
  /// Konstruktor mit Ganzzahl
  /*explicit*/ constexpr FestPunkt( int64_t i ) : w( i ) { }
  /// Konstruktor mit Ganzzahl
  /*explicit*/ constexpr FestPunkt( int i) : w( i) { }
  /// Konstruktor aus Float
  /*explicit*/ constexpr FestPunkt(float f) : w( int64_t( f * FloatZuFP)) { }
  /// Konstruktor aus Double
  /*explicit*/ constexpr FestPunkt(double d) : w( int64_t( d * DoubleZuFP)) { }
  /// Benamter Konstruktor zur Instanziierung aus Festkomma-Bits
  static FestPunkt FestPunktAusBits( uint64_t w) { FestPunkt fp; fp.w = w; return fp; }

  /// Kopieren
  FestPunkt(const FestPunkt&) = default;
  FestPunkt& operator = (const FestPunkt&) = default;

  /// Konvertierung zu Float
  constexpr float GetFloat() const { return float( int64_t( w)) * FPZuFloat; }
  /// Konvertierung zu Double
  constexpr double GetDouble() const { return double( int64_t( w)) * FPZuDouble; }

  /*constexpr*/ FestPunkt& operator += (FestPunkt fp) { w = (*this + fp).w; return *this; }
  /*constexpr*/ FestPunkt& operator -= (FestPunkt fp) { w = (*this - fp).w; return *this; }
  /*constexpr*/ FestPunkt& operator *= (FestPunkt fp) { w = (*this * fp).w; return *this; }
  /*constexpr*/ FestPunkt& operator /= (FestPunkt fp) { w = (*this / fp).w; return *this; }

  friend FestPunkt operator + (FestPunkt a, FestPunkt b) { return FestPunktAusBits( a.w + b.w); }
  friend FestPunkt operator - (FestPunkt a, FestPunkt b) { return FestPunktAusBits( a.w - b.w); }
  friend FestPunkt operator * (FestPunkt a, FestPunkt b) 
  {
    static_assert(AnzKommaBits < 32, "Mit 32 Nachkommabits geht der Übertrag verloren. Bessere Implementation nötig.");
    // A * B
    // = ((ah + al) * (bh + bl)) >> AnzKommaBits;
    // = ((ai*2^AKB + al) * (bi*2^AKB + bl)) >> AKB;
    // = ((ai*bi*2^2*AKB) + (ai*bl*2^AKB) + (bi*al*2^AKB) + (al*bl)) >> AKB;
    // = ((ai*bi*2^AKB) + (ai*bl) + (bi*al))   +   ((al*bl) >> AKB);
    uint64_t fa = std::abs( int64_t( a.w)), fb = std::abs( int64_t( b.w));
    bool ergebnisIstNegativ = ((a.w >> 63) ^ (b.w >> 63)) != 0;

    uint64_t al = fa & NachKommaMaske, bl = fb & NachKommaMaske;
    uint64_t nachkomma = (al * bl) >> AnzKommaBits;
    uint64_t ai = fa >> AnzKommaBits, bi = fb >> AnzKommaBits;
    uint64_t vorkomma = ((ai * bi) << AnzKommaBits) + (ai*bl) + (bi*al);
    uint64_t abserg = vorkomma + nachkomma;
    return FestPunktAusBits( uint64_t( ergebnisIstNegativ ? -int64_t( abserg) : abserg));
  }

  friend FestPunkt operator / (FestPunkt a, FestPunkt b) 
  {
    // A / B wie in der Schule:
    // A / B = x
    // A - B*x = rest
    // rest << AKB / B
    // = ah / B + al / B
    uint64_t dividend = std::abs( int64_t( a.w)), divisor = std::abs( int64_t( b.w));
    bool ergebnisIstNegativ = ((a.w >> 63) ^ (b.w >> 63)) != 0;

    // jetzt in Acht-Bit-Schritten jeweils Dividend rausschieben, dividieren, Ergebnis als 8 Nachkommabits eintragen
    static constexpr size_t NumBitsProSchritt = 8;
    static constexpr size_t NumSchritte = AnzKommaBits / NumBitsProSchritt;
    static_assert((AnzKommaBits % NumBitsProSchritt) == 0, "Schrittweise Division muss zu Anzahl Nachkommabits passen");
    uint64_t erg = dividend / divisor, div = erg;

    for( size_t i = 0; i < NumSchritte; ++i )
    {
      erg <<= NumBitsProSchritt;
      dividend = (dividend - div * divisor) << NumBitsProSchritt;
      div = dividend / divisor;
      erg += div;
    }
    return FestPunktAusBits( uint64_t( ergebnisIstNegativ ? -int64_t( erg) : erg));
  }

  friend bool operator == (FestPunkt a, FestPunkt b) { return a.w == b.w; }
  friend bool operator != ( FestPunkt a, FestPunkt b ) { return a.w != b.w; }
  friend bool operator <= ( FestPunkt a, FestPunkt b ) { return int64_t( a.w) <= int64_t( b.w); }
  friend bool operator >= ( FestPunkt a, FestPunkt b ) { return int64_t( a.w) >= int64_t( b.w); }
  friend bool operator < ( FestPunkt a, FestPunkt b ) { return int64_t( a.w) < int64_t( b.w); }
  friend bool operator > ( FestPunkt a, FestPunkt b ) { return int64_t( a.w) > int64_t( b.w); }
};
Ich würde das alles gerne constexpr-fähig machen, aber Visual Studio ist da ein bisschen hinterher :-( Nuja, die UnitTests:

Code: Alles auswählen

/// Testet die FestPunkt-Arithmetik
#include "UT_PCH.h"

TEST( Traum_FestPunkt, Konvertierung)
{
  const double zahlen[] = { 0.0, 1.0, 0.25, 1024.0, 3.1415926, -1.0, -0.00375, -25385024.37054, 826742735.26401 };
  for( auto z : zahlen )
  {
    FestPunkt fp( z);
    double maxdiff = pow( 2.0, std::max( -24.0, std::log2( z) - 48.0));
    EXPECT_NEAR( z, fp.GetDouble(), maxdiff);
  }
}

TEST( Traum_FestPunkt, PlusMinus)
{
  const FestPunkt zahlen[] = { 0.0, 1.0, 0.25, 1024.0, 3.1415926, -1.0, -0.00375, -25385024.37054, 826742735.26401 };
  for( auto z1 : zahlen )
  {
    for( auto z2 : zahlen )
    {
      auto ergfp = z1 + z2;
      double ergd = z1.GetDouble() + z2.GetDouble();
      double maxdiff = pow( 2.0, std::max( -24.0, std::log2( std::abs( ergd)) - 48.0));
      EXPECT_NEAR( ergfp.GetDouble(), ergd, maxdiff );
    }
  }

  for( auto z1 : zahlen )
  {
    for( auto z2 : zahlen )
    {
      auto ergfp = z1 - z2;
      double ergd = z1.GetDouble() - z2.GetDouble();
      double maxdiff = pow( 2.0, std::max( -24.0, std::log2( std::max( std::abs( z1.GetDouble()), std::abs( z2.GetDouble())) - 48.0)));
      EXPECT_NEAR( ergfp.GetDouble(), ergd, maxdiff );
    }
  }
}

TEST( Traum_FestPunkt, Multiplikation )
{
  const FestPunkt zahlen[] = { 0.0, 1.0, 0.25, 1024.0, 3.1415926, -1.0, -0.00375, -25385024.37054, 936572.58274, -372058.35693, 826742735.26401 };
  for( auto z1 : zahlen )
  {
    for( auto z2 : zahlen )
    {
      // unsere Festkomma hat maximal 40 Vorkommabits. Multiplikationen, die jenseits davon rauskommen, brauchen wir gar nicht erst prüfen
      double ergd = z1.GetDouble() * z2.GetDouble();
      if( std::log2( std::abs( ergd)) >= 39.0 )
        continue;
      auto ergfp = z1 * z2;
      double maxdiff = pow( 2.0, std::max( -24.0, std::log2( std::max( std::abs( z1.GetDouble()), std::abs( z2.GetDouble())) - 48.0)));
      EXPECT_NEAR( ergfp.GetDouble(), ergd, maxdiff);
    }
  }
}

TEST( Traum_FestPunkt, Division )
{
  const FestPunkt zahlen[] = { 1.0, 0.25, 1024.0, 3.1415926, -1.0, -0.00375, -25385024.37054, 936572.58274, -372058.35693, 826742735.26401, 1'834'735'462.49547 };
  for( auto z1 : zahlen )
  {
    for( auto z2 : zahlen )
    {
      // Division geht bis 40bit Divident bzw. 32bit Vorkomma Divisor
      double ergd = z1.GetDouble() / z2.GetDouble();
      auto ergfp = z1 / z2;
      double maxdiff = pow( 2.0, std::max( -24.0, std::log2( std::abs( ergd)) - 48.0));
      EXPECT_NEAR( ergfp.GetDouble(), ergd, maxdiff );
    }
  }
}

TEST( Traum_FestPunkt, Vergleiche)
{
  const FestPunkt zahlen[] = { 1.0, 0.25, 1024.0, 3.1415926, -1.0, -0.00375, -25385024.37054, 936572.58274, -372058.35693, 826742735.26401, 1'834'735'462.49547, 453'305'284'540.204274, -852'289'572'499.38583 };
  for( auto z1 : zahlen )
  {
    for( auto z2 : zahlen )
    {
      EXPECT_EQ( (z1 == z2), (z1.GetDouble() == z2.GetDouble()) );
      EXPECT_EQ( (z1 != z2), (z1.GetDouble() != z2.GetDouble()) );
      EXPECT_EQ( (z1 <= z2), (z1.GetDouble() <= z2.GetDouble()) );
      EXPECT_EQ( (z1 >= z2), (z1.GetDouble() >= z2.GetDouble()) );
      EXPECT_EQ( (z1 <  z2), (z1.GetDouble() <  z2.GetDouble()) );
      EXPECT_EQ( (z1 >  z2), (z1.GetDouble() >  z2.GetDouble()) );
    }
  }
}
Früher mal Dreamworlds. Früher mal Open Asset Import Library. Heutzutage nur noch so rumwursteln.
DerAlbi
Establishment
Beiträge: 269
Registriert: 20.05.2011, 05:37

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von DerAlbi »

Huhu :-)
darf ich fragen, wieso du auf so perverse megatypgröße wie 40:24 gehst?
ist das
a) um alles abzudecken
b)oder würdest du kleinere Fixpoints nutzen können, wenn du dir um die architektur mehr gedanken machen würdest?
c) willst du dir überhaupt um die Fixpoints gedanken machen müssen?

(interessiert mich wirklich: weil ich ja nun so eine fixpoint klasse habe die alles erschlagen kann... mich interessieren andere usecasses, um zu gucken, ob die auch erschlagen werden)
Es funktionieren bei mir leider nur 32bit-Fixpoints komplett. 64bit funktioniert nicht, solang es keinen offiziellen 128-bit typ gibt. :-(
Benutzeravatar
Schrompf
Moderator
Beiträge: 4838
Registriert: 25.02.2009, 23:44
Benutzertext: Lernt nur selten dazu
Echter Name: Thomas Ziegenhagen
Wohnort: Dresden
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Schrompf »

Ich brauche das halt, um alle Skalen meines VoxelSurvival-Landschaftsgenerators abzudecken. Der muss zum Einen die Grobstrukturen in Kilometergröße entwerfen und die dann aber auch bis in den Zentimeter-Bereich mit genügend hoher Sub-Voxel-Genauigkeit für z.B. leicht geneigte Strecken oder sowas zu haben. 16.16 war mir jedenfalls zu wenig Nachkomma und gleichzeitig viel zu wenig Reichweite für die großen Strukturen. Also ist die nächste Stufe halt 64bit. Und da steht in der Multiplikation ein Kommentar, wonach wegen der leicht vereinfachten Implementation 32bit Nachkomma zu groß wären. Also 24bit Nachkomma, Rest für Vorkomma.
Früher mal Dreamworlds. Früher mal Open Asset Import Library. Heutzutage nur noch so rumwursteln.
DerAlbi
Establishment
Beiträge: 269
Registriert: 20.05.2011, 05:37

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von DerAlbi »

Rein Konzeptuell: wie hast du das vorher mit floats gemacht? Weil das, was du beschreibst, klingt so, als könnte man es mit float, wenn man in Absolutkoordinaten rechnet, nicht erreichen.
Hast du die Blöcke lokal in einem Zahlenraum behalten und dann an die richtige Stelle transformiert? (wo es dann gut aussah, weil die schlechte auflösung weit weg ist..)
Benutzeravatar
Schrompf
Moderator
Beiträge: 4838
Registriert: 25.02.2009, 23:44
Benutzertext: Lernt nur selten dazu
Echter Name: Thomas Ziegenhagen
Wohnort: Dresden
Kontaktdaten:

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Schrompf »

Nachtrag: Sinus. Bzw. Cosinus, weil symmetrisch. Begonnen habe ich mit ner Taylor-Reihe, dann aber wie in spannenden Ausführungen empfohlen auf ein Polynom verallgemeinert und für die vier Stützstellen 0, 0.25π, 0.33π und 0.5π die Konstanten mit Wolfram Alpha lösen lassen. Die Ergebnis-Implementation ist auf 0,4‰ genau, also schon noch ein Stück entfernt von den Möglichkeiten eines 40.24-Festkomma-Formats, aber ausreichend. Und vor allem für meine Zwecke deterministisch.

Code:

Code: Alles auswählen

namespace FestPunktImpl
{
  constexpr FestPunkt EinsDurchZweiPi = FestPunkt::AusBits( 2670177); // (1.0f / 2*PI) * 2^24 = 2670176,85772

   /// Cosinus von Festkommawert in Kreisanteilen, also 0..1 anstatt 0..2PI
  inline constexpr FestPunkt cosinus( FestPunkt w )
  {
    static_assert(FestPunkt::AnzKommaBits == 24, "Hardcoded for 24bit Nachkomma");

    // Cosinus ist symmetrisch, also können wir einfach absolut arbeiten
    // HACK: (thom) std::abs() und Konsorten sind nicht constexpr, also müssen wir das abs() zu Fuß ausformulieren
    int64_t x = int64_t ( w.w) < 0 ? -int64_t( w.w) : int64_t( w.w);

    // Quadrant bestimmen als 0..3: 1 und 2 sind negativ, 1 und 3 sind außerdem rückwärts
    int64_t quadrant = (x >> 22) & 3;
    // Rest begrenzen auf 0..0.25
    x &= 0x00000000003fffff;
    x = (quadrant & 1) == 1 ? (0x400000 - x) : (x);

    // Hat außerdem den Vorteil, dass bei 24 Nachkommabits jetzt alle Zahlen in int32_t passen, also
    // können wir ohne große Magie zwei Zahlen in einem int64 multiplizieren
    constexpr int64_t ka = -331142797; // 2^24 * -19.7376487742810 == -331142796,818247581696
    constexpr int64_t kb = 1086454529; // 2^24 * 64.7577362626128 == 1086454528,9488876699648
    constexpr int64_t kc = -1330193214; // 2^23 * -79.2856939858653 == -1330193213,7107630850048
    int64_t qx = (x * x) >> 24;
    int64_t erg = 0x1000000 + ((ka * qx) >> 24);
    int64_t qqx = (qx * qx) >> 24;
    erg += (kb * qqx) >> 24;
    int64_t qqqx = (qx * qqx) >> 24;
    erg += (kc * qqqx) >> 24;

    // Quadranten anwenden - Rückwärts haben wir oben schon behandelt, müssen nur noch negieren
    erg = ((0b0110 >> quadrant) & 1) == 1 ? -erg : erg;
    return FestPunkt::AusBits( erg);
  }
}

inline constexpr FestPunkt sin(FestPunkt radian)
{
  // Normalisieren auf Einheitskreis, also 1 == 360°.
  auto w = radian * FestPunktImpl::EinsDurchZweiPi;
  return FestPunktImpl::cosinus( w - FestPunkt::AusBits( 0x400000));
}

inline constexpr FestPunkt cos(FestPunkt radian)
{
  // Normalisieren auf Einheitskreis, also 1 == 360°.
  auto w = radian * FestPunktImpl::EinsDurchZweiPi;
  return FestPunktImpl::cosinus( w);
}

/*
// Taylor um 0
f(x) = f0(0) + f1(0)*x/1 + f2(0)*x²/2 + f3(0)*x³/6 + f4(0)*x^4/24 + f5(0)*x^5/120 + f6(0)*x^6/720
f0(x) = f4(x) = cos(x) = 1
f1(x) = f5(x) = -sin(x) = 0
f2(x) = f6(x) = -cos(x) = -1
f3(x) = f7(x) = sin(x) = 0

-->
cos(x) = 1 + -1*x^2/2 + 1*x^4/24 + -1*x^6/720

angepasste Kurve mit drei Unbekannten
cos(x) = 1 + a*x^2 + b*x^4 + c*x^6

// anpassen an 3 Stützpunkte bei x (0..0.25) mit (0.5*0.25 -> 1/sqrt(2)), (2/3*0.25 -> 0.5) und (1*0.25 -> 0)
1/sqrt(2) = 1 + a*(1/8)^2  + b*(1/8)^4 + c*(1/8)^6
0.5       = 1 + a*(2/12)^2 + b*(2/12)^4 + c*(2/12)^6
0         = 1 + a*(1/4)^2  + b*(1/4)^4  + c*(1/4)^6

// und Wolfram Alpha antwortet:
a≈-19.7376487742810 and b≈64.7577362626128 and c≈-79.2856939858653
// die finale Formel lautet also
y = 1 - 19.7376487742810 * x^2 + 64.7577362626128 * x^4 - 79.2856939858653 * x^6

*/
Hab mir auch wieder einen Strauß UnitTests dazu geschrieben, die behaupten, dass der Co-Sinus im Rahmen der Genauigkeit passt, und das er relativ fix ist... ein paar Millionen davon pro Sekunde gehen schon. Habe aber nicht genau werkgebankt, in der Messschleife war auch noch der sin() der Standardlib drin.
Früher mal Dreamworlds. Früher mal Open Asset Import Library. Heutzutage nur noch so rumwursteln.
Spiele Programmierer
Establishment
Beiträge: 426
Registriert: 23.01.2013, 15:55

Re: [C++] kann ich mich auf die Float-Arithmetik verlassen?

Beitrag von Spiele Programmierer »

Es gibt übrigens auch die Möglichkeit die optimalen Koeffizienten für so ein Approximationspolynom eines gegebenen Grads automatisch so zu bestimmen, dass der maximale relative bzw. absolute Fehler minimal ist.

Ein kleines Program das ich in der Vergangenheit schon mal genutzt habe ist folgendes (die offizielle Seite scheint leider nicht mehr zu existieren):
Erklärung: https://web.archive.org/web/20160810075 ... ximations#
Doku: https://web.archive.org/web/20160827002 ... aths/remez
Download: https://web.archive.org/web/20161105141 ... s/lolremez
Antworten