C++Guns – RoboBlog blogging the bot

Was macht -ffast-math eigentlich und was bringt es für Vorteile?

Aus der Dokumentation:

This option is not turned on by any -O option besides -Ofast since it can result in incorrect output for programs that depend on an exact implementation of IEEE or ISO rules/specifications for math functions. It may, however, yield faster code for programs that do not require the guarantees of these specifications.

-ffast-math
Sets the options -fno-math-errno, -funsafe-math-optimizations, -ffinite-math-only, -fno-rounding-math, -fno-signaling-nans, -fcx-limited-range and -fexcess-precision=fast.

This option causes the preprocessor macro __FAST_MATH__ to be defined.

-fno-math-errno
Do not set errno after calling math functions that are executed with a single instruction, e.g., sqrt. A program that relies on IEEE exceptions for math error handling may want to use this flag for speed while maintaining IEEE arithmetic compatibility.

The default is -fmath-errno.

-funsafe-math-optimizations

Allow optimizations for floating-point arithmetic that (a) assume that arguments and results are valid and (b) may violate IEEE or ANSI standards. When used at link time, it may include libraries or startup files that change the default FPU control word or other similar optimizations.

Enables -fno-signed-zeros, -fno-trapping-math, -fassociative-math and -freciprocal-math.

-fno-signed-zeros

Allow optimizations for floating-point arithmetic that ignore the signedness of zero. IEEE arithmetic specifies the behavior of distinct +0.0 and -0.0 values, which then prohibits simplification of expressions such as x+0.0 or 0.0*x (even with -ffinite-math-only). This option implies that the sign of a zero result isn’t significant.

-fno-trapping-math

Compile code assuming that floating-point operations cannot generate user-visible traps. These traps include division by zero, overflow, underflow, inexact result and invalid operation. This option requires that -fno-signaling-nans be in effect. Setting this option may allow faster code if one relies on “non-stop” IEEE arithmetic, for example.

-fassociative-math

