The Fractal Geometry of Nature

(from Benoid B. Mandelbrot, 1983)

The Fractal Geometry of Nature

(from Benoid B. Mandelbrot, 1983)

Ende der 80er tauchten tolle Bilder auf, die mit einer einfachen Formel erzeugt wurden.

z=z^2+c

Mit dieser Formel und "quadratic dynamics" sollten Computer diese Bilder erzeugen können.

Ich hatte keine Idee wie das gehen sollte!

Computer Kurzweil III

(Spektrum der Wissenschaft, Sonderheft 8, Computer Kurzweil III, 1989)

Ein erster Lichtblick:

n <- 0
WHILE n < 100 und Betrag(z) < 2
  z = z^2 + c
  n <- n + 1
Färbe den Aufpunkt

für Programmiersprachen, die mit Komplexen Zahlen umgehen können.

Computer Kurzweil III

Variante für andere Programmiersprachen

n <- 0
WHILE n < 100 und x^2 + y^2  < 4
  xx  <- x^2 - y^2 + a
  y <- 2xy + b
  x <- xx
  n <- n+1
Färbe den Aufpunkt

Der Algorithmus enthält die Variablen x und y, die jedoch nichts mit den Bildschirm Koordinaten zu tun haben; n wird mit 0 initialisiert, x, y und xx nicht.

(Die Transformation der Bildschirmkoordinaten zu a und b wird getrennt beschrieben.)

Ohne FPU

