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

[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. :-(

[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.

[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..)

[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.

[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