Allow re-association of operands in series of floating-point operations. This violates the ISO C and C++ language standard by possibly changing computation result.
NOTE: re-ordering may change the sign of zero as well as ignore NaNs and inhibit or create underflow or overflow (and thus cannot be used on code that relies on rounding behavior like (x + 2**52) - 2**52. May also reorder floating-point comparisons and thus may not be used when ordered comparisons are required. This option requires that both -fno-signed-zeros and -fno-trapping-math be in effect. Moreover, it doesn’t make much sense with -frounding-math.

-freciprocal-math

Allow the reciprocal of a value to be used instead of dividing by the value if this enables optimizations. For example x / y can be replaced with x * (1/y), which is useful if (1/y) is subject to common subexpression elimination. Note that this loses precision and increases the number of flops operating on the value.

-ffinite-math-only

Allow optimizations for floating-point arithmetic that assume that arguments and results are not NaNs or +-Infs.

-frounding-math

Disable transformations and optimizations that assume default floating-point rounding behavior. This is round-to-zero for all floating point to integer conversions, and round-to-nearest for all other arithmetic truncations. This option should be specified for programs that change the FP rounding mode dynamically, or that may be executed with a non-default rounding mode. This option disables constant folding of floating-point expressions at compile time (which may be affected by rounding mode) and arithmetic transformations that are unsafe in the presence of sign-dependent rounding modes.

-f-signaling-nans

Compile code assuming that IEEE signaling NaNs may generate user-visible traps during floating-point operations. Setting this option disables optimizations that may change the number of exceptions visible with signaling NaNs. This option implies -ftrapping-math.

This option causes the preprocessor macro __SUPPORT_SNAN__ to be defined.

-fcx-limited-range

When enabled, this option states that a range reduction step is not needed when performing complex division. Also, there is no checking whether the result of a complex multiplication or division is NaN + I*NaN, with an attempt to rescue the situation in that case. The default is -fno-cx-limited-range, but is enabled by -ffast-math.

This option controls the default setting of the ISO C99 CX_LIMITED_RANGE pragma. Nevertheless, the option applies to all languages.

-fexcess-precision=style

This option allows further control over excess precision on machines where floating-point operations occur in a format with more precision or range than the IEEE standard and interchange floating-point types. By default, -fexcess-precision=fast is in effect; this means that operations may be carried out in a wider precision than the types specified in the source if that would result in faster code, and it is unpredictable when rounding to the types specified in the source code takes place. When compiling C, if -fexcess-precision=standard is specified then excess precision follows the rules specified in ISO C99; in particular, both casts and assignments cause values to be rounded to their semantic types (whereas -ffloat-store only affects assignments). This option is enabled by default for C if a strict conformance option such as -std=c99 is used. -ffast-math enables -fexcess-precision=fast by default regardless of whether a strict conformance option is used.

-fexcess-precision=standard is not implemented for languages other than C. On the x86, it has no effect if -mfpmath=sse or -mfpmath=sse+387 is specified; in the former case, IEEE semantics apply without excess precision, and in the latter, rounding is unpredictable.

Das ist eine ganze Menge. Viele Sachen klingen für mich, als wenn einfach jede Menge if/else Zweige raus genommen werden. Und ich denke diese Optimierungen wirken für double wie für float. Das erklärt noch nicht den Geschwindigkeitsunterschied im letzten Post. Da muss ich wohl mal in den Assembler Code schaun.

https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#Optimize-Options
https://gcc.gnu.org/wiki/FloatingPointMath

Wieder so ein C++ Blog, der aber in Wahrheit scheiß C Code verbreitet. Und dann erzählt er auch noch die Hälfte und macht unvollständige Angaben zu seinen Tests.

Es geht in dem Artikel darum, dass viele denken, dass floats kleiner und schneller als doubles sein. Kleiner sind sie, ja. Verbrauchen nur 4 Byte statt 8 Byte. Aber schneller sind sie nicht unbedingt. Sollten doch unsere heutigen CPUs mit genügen Hardware ausgestattet sein, um auch 64 Bit Zahlen genauso schnell zu verarbeiten. Na, wir werden sehen...

Endergebnis für ungeduldige:
Für kleine Datenmengen und IEEE Gleitkommazahlen verhalten sich bei meinen Tests floats und doubles gleich.
Verzichtet man auf Standard Zahlen und aktiviert fancy Optimierungen, wird man mit einem 1.5 bis 3.5-fachen schnelleren Programm belohnt, welches falsche Ergebnisse liefern kann.
Für große Mengen sieht die Sache anders aus. Dann wird der RAM Bus zum Flaschenhals.

Im ersten Teil übernehme ich sein Code und versuche die selben Ergebnisse zu erzielen. Im zweiten Teil schreibe ich das ganze nach C++ um. Leider gibt er nicht an, wie groß das Array sein soll und wie viele Iterationen gerechnet werden sollen. Also hab ich mal irgendwas eingestellt um ein paar Sekunden Laufzeit zu bekommen. Wenn man wie im original Artikel nur unter einer Sekunde misst, ist die Varianz zu hoch, und man weiß gar nichts.

Testfall:
AMD Athlon(tm) 64 X2 Dual Core Processor 4800+
32Bit Debian
g++ 6.3.0
ARR_SIZE = 200000
NUM_ITER = 20000

Die Standardeinstellung für 64 Bit Systeme mit Optimierung ist SSE. Ich werde daher explizit auf den floating-point Coprozessor mit der Option -mfpmath=387 umschalten. Für SSE werde ich die Option -mfpmath=sse -msse2 benutzen.

Die Option -ffast-math kann übrigens falsche Ergebnisse liefern. Die exakten Vorgaben der IEEE Spezifikationen werden nicht eingehalten. Dafür sollte man ein schnelleres Programm haben.

Mittelwert von drei Berechnungen pro Einstellung. Alle Angaben in Sekunden.

Also, dasss floats schneller sind als doubles, oder dass doubles schneller sind als floats kann ich in diesem Testfall nicht bestätigen. Bei mir sind sie im Grunde gleich schnell. Nur die Option -ffast-math bewirkt eine Geschwindigkeitsverbesserung um das 3-fache bzw. 1.5-fache.

Aber das kann von CPU zu CPU anders sein. Vor allen von microshit zu GNU.

Testfall 2:
Intel(R) Atom(TM) CPU N270 @ 1.60GHz
32 Bit Debian
g++ 4.9.2
ARR_SIZE = 200000
NUM_ITER = 20000

Mittelwert von drei Berechnungen pro Einstellung. Alle Angaben in Sekunden.

Durch den Sparstrumpf Prozessor zeigt sich insgesamt eine höhere Laufzeit. Floats und doubles sind wieder in etwa gleich schnell. Nur die -ffast-math Option macht einen Unterschied. Für floats eine Verbesserung der Laufzeit wieder um das 3-fache. Und für doubles sehe ich sogar eine Verschlechterung um das 1.5-fache. Das müsste man jetzt versuchen zu erklären...

Testfall3:
AMD FX(tm)-9590 Eight-Core Processor
64 bit Debian
g++ 4.9.2
ARR_SIZE = 200000
NUM_ITER = 20000

Mittelwert von drei Berechnungen pro Einstellung. Alle Angaben in Sekunden.

Die Laufzeit fällt insgesamt, verursacht durch mehr GHz Power. Bis auf -ffast-math sind floats und doubles wieder gleich schnell. Der Geschwindigkeitsvorteil für fast-math liegt wieder bei über 3.5-fach für floats und über 1.5-fach für doubles.

Endergebnis:
Für kleine Datenmengen und IEEE Gleitkommazahlen verhalten sich floats und doubles gleich.
Verzichtet man auf Standard Zahlen und aktiviert fancy Optimierungen, wird man mit einem 1.5 bis 3.5-fachen schnelleren Programm belohnt, welches falsche Ergebnisse liefern kann.
Für große Mengen sieht die Sache anders aus. Dann wird der RAM Bus zum Flaschenhals.

Warum verhält sich die Geschwindigkeit mit -ffast-math so unterschiedlich? Das muss man genauer untersuchen.

http://www.bfilipek.com/2012/05/float-vs-double.html

As I told you in my last post, Fortran sucks! This is the C++ Version and its only half as long.

#include <iostream>
#include <limits>
#include <cmath>
using namespace std;

template<typename Func> 
void test(Func func, string funcname, double a, double b, double x) {
  double res = func(a,b);
  cout << funcname << "(" << a << "," << b << ") = " << res << "\n";
  
  if(isnan(x)) {
    if(not isnan(res)) {
      cout << "But should be " << x << "\n";
    }
  } else if(res != x) {
    cout << "But should be " << x << "\n";
  }
}

int main() {
  double NaN = numeric_limits<double>::quiet_NaN();
  test(min<double>, "min", -3.0, NaN, -3.0);
  test(min<double>, "min", NAN, -3.0, -3.0);
  test(min<double>, "min", NAN,  NAN, NAN);
  
  test(max<double>, "max", -3.0, NAN, -3.0);
  test(max<double>, "max", NAN, -3.0, -3.0);
  test(max<double>, "max", NAN,  NAN, NAN);
}

The ugly procedure pointer is replaced with a simple template. But std::min() and std::max() are also template functions. So we must give the compiler a hint which one we want. We want std::min() for doubles. The code is fine, nothing more to say. Except it dosent work as expected ;)