Pascal Quelltext von 1991. (Quelle: wahrscheinlich c't?)
Großartig für AT286 ohne 287 Coprozessor

{$B-,R-,I-,S-,X+}
USES Graph,Crt;

CONST
      MaxTiefe    =  100;      
      Ymin        = -1200;
      Ymax        =  1200;
      Xmin        = -2000;
      Xmax        =  650;
      Norm        =  1000;

VAR
      MaxX,MaxY,HalbY : INTEGER;

FUNCTION iterate(nx, ny : LONGINT) : LONGINT;

VAR
    tiefe   : INTEGER;
    X,Y,Xalt,
    Xp,Yp   : LONGINT;

BEGIN
    Xp := ((Xmax - Xmin) * nx) DIV (MaxX -1) + Xmin;
    Yp := ((Ymax - Ymin) * ny) DIV (MaxY -1) + Ymin;
    tiefe := 0;
    x := 0;
    y := 0;
    WHILE (tiefe < MaxTiefe) AND ((X * X + Y * Y) < 4000000) DO BEGIN
            Xalt := x;
            X := (X * X - Y * Y) DIV NORM +Xp;
            Y := (2 * Xalt * Y) DIV NORM +Yp;
            INC(tiefe);
      END;
      iterate := tiefe;
END;

VAR
      Gd, Gm        : INTEGER;
      nx,ny,i       : INTEGER;

BEGIN
      Gd := VGA;
      Gm := VGAHi;
      InitGraph(Gd,Gm,'\ps\tp\bgi');
      MaxX := GetMaxX;
      MaxY := GetMaxY;
      HalbY := MaxY DIV 2;
      nx := 0;
      WHILE nx < MaxX DO BEGIN
            ny := 0;
            WHILE ny < HalbY DO BEGIN
                  i := Iterate(nx,ny) MOD 16 + 1;
                  PutPixel(nX,nY,i);
                  PutPixel(nX,MaxY-nY-4,i);
                  INC(ny);
            END;
            INC(nx);
            IF KeyPressed THEN nx := MaxX;
      END;
      ReadKey;
      CloseGraph;
END.

Das Mandelbrot Kata

''Setze den Mandelbrot Algorithmus in der neuen Programmiersprache um.''

(Mein persönliches Kata, für den ersten Kontakt mit einer neuen Programmiersprache)

Mandelbrot Kata

Mandelbrot Kata: lua

start = os.clock()

disp_width = 80
disp_height  = 25
scale_real = 3.0 / disp_width
scale_imag = 3.0 / disp_height
line = ""

for j = 1, disp_height do
  line = ""
  for i = 1, disp_width  do
    x = (i * scale_real) - 1.8
    y = (j * scale_imag) - 1.5
    u = 0.0
    v = 0.0
    u2 = 0.0
    v2 = 0.0
    iter = 0

    while (u2 + v2 < 4.0) and (iter < 100) do
      v = 2 * v * u + y
      u = u2 - v2 + x
      u2 = u*u
      v2 = v*v
      iter = iter + 1
    end

    if iter < 10 then c = ' '
    elseif iter < 20 then c = '.'
    elseif iter < 30 then c = 'o'
    elseif iter < 50 then c = '*'
    elseif iter < 70 then c = 'x'
    elseif iter < 90 then c = 'X'
    then c = '@'
    end

    line = line .. c

  end
  print(line)
  line = ""

end

ende = os.clock()
print("Dauer:" .. ende - start .. "sek.")

Ziel des Mandelbrot Katas

Lernziele:

• Variablen / Typen
• Schleifen
• Verzweigungen
• Einfach Ausgaben
• Laufzeitumgebungsaufrufe (Zeit)
• (Funktionen)

Nebenwirkungen

Wenn man die gleiche Aufgabe mit verschieden Programmiersprachen löst, kann man diese vergleichen.

Computerkurzweil in C++

Eine Implementierung in C++
(basierend auf dem Artikel von 1989)

n <- 0
WHILE n < 100 und Betrag(z) < 2
  z = z^2 + c
  n <- n + 1
Faerbe den Aufpunkt

int pixel(double re, double im)
{
  std::complex<double> z{ 0. }, c{ re, im };
  int n{ 0 };
  while (n < 100 && std::abs(z) < 2.0) {
    z = z * z + c;
    ++n;
  }
  return n;
}

Berechne ein ganzes Bild

Bildschirm-Koordinaten in eine Komplexe Zahlenebene transponieren.

    const int disp_width = 179, disp_height = 80;
    const double start_real = -1.8, start_imag = -1.5;
    const double scale_real = 3.0 / disp_width;
    const double scale_imag = 3.0 / disp_height;

    for (int j = 1; j < disp_height; ++j) {
      for (int i = 1; i < disp_width; ++i) {
        const double re = start_real + i * scale_real;
        const double im = start_imag + j * scale_imag;
        const int iter = pixel(re, im);
        cout << transform(iter);
      }
      cout << '\n';
    }

Bildschirmauflösung 178x79
Wertebereich: [-1.8..1.2; -1.5..1.5] 

i,j: Bildschirmkoordinaten
re,im: Korrdinaten im einer Komplexen Ebende

Transform

Das Ergebnis muss in "Farben" umgewandelt werden.

  char transform(int iter)
  {
    if (iter < 10) return ' ';
    if (iter < 20) return '.';
    if (iter < 30) return 'o';
    if (iter < 50) return '*';
    if (iter < 70) return 'x';
    if (iter < 90) return 'X';
    return '@';
  }

Demo

Projekt: demo

Datei: SimpleExample.h

(check test.cpp)

Demo

Ohne komplexe Zahlen

2 alternative Lösungen ohne std::complex<T>

   // (1 von 1989)
   float u{ 0 }, v{ 0 };
   int iter{ 0 };

   while (u*u + v*v < 4.f && iter < 100) {
     const float uu = v*v - u*u + re;
     u = 2.f * u*v + im;
     v = uu;
     ++iter;
   }

   // (2 intel, ...)
   float u{ 0 }, v{ 0 }, u2{ 0 }, v2{ 0 };
   int iter{ 0 };

   while (u2 + v2 < 4.f && iter < 100) {
     v = 2.f * v * u + im;
     u = u2 - v2 + re;
     u2 = u*u;
     v2 = v*v;
     ++iter;
   }

Vergleich: std::complex vs. float

Windows 10 - Visual Studio (64bit)
Intel(R) Core(TM) i5-5200U CPU @ 2.20GHz

Vergleich: std::complex vs. float

Mint Linux 17.2 - GCC 4.8.5
Intel(R) Core(TM) i5-5200U CPU @ 2.20GHz

Mathematik der Komplexen Zahlen

für C++-Programmierer

Vereinfachte Klasse für Komplexe Zahlen

template<typename FLOAT_TYPE>
struct Complex {

    typedef FLOAT_TYPE value_type;

    FLOAT_TYPE re_;
    FLOAT_TYPE im_;

    FLOAT_TYPE real() const { return re_; }
    FLOAT_TYPE imag() const { return im_; }
};

Complex<double> z = { 0. }, c = { re, im };

Operatoren

Mathematische Operatoren

template<typename FLOAT_TYPE>
inline Complex<FLOAT_TYPE> operator+(const Complex<FLOAT_TYPE>& l, 
                                     const Complex<FLOAT_TYPE>& r)
{
  return Complex<FLOAT_TYPE>{l.re_ + r.re_, 
                             l.im_ + r.im_ };
}

template<typename FLOAT_TYPE>
inline Complex<FLOAT_TYPE> operator*(const Complex<FLOAT_TYPE>& l, 
                                     const Complex<FLOAT_TYPE>& r)
{
  return Complex<FLOAT_TYPE>{ l.re_ * r.re_ - l.im_ * r.im_, 
                              l.re_ * r.im_ + l.im_ * r.re_ };
}

z = z * z + c;

std::abs

Absoluter Betrag.

template<typename T>
T abs(const Complex<T>& z)
{
  return sqrt(z.re_*z.re_ + z.im_*z.im_);
}

while (abs(z) < 2.0 && iter < 100) {
...
}

Beispiel Rechnungen

mit dem Mandelbrot Algorithmus (z = z * z + c)

3 kleine Beispiele für:

• 0, 0i
• -2, 0i
• -1, 0i

(Koordinaten auf der Komplexen Zahlenebene)

Berechne (0, 0i)

Startwert c = (0, 0i)

Complex z{ 0 }, c{ 0.0, 0.0 };
while (std::abs(z) < 2.0 && iter < 100) {
  z = z * z + c;
  ++iter;
}

z = (0, 0i)*(0, 0i) + (0, 0i) = (0, 0i)
bei allen Berechnungen

=> Endlosschleife ohne die Abbruchbedingung iter<100

Berechne (-2, 0i)

Startwert c = (-2, 0i)

Complex z{ 0 }, c{ -2.0, 0.0 };
while (std::abs(z) < 2.0 && iter < 100) {
  z = z * z + c;
  ++iter;
}

z = (0, 0i)*(0, 0i) + (-2, 0i) = (-2, 0i)
abs(-2, 0i) = 2

abs(z) < 2.0 ist im ersten Durchgang erfüllt.

iter == 1

Pythagoras als Kreisformel

In der ersten Iteration erhält z den Wert von c.

Die Kreisformel definiert einen Grenzwert um die Mandelbrotmenge:

Berechne (-1, 0i)

z' = z*z + (-1, 0i)
z^2 = z*z

Der Wert von z wechselt zwischen (0, 0i) und (-1, 0i) und erreicht nie abs(z) >= 2.

=> iter == 100

Regionen des Chaos

Blau: 1 Iteration Rot: 100 Iterationen

Grenzwerte

while (std::abs(z) < 2.0 && iter < 100) {
  z = z * z + c;
  ++iter;
}

Was passiert, wenn man die Konstanten 100 (maximale Iterationen) und 2.0 (Escape Radius) anpasst?

Mehr Informationen

http://mathworld.wolfram.com/MandelbrotSet.html

https://de.wikipedia.org/wiki/Mandelbrot-Menge

https://en.wikipedia.org/wiki/Mandelbrot_set

Optimierungen

Neben den schönen Bildern gilt die Hauptfrage bei der Mandelbrotmengen-Berechnen nach der Performance. Wie kann man die Bilder möglichst schnell berechnen?

• while Schleife
• loop unroll
• Literale
• Schritte berechnen
• Aus den Schleifen herausnehmen
• teure Operationen (e.g. sqrt, abs)

Optimierung der Schleife

    ...
    std::complex z{0.}, c{x,y};
    int iter{0};    

    while ((std::abs(z) < 2.0) && (iter < 100)) {
        z = z * z + c;
        ++iter;
    }
    ...

while -Abfrage ist immer true beim ersten Test!

Optimierung der Schleife

// `while` - loop
int pixel(double x, double y) 
{
    std::complex z{0.}, c{x,y};
    int iter{0};    
    while ((std::abs(z) < 2.0) && (iter < 100)) {
        z = z * z + c;
        ++iter;
    }        
    return iter;
}

// `do` - loop
int pixel(double x, double y) 
{
    std::complex z{0.}, c{x,y};
    int iter{0};    
    do {
        z = z * z + c;
        ++iter;
    } while ((std::abs(z) < 2.0) && (iter < 100));

    return iter;
}

Optimierung der Schleife

Wie kann die Optimierung gemessen werden?

• Blauer Bereich
  • 1 statt 2 Abfragen
• Roter Bereich
  • 100 statt 101 Abfragen

Optimierung der Schleife

Alle Werte in MPixel/s
GCC 4.8.5 auf Mint Linux 13.2, VS2015 unter Windows 10

Schleifen ausrollen

Was bringt es Schleifen auszurollen? (loop unrolling)

Demo

Loop unrolling

Alle Werte in MPixel/s
GCC auf Mint Linux 13.2, VS2015 unter Windows 10

Literale

  template<typename FLOAT_TYPE, int MAX_ITER>
  struct Pixel
  ...
    do {
        v = 2.0 * v * u + im;
        u = u2 - v2 + re;
        u2 = u*u;
        v2 = v*v;
        ++iter;
      } while (u2 + v2 < 4.0 && iter < MAX_ITER);
  ...

Welches Literal soll für "2.0" verwendet werden?

• "2"
• "2.f"
• "2.0"

Benchmark

Raus aus der Schleife

Alles aus der Schleife heraus ziehen was möglich ist!

for (int j = 1; j < disp_height; ++j) {
  for (int i = 1; i < disp_width; ++i) {
    const double re = i * scale_real - 1.8;
    const double im = j * scale_imag - 1.5; // (A)
    const int iter = pixel(re, im);
    ...
  }
}

for (int j = 1; j < disp_height; ++j) {
  const double im = j * scale_imag - 1.5;  // (A)
  for (int i = 1; i < disp_width; ++i) {
    const double re = i * scale_real - 1.8; 
    const int iter = pixel(re, im);
    ...
  }
}

Schritt Berechnung

Startwert vs. Addition

for (int j = 1; j < disp_height; ++j) {
  const double im = y1 + j * scale_imag; // (1)
  for (int i = 1; i < disp_width; ++i) {
    const double re = x1 + i * scale_real; // (2)
    const int iter = pixel(re, im);
    ...
  }
}

double im = y1; // (1a)
for (int j = 1; j < disp_height; ++j) {
  im += scale_imag; // (1b)
  double re = x1;   // (2a)
  for (int i = 1; i < disp_width; ++i) {
    re += scale_real; // (2b)
    const int iter = pixel(re, im);
    ...
  }
}

Schritt Berechnung

Von Startwert neu berechnen: Der Fehler bleibt konstant / in der gleichen Größenordnung.

for(int j = 0; ...) {
    const float im = j * scale_imag + starty;
    ...
}

Bei der Addition summieren sich die numerischen Fehler auf!

float im = starty;
for(int j = 0; ...) {
    im += scale_imag;
    ...
}

Schritt Berechnung

• Die Addition ist schneller als die Index Berechnung
• Die Addition summiert Fehler auf.
  • Probleme bei der Unit tests und textueller Darstellung

Vermeide (teure) Operationen

z.B. abs() bzw. sqrt():

do {
...
} while ((std::abs(z) < 2.0) && (iter < 100));

do {
...
} while ((z.imag() * z.imag() + z.real() * z.real() < 4.0) && (iter < 100));

Vermeide (teure) Operationen

Vermeide (teure) Operationen

std::sqrt(re*re+im*im) >= 2.0 ist merklich teuer als re*re+im*im >= 4.0

std::hypot() ist teurer als std::sqrt() wegen zusätzlicher Checks.

std::abs(std::complex) ist teurer als std::hypot() wegen zusätzlicher NaN, Unendlich und Stabilitäts-Checks.

Datentypen

Vergleiche verschiedene Datentypen

Warum ist int so langsam?

Wieso ist int langsamer als Fließkomma-Typen?

Warum ist int so langsam?

int ist nicht das Problem - es is die Division-Korrektur!

Warum ist long double so langsam (auf Linux)?

Der GCC erzeugt für alle Datentypen schnelleren Code - außer für long double. Warum?

Warum ist long double so langsam (auf Linux)?

Antwort A:

Linux: sizeof(long double) == 16

Windows (Visual Studio) sizeof(long double) == 8

Kein SIMD für long double

Antwort B:

SSE unterstützt nur 4 und 8 Byte Fließkommatypen.

long double muss auf der FPU gerechnet werden (Linux).

Die SSE Recheneinheiten sind schneller als die FPU.

Ist std::complex langsam?

• http://www.lomont.org/Math/Papers/2011/Intro%20to%20Intel%20AVX-Final.pdf

Kein Problem - AVX ist viel schneller

• CppCon 2015: André Bergner “Faster Complex Numbers”
• http://scicomp.stackexchange.com/questions/1998/why-would-a-computational-scientist-need-to-implement-their-own-version-of-std

Nutze GCC mit -ffast-math

Complex types

Alle Werte in MPixel/s

Zusammenfassung

• Es gibt keine Zaubertricks um den Compiler um Längen zu schlagen.
  • Loop unroll, ...
• Es gibt viele Möglichkeiten den Compiler auszubremsen.
  • Float-Literal, ...
• Die Wahl des richtigen Datentyps ist entscheidend.
  • siehe SIMD
• Eine Klasse / Bibliothek mit einem anderen Schwerpunkt kann das Programm wesentlich verlangsamen
  • std::complex / NaN-Prüfungen

Multithreading

Siehe:

• "The free lunch is over" oder
• Welcome to the Jungle

beide: Herb Sutter.

• OpenMP
• Intel TBB
• Performance Vergleich / Analyse

OpenMP

Äußere Schleife parallelisieren:

#pragma omp parallel for
for (int j = 0; j < height; ++j) {
    const FLOAT_TYPE y = (j * dy) + y1;
    for (int i = 0; i < width; ++i) {
        const FLOAT_TYPE x = (i * dx) + x1;
        auto iter = pixel(x, y);
        ...
    }
}

• #pragma omp OpenMP Befehl
• parallel for Parallele for-Schleife, muss vor einer for Anweisung stehen

Intel TBB

Äußere Schleife durch tbb::parallel_for ersetzen.

tbb::parallel_for(int(0), int(height), [&](int j) {
    const FLOAT_TYPE y = (j * dy) + y1;
    for (int i = 0; i < width; ++i) {
        const FLOAT_TYPE x = (i * dx) + x1;
        const auto iter = pixel(x, y);
        ...
    }
});

Der Schleifen-body (die innere Schleife) wird als Lambda formuliert.

Intel TBB vs OpenMP

OpenMP

#pragma omp parallel for
for (int j = 0; j < height; ++j) {
    ...
}

Intel TBB

tbb::parallel_for(int(0), int(height), [&](int j) {
    ...
});

Der innere Teil (...) ist in beiden Varianten gleich.

Performance Vergleich

Visual Studio 2015.

Erste Analyse

2 Ergebnisse:

1. Die Beschleunigung ist nicht Factor 4 (12.6 x 4 = 50.4)
2. OpenMP ist scheinbar viel langsamer als TBB/std::futures

Multi-Core Archithektur

• Was bedeutet Hyper-Threading?
• Was bedeutet Intel Turbo Boost?

Hyper-Threading Tests

Intel Core i7-2860QM 2.5 GHz 4 Kerne, HyperThreading (8 Hardware-Threads)

HT: Hyper-Threading
TB: Intel Turbo Boost
Single threaded vs TBB

Warum ist OpenMP so langsam

OpenMP verwendet als default(*) einen statischen Scheduler, dazu wird der Bereich in gleich große Blöcke aufgeteilt.

Der dynamischen Scheduler verteilt die Arbeit "gleichmäßig" auf die Kerne.

#pragma omp parallel for schedule(dynamic)

(*) unspezifiziert. (aber alle mir bekannten Compiler nutzen den statischen Scheduler als default Scheduler)

Warum ist OpenMP so langsam

Performance Vergleich

Benchmarks and Multithreading

Demo

1. Erzeugen des Thread pools
2. Temperaturabhängige Messungen

Mehr C++

• Rekursion
• Ranges

Rekursion

Schleifen und Rekursion sind das Gleiche - Man kann sie einfach ineinander überführen.

aus der Sicht eines Mathematikers

Kann man den Mandelbrot Mengen Algorithmus effizient rekursiv programmieren?

Ein erster Versuch

int pixel_detail(int iter, COMPLEX_TYPE z, 
                 const COMPLEX_TYPE c, const int max_iter) 
{
  if (iter < max_iter) {
    z = z * z + c;
    if ((z.imag() * z.imag() + z.real() * z.real() < FLOAT_TYPE(4.0)))
      return pixel_detail(iter + 1, z, c, max_iter);
  }
  return iter;
}

Da war doch noch was Endrekursion bzw. tail recursive optimization.

Endrekursion

Eine rekursive Funktion f ist endrekursiv (englisch tail recursive; auch endständig rekursiv, iterativ rekursiv, repetitiv rekursiv), wenn der rekursive Funktionsaufruf die letzte Aktion zur Berechnung von f ist. Vorteil dieser Funktionsdefinition ist, dass kein zusätzlicher Speicherplatz zur Verwaltung der Rekursion benötigt wird.

Source Wikipedia: https://de.wikipedia.org/wiki/Endrekursion

Endrekursion

Rekursive Funktionen sollten vom Compiler effektiv optimiert werden, wenn:

• der rekursive Aufruf in der letzten Zeile erfolgt und
• alle lokalen Variablen als Parameter übergeben werden.

Tail Recursive Optimization

int pixel_detail(int iter, COMPLEX_TYPE z, 
                 const COMPLEX_TYPE c, const int max_iter) 
{
  if (iter == max_iter) 
    return max_iter;
  z = z * z + c;
  if ((z.imag() * z.imag() + z.real() * z.real() >= FLOAT_TYPE(4.0)))
    return iter;

  return inner(iter + 1, z, c, max_iter);
}

Keine lokalen Variablen,
die letzte Zeile enthält den rekursiven Aufruf.

Compiler-Vergleich

mandel::Complex statt std::complex

Zusammenfassung

tail recursive optimization

Nette Idee - kann funktionieren - muss aber nicht

Keine Schleifen

Man soll in C++ keine Schleifen selbst programmieren, sondern (STL-)Algorithmen verwenden!

Für den Mandelbrot Algorithmus ist das schwierig, da die Anzahl der Iterationen nicht vorhersehbar ist.

Aber die for-Schleifen über alle Pixel (den Bildschirm) sollte funktionieren.

Berechne ein ganzes Bild (Wdh.)

Überarbeitetes Beispiel:

    const int disp_width = 179, disp_height = 80;
    const double start_real = -1.8, start_imag = -1.5;
    const double scale_real = 3.0 / disp_width;
    const double scale_imag = 3.0 / disp_height;

    for (int j = 1; j < disp_height; ++j) {
      for (int i = 1; i < disp_width; ++i) {
        const double re = start_real + i * scale_real;
        const double im = start_imag + j * scale_imag;
        const int iter = pixel(re, im);
        cout << transform(iter);
      }
      cout << '\n';
    }

Bildschirmauflösung 178x79
Wertebereich: [-1.8..1.2; -1.5..1.5] 

i,j: Bildschirmkoordinaten
re,im: Korrdinaten im einer Komplexen Ebende

Range über alle Pixel

Randbedingungen

• Ausschnitt der Komplexen Zahlenebene.
• Ein Pixel ist ein Wert der Komplexen Zahlenebene.
• Der Ergebnis wird in einem std::vector gespeichert.

(logische) Erweiterung von boost::range::iterator_range bzw. irange

template<class Integer, class StepSize>
iterator_range< range_detail::integer_iterator_with_step<Integer, StepSize> >
irange(Integer first, Integer last, StepSize step_size);

parallel stl / array_view

2D Integer range from array_view proposal

std::experimental::D4087::bounds<2> bnd{ width, height };
std::experimental::parallel::for_each(
    std::experimental::parallel::par, std::begin(bnd), std::end(bnd), 
        [&](std::experimental::D4087::index<2> idx) {

        const auto i = idx[0];
        const auto j = idx[1];
        const double re = start_real + i * scale_real;
        const double im = start_imag + j * scale_imag;
        const auto iter = pixel(val.real(), val.imag());
        const color = transform(iter); // num iterationen -> color
        ...
}

Parallel STL on CodePlex bzw. das array_view proposal:
Multidimensional bounds, offset and array_view, revision 7

Beispiel

ComplexRange2D<complex<double>> range{ -1.0, 1.0, -2.0, 2.0, 3, 5 };
for (const auto v : range) {
    cout << v << " ";
}
cout << endl;

Wertebereich: [-1.0..1.0; -2.0..2.0]
Bildschirmauflösung 3x5 Pixel

(-1,-2) (0,-2) (1,-2) (-1,-1) (0,-1) (1,-1) (-1,0) (0,0) (1,0) (-1,1) (0,1) (1,1) (-1,2) (0,2) (1,2)

Mandelbrot ala boost::irange

using cmplx = complex<double>;
Color t(int iteration);  // num iterationen -> color
vector<Color> picture(width*height);
...
ComplexRange2D<cmplx> range(
    start_real, end_real,  
    start_imag, end_imag, 
    width, height);

transform(begin(range), end(range), begin(picture), 
    [&pixel, &t](cmplx val) {
        const auto iter = pixel(val.real(), val.imag());
        return t(iter);
});

Alles ausser SIMD

Warum kein SIMD?

SIMD ist zu wichtig (für die Mandelbrot-Berechnung) um in einem Unterabschnitt abgehandelt zu werden.

Teaser:

• Rust ist schneller als C++
  benchmarkgame, mandelbrot
  
• Benchmark

The Computer Language Benchmarks Game

http://benchmarksgame.alioth.debian.org/u64q/performance.php?test=mandelbrot

SIMD Performance Vergleich

Visual Studio 2015 | i5-5200U CPU @ 2.20GHz
c++ vector class by Agner Fog
(*) Ergebnis unerwartet