Also Leapfrog ist voll cool. Es berechnet nicht Ort und Geschwindigkeit zum selben Simulationszeitpunkt, sondern versetzt. Frei nach dem Quantenprinzip, dass nicht Ort und Geschwindigkeit gleichzeitig mit beliebiger Genauigkeit gemessen werden kann.

Das Beispiel ist im ekligen C Code überhaupt geschrieben. Zeit es umzuschreiben.
Interessanterweise ist die C++ Version mit t_max=1000000 mit 6.8sec statt 8.4sec sogar um 20% schneller. Ich denke es liegt an den Schleifen die besser optimiert werden können, da die Anzahl der Iterationen fix zur Compilezeit ist. Und diese Information durch die constexpr size() Funktion des std::array Containers auch dem Compiler bekannt ist. Im Gegenzug zur C Version. Wo die Anzahl zwar auch fix zur Compilezeit ist, aber durch eine Laufzeit Variable n im Programm umher gereicht wird. Sonst ist der Code ja gleich.

Anbei die erste C++ Version und der C Code als Vergleich. Mehr wird folgen

inverse_square.c

inverse_square_part1.cpp

[1] https://de.wikipedia.org/wiki/Leapfrog-Verfahren
[2] https://www.physics.drexel.edu/students/courses/Comp_Phys/Integrators/leapfrog/

Habt ihr euch schon immer einmal gefragt, welche Bedingung in einer if() denn letztlich wahr wurde? Hier ist ein einfacher Weg das zu protokollieren.

inline bool test(const bool b, bool& p) {
  p = b;
  return b;
};

void func() {
  double x = 7;

// I thinks this is the best
    bool p1 = x>5;
    bool p2 = x<10;
    if(p1 and p2) {

// Or choose from this
//  bool p1 = false, p2 = false;
//  if(p=x>5 and p2=x<10) {
//  if(test(x>5, p1) and test(x<10, p2)) {
    if(p1) {
      cout << "p1 was true\n";
    }
    if(p2) {
      cout << "p2 was true\n";
    }
  }
}

Thanks to Juan Jimenez.

Diesen Post möchte ich mit einem kürzlich aufgetretenen Problem in Fortran beginnen.
In Fortran gibt es das forall keyword. Damit ist es möglich Schleifen zu bauen, wie die do Schleife, nur dass es möglich ist, forall Schleifen parallel abzuarbeiten.
Eine Vorraussetzung dafür ist folgende: Wenn data der Vector ist, auf den parallel gearbeitet werden soll, dann müssen die Operationen auf data(i) und data(j) mit i != j unabhängig voneinander sein.
Das ist, für das Beispiel, welches ich hier vorstellen möchte, der Fall.

Beispiel: Jeder Punkt im Array soll transformiert, und in einem neuen Array abgespeichert werden. Hier ist die erste native Implementation:

program LLL
  use naaaah
  implicit none
  integer :: i
  type(Point2D), allocatable :: data(:), data2(:)
  real :: x
  
  allocate(data(100))
  allocate(data2(100))
  
  forall(i=1:size(data))
    x = data(i)%x * cos(data(i)%y)
    data2(i)%x = x
    data2(i)%y = data(i)%y
  end forall
end program

Und die Compiler Fehlermeldung

x = data(i)%x * cos(data(i)%y)
Warning: The FORALL with index ‘i’ is not used on the left side of the assignment at (1) and so might cause multiple assignment to this object

Das ist vollkommen richtig, WENN die forall Schleife parallel ausgeführt wird, so wird parallel auf x geschrieben. Aber x ist nur ein scalar. Daher der Fehler. Abhilfe würde schaffen, x als Array zu declarieren. Aber hier gehen bei HPC sofort alle roten Alarmglocken an und die Notbremse wird gezogen.
Okay.
Untersuchen wir jetzt das Stück Code mal nach dem "Was aber nicht Wie" Prinzip. WAS soll der Code machen? Hm keine Ahnung. WIE macht er es? Er multipliziert x in Abhängigkeit von y. Okay das ist alles, was man wissen muss, um den Code umzuschreiben. Nennen wir die neue Funktion, welche obige Multiplikation ausführt als "toLLL". Damit sieht der Code so aus:

pure function toLLL(point) result(LLLPoint)
  implicit none
  type(Point2D), intent(in) :: point
  type(Point2D) LLLpoint
  
  LLLpoint%x = point%x * cos(point%y*DEGRAD)
  LLLpoint%y = point%y
end function

forall(i=1:size(data))
  data2(i) = toLLL(data(i))
end forall

Nun compiliert der Code ohne Fehler! Warum? Weil die vermeintliche Variable x nicht mehr temporär vorkommt. Lass es mich nochmal verdeutlichen. Im alten Code hatten wir so etwas:

xTemp = xalt*fac;
xneu = xTemp

Und im neuen Code haben wir so etwas

xneu= xalt*fac;

Es gibt also keine temporäre Zwischenvariablen die irgendwie von Hand verwaltet werden müsste. PLUS, der Code ist lesbarer geworden. PLUS der Code ist performanter, da er u.U. parallel ausgeführt werden kann. Ich will nicht behaupten, dass durch das Einführen einer neuen Funktion der Code mehr Richtung Funktionale Programmierung geht. Aber er ist definitiv besser. Das wäre es erst, wenn die Funktion toLLL als Funktionsparamter mit übergeben würde. Hier ist die Grenze von Fortran erreicht. Ohne Templates muss für jede generische Variante von Point2D eine neue Funktion toLLL geschrieben werden. Und das inlinen über mehrere Compilations Units hinweg ist auch nicht möglich.

In C++ gibt es natürlich die selbe Art von Problem, aber sie können leicht gelöst werden. Von Hand geschriebene Schleifen funktionieren natürlich immer. Aber, sobald Parallelisierung in's Spiel kommt. Egal auf welcher Art und Weise auch immer, ist es nicht mehr so einfach wie früher. Parallelisierung bedeutet auf viele Kleinigkeiten achten, auf die man vorher nicht achten musste. Aber wir leben ja nicht im Mittelalter. Wie haben Programme, welche die Checks für uns durchführen. Wir müssen sie nur benutzen.
Ein Beispiel dafür wäre std::transform Hier wird keine Schleife mehr explizit aufgeschrieben. Es wird nur noch gesagt WAS getan werden sollen. Das transformieren von data mit Funktion toLLL nach data2.
Hier das komplette C++ Programm

#include <vector>
#include <algorithm>
#include <cmath>
#include <iostream>

using namespace std;

constexpr const double DEGRAD = M_PI/180.0;

struct Point2D {
  double x, y;
};

Point2D toLLL(const Point2D& point) {
  return {point.x*cos(point.y*DEGRAD), point.y};
}
  

int main() {
  vector<Point2D> data(180), data2(180);
  for(size_t i=0; i < 180; ++i) {
    data[i].x = 180;
    data[i].y = -90.0+i;
  }
  
  transform(std::execution::par, begin(data), end(data), begin(data2), toLLL);
  
  for(const Point2D& p : data2) {
    cout << p.x << " " << p.y << "\n";
  }
}

Der ersten Parameter von transform bestimmt, ob der Datensatz parallel oder sequentiell abgearbeitet werden soll. Das Einführen der Funktion toLLL verhindert temporäre Variablen, und die Funktion wird als Parameter übergeben. Das generalisieren der Funktion toLLL auf einem Point Datentyp mit beliebiger Dimension und fundamentalen Type (int,float) ist nun leicht. Beispiel:

template<class Point>
Point toLLL(const Point& point) {
  Point LLLpoint = point;
  LLLpoint[0] *= cos(LLLpoint[1]*DEGRAD);
  return LLLpoint;
}

Die Idee ist von Stephen C. Dewhurst [1] die Abkürzung von mir

FVCRTTPPP - Fancy Variadic Curiously Recurring Template Template Parameter Pack Pattern

Praxisbeispiel:

using FVCRTTPPPScalar = Scalar<int, Ctr, Eq, Rel, Stream>; ! A countable, equal compareable, relational-able, streamable, scalar integer type.
FVCRTTPPPScalar x = 0;

[1] http://stevedewhurst.com/once_weakly/once-weakly20170328/once-weakly20170328.pdf

Der normale dynamic programing Zuker Algorithmus von Faltungen von Sequenzen hat eine Laufzeit von O(n^4), mit zuätzlichen Cache noch O(n^3). Mit der allgemeinen Four-Russians Speedup Methode für dynamic programing Algorithmen ist eine Laufzeit von O(n^3/log n) möglich. Für kleines n also absolut irrelevant. Für großes n halbiert, oder drittelt sogar die Laufzeit. Was dann aber immer noch unheimlich zu lange ist.
Aber was ist die Idee dahinter?
Als erstes wird festgestellt, dass es nicht möglich ist, die Four-Russians Technik auf die klassische Implementation des Zuker anzuwenden. Da die Auswertungsreihenfolge der Matrizen nicht passt. Also wird die Reihenfolge erst mal umgestellt. Erst kommt eine Schleife die nicht von der derzeitigen Spalte j in der Matrix abhängt. Und dann kommt eine 3fach Schleife die abhängt.
In der innersten Schleife findet nun die Optimierung statt. Die Schleife sieht so aus:

S(i,j)=max( S(i,j), S(i,k-1)+S(k,j) )

Für fixes i,j iteriert nur k von i+1 zu j-1. Die Idee ist nun, Indizes zu Gruppe der Größe q zusammen zu fassen. So dass die innerste Schleife nur noch n/q oft läuft.
Nun ja wenn sie meinen. Um das Maximum zu ermitteln, muss man so oder so jeden Wert in der Matrix in diesem Bereich anfassen. Wie die das verhindern wollen, das lese ich aus dem Paper nicht raus.
Zeitverschwendung.

[1] A Simple, Practical and Complete Time Algorithm for RNA Folding Using the Four-Russians Speedup - Yelena Frid and Dan Gusfield

Für die, die es immer noch nicht kapiert haben: Hier noch ein einfachereres Beispiel. Temporäre Objekte, alles klar? Speichert keine Referenzen, es gibt kein Grund das zu tun. Sonst schießt ich euch in den Fuß. In jede Zehe einzeln.

bool alive = true;

struct Foo {
  Foo() { }
  ~Foo() { alive = false; }
};

struct Bar {
  const Foo& dangling;
  
  Bar(const Foo& foo) : dangling(foo) { }
  
  void func() {
    if(alive) {
      cout << "dangling is okay\n";
    } else {
      cout << "dangling is dangling\n";
    }
  }
};

int main() {
  Bar bar{Foo{}};
  bar.func();
  return 0;
}