Certains langages de programmation ont un support natif pour les nombres complexes. Malheureusement ce n'est pas le cas de Java (et par conséquent de Processing). C'est pourquoi nous allons développer une classe permettant de faire des calculs avec des nombres complexes. Notre point de départ est :

class Complexe {
  float a; // partie réelle
  float b; // partie imaginaire
}

Exercice 20.1. Ajoutez quelques constructeurs à votre classe :

• Un constructeur qui crée un nombre complexe dont les parties réelle et imaginaire sont données en argument ;
• Un constructeur qui crée un nombre complexe dont la partie réelle est donnée en argument et la partie imaginaire est 0 ;
• Un constructeur sans argument (constructeur par défaut) qui crée le nombre complexe 0 ;
• Un constructeur qui crée un nombre complexe qui est égal à un autre nombre complexe donné en argument (constructeur par recopie).

Essayons de créer et afficher un nombre complexe :

Complexe z = new Complexe(2, 3.5);
println(z);

Malheureusement nous allons voir quelque chose qui ressemble à :

sketch_180323a$Complexe@72e76b50

Mais heureusement si on ajoute une méthode toString() à une classe, on peut afficher les objets sous la forme qu'on souhaite.

Exercice 20.2. Ajoutez la méthode suivante à votre classe :

String toString() {
  return a + "+" + b + "i";
}

Réessayez l'affichage. Améliorez cette méthode pour qu'elle produise une écriture plus lisible dans certains cas particuliers. Par exemple 0.0 à la place de 0.0+0.0i, 1.0-i à la place de 1.0+-1.0i etc. Testez votre méthode avec le code suivant :

Complexe[] z = new Complexe[10];
int i = 0;
for (int a = 0; a <= 1; a++) {
  for (int b = -2; b <= 2; b++) {
    z[i] = new Complexe(a, b);
    i++;
  }
}

for (i = 0; i < 10; i++) {
  println(z[i]);
}

Il doit afficher

-2.0i
-i
0.0
i
2.0i
1.0-2.0i
1.0-i
1.0
1.0+i
1.0+2.0i

Exercice 20.3. On peut accéder aux parties réelle et imaginaire d'un nombre z en utilisant z.a et z.b. Mais souvent accéder directement aux variables d'un objet n'est pas une bonne idée. Java a même un mécanisme pour empêcher cela :

class Complexe {
  private float a;
  private float b;
}

Déclarer les variables comme private ne vous empêchera pas d'y accéder en Processing, mais il faut s'habituer à éviter cette pratique. Déclarez deux méthodes re() et im() qui renvoient les parties réelle et imaginaire (accesseurs). On utilisera ces méthodes à la place de l'accès direct aux champs :

Complexe z = new Complexe(2, 3);
float a = z.re();
println("La partie réelle de z est", a);
float b = z.im();
println("La partie imaginaire de z est", b);

Exercice 20.4. Écrivez une méthode boolean egal(Complexe z) qui compare 2 nombres complexes. Deux nombres complexes sont égaux si et seulement si leurs parties réelles sont égales et leurs parties imaginaires sont égales. Testez avec le code suivant :

Complexe z1 = new Complexe(1.5, 2);
Complexe z2 = new Complexe(z1);
Complexe z3 = z1;

println(z1.egal(z1));
println(z1 == z1);
println(z1.egal(z2));
println(z1 == z2);
println(z1.egal(z3));
println(z1 == z3);

Expliquez les résultats.

Exercice 20.5. Écrivez une méthode boolean egal(float c) qui compare un nombre complexe et un nombre réel. Un nombre complexe z est égal à un nombre réel c si la partie réelle de z est égale à c et la partie imaginaire de z est nulle.

Complexe z1 = new Complexe(2);
Complexe z2 = new Complexe(z1.re(), 1);
println(z1.egal(2)); // true
println(z2.egal(2)); // false

Pour calculer la somme de deux nombres complexes, on ajoute à notre classe la méthode suivante :

Complexe plus(Complexe z) {
  float aSomme = a + z.a;
  float bSomme = b + z.b;
  Complexe somme = new Complexe(aSomme, bSomme);
  return somme;
}

Notez bien que cette méthode crée un nouvel objet qui est la somme de l'objet qui l'appelle (argument implicite) et l'objet passé en argument (argument explicite). Ensuite elle renvoie ce nouvel objet. Voici comment on peut l'utiliser.

Complexe z1 = new Complexe(1.5, 2);
Complexe z2 = new Complexe(2.5, -2);
Complexe z3 = z1.plus(z2);
println(z1, " + ", z2, " = ", z3); // 4.0

// La somme doit être commutative
if (z2.plus(z3).egal(z3.plus(z2))) {
  println("Tout va bien !");
} else {
  println("Nous avons un sérieux problème");
}

Notez aussi qu'on peut contracter notre méthode en une seule ligne :

Complexe plus(Complexe z) {
  return new Complexe(a + z.a, b + z.b);
}

Exercice 20.6. Sur le même modèle, écrivez une méthode qui calcule la somme d'un nombre complexe et un nombre réel.

Complexe z1 = new Complexe(1.5, 2);
Complexe z2 = new Complexe(2.5, 1);
Complexe z3 = z1.plus(z2);
println(z3); // 4.0 + 3.0i
Complexe z4 = z3.plus(2);
println(z4); // 6.0 + 3.0i

Exercice 20.7. Attaquez-vous à la soustraction et la multiplication (par nombre complexe et par nombre réel). Testez.

Chacune de ces méthodes renvoie un objet Complexe. L'avantage de cette approche est que vous pouvez écrire des expressions arithmétiques en enchaînant les appels. Par exemple

z₄ = z₁(3z₂ + z₁z₃)

peut s'écrire ainsi :

z4 = z1.mult(z2.mult(3).plus(z1.mult(z3)));

Exercice 20.8. Que signifie l'expression suivante ?

z1.mult(z2).mult(5).moins(4).mult(z1.mult(z1).plus(z2.mult(2)))

Exercice 20.9. Écrivez des méthodes qui calculent le module et le conjugué d'un nombre complexe. Testez.

Complexe z = new Complexe(1, -2);

// Le conjugué du conjugué de z doit être égal à z
if (z.conj().conj().egal(z)) {
  println("Tout va bien !");
} else {
  println("Nous avons un sérieux problème");
}

// z multiplié par son conjugé est égal au carré du module de z
if (z.mult(z.conj()).egal(z.mod() * z.mod())) {
  println("Tout va bien !");
} else {
  println("Nous avons un sérieux problème");
}  

Source xkcd

Exercice 20.10. Implantez la division par un nombre complexe et par un nombre réel.

Complexe z1 = new Complexe(1);
Complexe z2 = new Complexe(1, 1);
println(z1.div(z2)); // 0.5-0.5i

Exercice 20.11. (facultatif) Écrivez une méthode qui renvoie l'argument d'un nombre complexe et une méthode qui élève un nombre complexe à une puissance entière.

Complexe z = new Complexe(1, 1);
println(z.puissance(5)); // -4.0-4.0i (à perte de précision près)

Après avoir bien travaillé, nous pouvons enfin commencer à nous amuser. Rappelons qu'un nombre complexe peut être représenté par un point dont les coordonnées sont la partie réelle et la partie imaginaire du nombre. Nous allons utiliser les nombre complexes pour dessiner des fractales, ces objets mathématiques fascinants dont une des caractéristiques principales est l’auto-similarité.

Source chainsawsuit

Dans un premier temps nous allons nous intéresser au comportement de la suite

z₀ = 0

z_n+1 = (z_n)² + c

Plus précisément, nous chercherons à savoir si elle diverge ou reste bornée pour différentes valeurs de c. Nous allons tracer la trajectoire de cette suite. On ne peut pas itérer jusqu'à l'infini, c'est pourquoi nous allons calculer seulement les premiers n_max membres. On peut montrer que si au bout d'un moment la suite sort du cercle de rayon 2 centé à l'origine, elle est divergente. On considère qu'elle ne l'est pas si au bout de n_max itérations elle reste dans le cercle. Voici donc notre plan :

• Setup :
  • Choisir un c aléatoire
  • n = 0
  • z₀ = 0
• Draw
  • Calculer z_n+1
  • Tracer une ligne entre z_n et z_n+1
  • Incrémenter n
  • Si n dépasse n_max ou z_n sort du cercle de rayon 2, s'arrêter

Pour établir la correspondance entre la zone du plan complexe qui nous intéresse et les coordonnées des pixels de la fenêtre, nous allons faire appel à notre vieille amie, la fonction map()

Exemple 20.1. Escape game

// rectangle du plan complexe qu'on veut afficher sur la fenêtre
Complexe zMin = new Complexe(-2.1, -2.1);
Complexe zMax = new Complexe(2.1, 2.1);

int nMax = 100; // nombre maximal d'itérations

Complexe c;
Complexe z;
int n;

int xFenetre(float re) {
  return int(map(re, zMin.re(), zMax.re(), 0, width));
}

int yFenetre(float im) {
  return int(map(im, zMin.im(), zMax.im(), height, 0));
}

void setup() {
  size(400, 400);
  background(255);

  // dessiner les axes
  int x0 = xFenetre(0);
  int y0 = yFenetre(0);
  stroke(0);
  line(0, y0, width, y0);
  line(x0, 0, x0, height);

  // dessiner un cercle de rayon 2
  int x2 = xFenetre(2);
  int d = 2 * (x2 - x0);
  noFill();
  ellipse(x0, y0, d, d);

  // Générer un c aléatoire
  c = new Complexe(random(-1, 1), random(-1, 1));
  fill(0);
  textAlign(LEFT, TOP);
  textSize(12);
  text("c = " + c.toString(), 0, 0);

  // z(0) = 0
  n = 0;
  z = new Complexe();
}

void draw() {
  // z(n+1) = z(n)^2 + c
  Complexe zSuiv = z.mult(z).plus(c);

  // tracer une ligne entre z(n) et z(n+1)
  stroke(255, 0, 0);
  line(xFenetre(z.re()), yFenetre(z.im()), xFenetre(zSuiv.re()), yFenetre(zSuiv.im()));

  n++;
  z = zSuiv;

  if (n == nMax || z.mod() > 2) {
    println(n);
    noLoop();
  }
}

En exécutant ce sketch plusieurs fois, on se rend compte que pour certaines valeurs de c la suite quitte très rapidement le cercle :

Pour d'autres valeurs elle met plus de temps pour s'échapper :

Et pour certaines valeurs de c la suite reste dans le cercle au bout de 100 membres et elle n'a pas l'air de vouloir le quitter.

L'ensemble de Mandelbrot est l'un des exemples les plus connus de visualisation mathématique. La définition de cet ensemble est la suivante :

L'ensemble de Mandelbrot est l'ensemble de nombres complexes c pour lesquels la suite

z₀ = 0

z_n+1 = (z_n)² + c

ne diverge pas.

La première visualisation de cet ensemble a été publiée par Robert W. Brooks et Peter Matelski en 1978 :

Source Wikipédia

40 ans plus tard, nous avons des moyens à faire des visualisations en plus haute résolution. Voici comment on va s'y prendre.

• Pour chaque pixel (x, y) de la fenêtre
  • Calculer le nombre complexe c correspondant à (x, y)
  • Itérer la suite z_n avec ce c pour déterminer si elle est divergente ou pas
  • Si la suite n'est pas divergente, colorer le pixel en noir

Exemple 20.2. L'ensemble de Mandelbrot en noir et blanc

// rectangle du plan complexe qu'on veut afficher sur la fenêtre
Complexe zMin = new Complexe(-2.5, -1);
Complexe zMax = new Complexe(1, 1);

int nMax = 100; // nombre maximal d'itérations

void setup() {
  size(840, 480);
  dessinerMandelbrot();
}

void dessinerMandelbrot() {
  background(255);
  stroke(0);
  // Pour chaque pixel (x, y) de la fenêtre
  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      // Calculer le nombre complexe correspondant à (x, y)
      Complexe c = zPlan(x, y);
      // Itérer la suite
      int n = itererSuite(c);
      // Si la suite n'a pas divergé, colorer le pixel en noir
      if (n == nMax) {
        point(x, y);
      }
    }
  }
}

// Transforme un pixel de la fenêtre en point du plan complexe
Complexe zPlan(int x, int y) {
  float re = map(x, 0, width, zMin.re(), zMax.re());
  float im = map(y, 0, height, zMin.im(), zMax.im());
  return new Complexe(re, im);
}

// Itère la suite z(0) = 0, z(n+1) = f(z(n)) pendant nMax itérations
// ou jusqu'à ce qu'elle sorte du cercle de rayon 2
// Renvoie le nombre d'itérations effectuées
int itererSuite(Complexe c) {
  Complexe z = new Complexe();
  int n = 0;
  while (n < nMax && z.mod() < 2) {
    z = f(z, c);
    n++;
  }
  return n;
}

Complexe f(Complexe z, Complexe c) {
  return z.mult(z).plus(c);
}

Exercice 20.12. Pour les points c proches de la frontière mais n'appartenant pas à l'ensemble la séquence met beaucoup de temps pour s'échapper. Ainsi plus nMax sera grand, plus le rendu sera précis, mais plus les calculs seront lents. Testez différentes valeurs de nMax et observez les différences de rendu.

Exercice 20.13. Pour obtenir un rendu plus joli, on peut colorer les pixels n'appartenant pas à l'ensemble en fonction du nombre d'itérations que la suite correspondante a mis pour s'échapper. Il y a plusieurs façons de faire cela. Une façon très simple est de transformer le nombre d'itérations (entre 0 et nMax) en niveau de gris (entre 0 et 255).

float gris = map(n, 0, nMax, 0, 255);

Utilisez cette méthode pour colorer les pixels à l'extérieur de l'ensemble.

Exercice 20.14. Le gris c'est bien, mais les autres couleurs sont encore mieux. On peut utiliser le nombre d'itérations pour déterminer la teinte de la couleur des pixels. Utilisez

colorMode(HSB, nMax, 1, 1);

et ensuite

stroke(n, 1, 1);

pour les pixels n'appartenant pas à l'ensemble.

Même à l’œil nu on peut constater que le rendu de l'image n'est pas immédiat et prend un certain temps.

Exercice 20.15. Utilisiez millis() pour mesurer le temps d'exécution de la fonction dessinerMandelbrot().

Sur notre machine le temps de rendu est environ 1.5 secondes. Cela va poser des problèmes si on veut animer notre sketch pour zoomer sur différentes régions de l'ensemble par exemple. Avant d'aller plus loin, prenons un peu de temps pour optimiser notre code.

Utiliser la fonction point() pour colorer les pixels est pratique mais pas très efficace. Heureusement Processing nous donne accès à un tableau pixels qui stocke les couleurs de tous les pixels de la fenêtre. Travailler directement avec ce tableau est beaucoup plus rapide.

Les pixels vivent en deux dimensions mais pour des raisons d'efficacité Processing les stocke en ordre linéaire dans un tableau à une dimension.

Pour trouver l'indice d'un pixel (x, y) dans le tableau pixels il faut calculer le nombre de pixels dans les lignes précédentes (y * width) et ajouter le nombre de pixels avant lui dans sa ligne (x) :

int i = y * width + x;

Par exemple l'indice du pixel (3, 2) de la figure ci-dessus sera 2 * 5 + 3 = 13.

Avant d'utiliser le tableau de pixels, il faut le charger en appelant la fonction loadPixels(). Après avoir fini, il faut signaler à Processing qu'il peut les dessiner en appelant la fonction updatePixels().

Exemple 20.3. Travailler avec les pixels

void setup() {
  size(400, 400);
}

void draw() {
  // Processing, passe-moi les pixels, j'ai besoin de leur parler
  loadPixels();

  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      float d = dist(x, y, mouseX, mouseY);
      int i = y * width + x;
      pixels[i] = color(d);
    }
  }

  // J'ai fini, tu peux les récupérer et les dessiner
  updatePixels();
}

Dans le cas où on parcourt les pixels dans l'ordre dans lequel ils stockés dans le tableau, on peut se passer de la formule y * width + x :

int i = 0;
for (int y = 0; y < height; y++) {
  for (int x = 0; x < width; x++) {
    float d = dist(x, y, mouseX, mouseY);
    pixels[i] = color(d);
    i++;
  }
}

Exercice 20.16. Reprenez le code qui dessine l'ensemble de Mandelbrot. Utilisez directement le tableau pixels au lieu de la fonction point(). Mesurez à nouveau le temps de rendu. Combien avez-vous gagné ?

Exercice 20.17. Pour vérifier si la séquence est sortie du cercle de rayon 2, on utilise la méthode mod() de la classe Complexe. Cette méthode calcule une racine carrée, une opération relativement coûteuse. On peut éviter ce calcul en ajoutant à notre classe Complexe une méthode mod2() qui calcule le carrée du module et qui n'a pas besoin de calculer la racine carrée. Ensuite on peut remplacer le test de divergence par z.mod2() < 4. Cette optimisation nous permettra de gratter encore quelques millisecondes précieuses. Implantez-la et mesurez le gain de temps.

On va rendre notre sketch interactif. L'idée est qu'utilisateur puisse zoomer / dézoomer et se déplacer dans le plan complexe pour explorer l'ensemble de Mandelbrot. Mais avant d'attaquer cette tâche, prenons un peu de temps pour faire du propre. Mettons notre code dans une classe Mandelbrot. Elle regroupe les fonctions que nous avons développé jusqu'à maintenant. En plus, on profite de l'occasion pour s'assurer que l'échelle sera la même sur les deux axes et que l'image ne sera pas déformée quelle que soit la taille de la fenêtre. Ainsi notre programme principal tiendra en quelques lignes.

Exemple 20.4. Mandelbrot OO

Mandelbrot mandel;

void setup() {
  size(840, 480);
  mandel = new Mandelbrot(-2.5, -1, 3.5, 100);
  int start = millis();
  mandel.dessiner();
  int stop = millis();
  println("Image rendue en", stop - start, "millisecondes");
}

class Mandelbrot {
  // rectangle du plan complexe qu'on veut afficher sur la fenêtre
  Complexe zMin;
  Complexe zMax;

  int nMax; // nombre maximal d'itérations

  Mandelbrot(float reMin, float imMin, float largeur, int nMax) {
    float reMax = reMin + largeur;
    // Calculé en fonction de la taille de la fenêtre
    // pour garder la même échelle sur les deux axes
    float imMax = imMin + largeur * height / width;

    zMin = new Complexe(reMin, imMin);
    zMax = new Complexe(reMax, imMax);
    this.nMax = nMax;
  }

  void dessiner() {
    colorMode(HSB, nMax, 1, 1);
    loadPixels();
    // Pour chaque pixel (x, y) de la fenêtre
    int i = 0;
    for (int y = 0; y < height; y++) {
      for (int x = 0; x < width; x++) {
        // Calculer le nombre complexe correspondant à (x, y)
        Complexe c = zPlan(x, y);
        // Itérer la suite
        int n = itererSuite(c);
        // Si la suite n'a pas divergé, colorer le pixel en noir
        // sinon, colorer en fonction du nombre d'itérations
        if (n == nMax) {
          pixels[i] = color(0);
        } else {
          pixels[i] = color(n, 1, 1);
        }
        i++;
      }
    }
    updatePixels();
  }

  // Transforme un pixel de la fenêtre en point du plan complexe
  Complexe zPlan(int x, int y) {
    float re = map(x, 0, width, zMin.re(), zMax.re());
    float im = map(y, 0, height, zMin.im(), zMax.im());
    return new Complexe(re, im);
  }

  // Itère la suite z(0) = 0, z(n+1) = f(z(n)) pendant nMax itérations
  // ou jusqu'à ce qu'elle sorte du cercle de rayon 2
  // Renvoie le nombre d'itérations effectuées
  int itererSuite(Complexe c) {
    Complexe z = new Complexe();
    int n = 0;
    while (n < nMax && z.mod2() < 4) {
      z = f(z, c);
      n++;
    }
    return n;
  }

  Complexe f(Complexe z, Complexe c) {
    return z.mult(z).plus(c);
  }
}

Maintenant que notre code est bien organisé, on va ajouter à la classe Mandelbrot des méthodes permettant de zoomer autour d'un point et de se déplacer dans le plan complexe.

void zoom(int x, int y, float factor) {
  Complexe z = zPlan(x, y);
  zMin = z.plus(zMin.moins(z).div(factor));
  zMax = z.plus(zMax.moins(z).div(factor));
}

void pan(int px, int py, int x, int y) {
  Complexe delta = zPlan(px, py).moins(zPlan(x, y));
  zMin = zMin.plus(delta);
  zMax = zMax.plus(delta);
}

Dans le sketch principal nous ajoutons des fonctions qui écoutent les éventements provenant de la souris et appellent ces méthodes.

Exemple 20.5. Explorer l'ensemble de Mandelbrot

Mandelbrot mandel;
int pressX, pressY;

void setup() {
  size(840, 480);
  mandel = new Mandelbrot(-2.5, -1, 3.5, 100);
  mandel.dessiner();
}

void draw() {
}

void mouseWheel(MouseEvent event) {
  mandel.zoom(mouseX, mouseY, 1 - event.getCount() / 10.0);
  mandel.dessiner();
}

void mousePressed() {
  pressX = mouseX;
  pressY = mouseY;
}

void mouseReleased() {
  mandel.pan(pressX, pressY, mouseX, mouseY);
  mandel.dessiner();
}

Maintenant on peut zoomer / dézoomer avec la molette de la souris et se déplacer en faisant « glisser-déposer ». On laisse au lecteur le soin de décortiquer les détails du code.

Exercice 20.18. Essayez de générer d'autres ensembles en gardant les mêmes règles de jeu mais en changeant légèrement la séquence qu'on itère. Voici quelques pistes :

• z_n+1 = (z_n)³ + c
• z_n+1 = (z̅_n)² + c
• z_n+1 = (|Re(z_n)| + i|Im(z_n)|)² + c

On commence toujours par z₀ = 0. Explorez les ensembles obtenus à partir de ces séquences. Saurez-vous trouver ceci :

Une autre fractale intéressante est l'ensemble de Julia. On définit l'ensemble de Mandelbrot en fixant z₀ à 0 et en faisant varier c. Pour les ensembles de Julia on fait l'inverse, on fixe c et on fait varier z₀.

L'ensemble de Julia associé à une constante c est l'ensemble de nombres complexes z₀ pour lesquels la suite

z_n+1 = (z_n)² + c

ne diverge pas.

Pour dessiner les ensembles de Julia on modifie légèrement notre classe :

• On ajoute des attributs julia et juliaC. Le premier dit si on doit dessiner l'ensemble de Julia ou celui de Mandelbrot. Le second stocke la constante associée à l'ensemble de Julia.

• On fait notre méthode itererSuite() plus générale en ajoutant un argument z0.

• On ajoute une méthode switchJulia() permettant de changer de Mandelbrot à Julia et vice versa.

• On modifie notre méthode dessiner() pour qu'elle dessine l'un des deux ensembles en fonction de la valeur de l'attribut julia.

On modifie le sketch principal pour que l'utilisateur puisse alterner entre les deux ensembles avec le bouton droit de la souris. Lorsque l'on passe de mode « Mandelbrot » en mode « Julia », la constante associe à l'ensemble de Julia correspond à la position de la souris.

Exemple 20.6. Ensembles de Julia

Mandelbrot mandel;
int pressX, pressY;

void setup() {
  size(840, 480);
  mandel = new Mandelbrot(-2.5, -1, 3.5, 100);
  mandel.dessiner();
}

void draw() {
}

void mouseWheel(MouseEvent event) {
  mandel.zoom(mouseX, mouseY, 1 - event.getCount() / 10.0);
  mandel.dessiner();
}

void mousePressed() {
  if (mouseButton == LEFT) {
    pressX = mouseX;
    pressY = mouseY;
  }
}

void mouseReleased() {
  if (mouseButton == LEFT) {
    mandel.pan(pressX, pressY, mouseX, mouseY);
  } else if (mouseButton == RIGHT) {
    mandel.switchJulia(mouseX, mouseY);
  }
  mandel.dessiner();
}

class Mandelbrot {
  // rectangle du plan complexe qu'on veut afficher sur la fenêtre
  Complexe zMin;
  Complexe zMax;

  int nMax; // nombre maximal d'itérations

  boolean julia; // dessinr l'ensemble de Julia ou de Mandelbrot ?
  Complexe juliaC; // la constante c associé à l'ensemble de Julia

  Mandelbrot(float reMin, float imMin, float largeur, int nMax) {
    float reMax = reMin + largeur;
    // Calculé en fonction de la taille de la fenêtre
    // pour garder la même échelle sur les deux axes
    float imMax = imMin + largeur * height / width;

    zMin = new Complexe(reMin, imMin);
    zMax = new Complexe(reMax, imMax);
    this.nMax = nMax;

    julia = false;
    juliaC = new Complexe();
  }

  void dessiner() {
    colorMode(HSB, nMax, 1, 1);
    Complexe zero = new Complexe();
    loadPixels();
    // Pour chaque pixel (x, y) de la fenêtre
    int i = 0;
    for (int y = 0; y < height; y++) {
      for (int x = 0; x < width; x++) {
        // Calculer le nombre complexe correspondant à (x, y)
        Complexe z = zPlan(x, y);
        // Itérer la suite
        int n;
        if (julia) {
          n = itererSuite(z, juliaC);
        } else {
          n = itererSuite(zero, z);
        }
        // Si la suite n'a pas divergé, colorer le pixel en noir
        // sinon, colorer en fonction du nombre d'itérations
        if (n == nMax) {
          pixels[i] = color(0);
        } else {
          pixels[i] = color(n, 1, 1);
        }
        i++;
      }
    }
    updatePixels();
  }

  // Transforme un pixel de la fenêtre en point du plan complexe
  Complexe zPlan(int x, int y) {
    float re = map(x, 0, width, zMin.re(), zMax.re());
    float im = map(y, 0, height, zMin.im(), zMax.im());
    return new Complexe(re, im);
  }

  // Itère la suite z(0) = z0, z(n+1) = f(z(n)) pendant nMax itérations
  // ou jusqu'à ce qu'elle sorte du cercle de rayon 2
  // Renvoie le nombre d'itérations effectuées
  int itererSuite(Complexe z0, Complexe c) {
    Complexe z = z0;
    int n = 0;
    while (n < nMax && z.mod2() < 4) {
      z = f(z, c);
      n++;
    }
    return n;
  }

  Complexe f(Complexe z, Complexe c) {
    return z.mult(z).plus(c);
  }

  void zoom(int x, int y, float factor) {
    Complexe z = zPlan(x, y);
    zMin = z.plus(zMin.moins(z).div(factor));
    zMax = z.plus(zMax.moins(z).div(factor));
  }

  void pan(int px, int py, int x, int y) {
    Complexe delta = zPlan(px, py).moins(zPlan(x, y));
    zMin = zMin.plus(delta);
    zMax = zMax.plus(delta);
  }

  void switchJulia(int x, int y) {
    if (!julia) {
      juliaC = zPlan(x, y);
    }
    julia = !julia;
  }
}

Ensemble de Julia pour c = -0.8400188+0.22430387i

En itérant des suites de nombres complexes, on peut obtenir d'autres formes fractales intéressantes. Considérons la fonction

f(k,z) = a_kz + b_kz̅ + c_k

où k est un entier entre 0 et 3, z est une variable complexe et a_k, b_k et c_k sont des constantes complexes.

On construit une suite z_n ainsi :

• On commence par z₀ = 0
• Pour les membres suivants :
  • On choisit au hasard un k avec probabilité p_k (où p_k sont des constantes, p₀ + p₁ + p₂ + p₃ = 1)
  • z_n+1 = f(k,z_n)

Exercice 20.19. Tracez la la suite z_n en s'inspirant de l'exemple 20.1. Dessinez seulement des points qui correspondent aux membres de la suite, pas les lignes entre les membres consécutifs. Continuez jusqu'à l'infini, ne vous limitez pas à un nombre maximal d'itérations ni dans un cercle. Dessinez le carré du plan complexe délimité par -5.1-0.1i et 5.1+10.1i. Utilisez les coefficients a_k, b_k, c_k et les probabilités p_k suivantes :

Complexe[] a = {
  new Complexe( 0.080,  0.000),
  new Complexe( 0.850, -0.040),
  new Complexe( 0.210,  0.245),
  new Complexe( 0.045, -0.010)
};

Complexe[] b = {
  new Complexe(-0.080,  0.000),
  new Complexe( 0.000,  0.000),
  new Complexe(-0.010, -0.015),
  new Complexe(-0.195,  0.270)
};

Complexe[] c = {
  new Complexe( 0.000,  0.000),
  new Complexe( 0.000,  1.600),
  new Complexe( 0.000,  1.600),
  new Complexe( 0.000,  0.440)
};

float p[] = {0.01, 0.85, 0.07, 0.07};

Si vous avez bien travaillé, votre récompense sera une jolie fougère qui apparaît progressivement devant vos yeux émerveillés.

Exercice 20.20. Jouez avec les différentes constantes magiques. Pouvez-vous générer d'autres formes intéressantes ?

À cause de la grande quantité de calculs nécessaire pour rendre chaque pixel, on sent que notre balade interactive dans les ensembles de Mandelbrot et Julia n'est pas très fluide. Pour quantifier cette notion de fluidité, faisons un « film » qui zoome sur une région intéressante de l'ensemble de Mandelbrot. Au bout de 100 images rendues nous allons nous arrêter et afficher la vitesse de rendu en images par seconde (fps).

Exemple 20.7. Mesure de la vitesse du rendu.

Mandelbrot mandel;

void setup() {
  size(840, 480);
  mandel = new Mandelbrot(-2.5, -1, 3.5, 100);
}

void draw() {
  mandel.zoom(424, 193, 1.1);
  mandel.dessiner();
  if (frameCount == 100) {
    float fps = frameCount * 1000.0 / millis();
    println(fps, "fps");
    noLoop();
  }
}

Même si sur notre système l'utilisation du tableau des pixels à la place de la fonction point() a permis de passer de 0.8 fps à 4 fps (soit une accélération x5), on est très loin des 60 fps auxquelles Processing nous a habitué. Que peut-on faire pour rendre notre code plus rapide ?

Pour mettre en évidence une source potentielle d'amélioration, faisons une petite expérience. Comptons le nombre d'objets de notre classe Complexe crées pour rendre une image. Une façon simple de faire ceci est de déclarer un compteur qu'on va incrémenter dans le constructeur.

class Complexe {
  ...
  Complexe(float re, float im) {
    a = re;
    b = im;
    nbObjets++;
  }
  ...
}

int nbObjets;

void setup() {
  size(840, 480);
  Mandelbrot mandel = new Mandelbrot(-2.5, -1, 3.5, 100);
  nbObjets = 0;
  mandel.dessiner();
  println(nbObjets, "objets crées");
}

On peut constater qu'on a crée plus de 20 millions d'objets (21 502 877 pour être précis) juste pour une image ! Les objets résident dans une zone spéciale de la mémoire appelée tas. La gestion de cette mémoire n'est pas gratuite. À chaque fois qu'on utilise new le système doit trouver un endroit libre de la taille nécessaire pour stocker notre objet. Quand un objet n'est plus utilisé, le système doit recycler la zone qu'il utilise et la rendre disponible pour d'autres objets (et heureusement, parce que sinon avec notre appétit on aurait saturé le tas très rapidement). Tout cela prend un certain temps, certes très petit, mais multiplié par 20 millions, ça commence à peser. Les variables locales habitent dans une autre zone de la mémoire appelée pile. La gestion de la pile ne coûte presque rien car la zone allouée ou libérée se trouve toujours à son « sommet » (d'où le nom).

Développer la classe Complexe nous a beaucoup appris sur la programmation orienté objet. Par la suite cette classe a énormément facilité le développement du code autour des ensembles de Mandelbrot et Julia. Elle nous a permis de passer à un niveau plus haut d'abstraction et de manipuler les nombres complexes sans se soucier des détails de différentes opérations (de la même façon qu'on fait des opérations avec des nombres réels sans devoir manipuler les bits individuels dans leur représentation binaire). Imaginons un instant la vie dure dans un monde sans la classe Complexe. Au lieu de

Complexe z;

on aurait utilisé quelque chose comme

float zRe, zIm;

et la simple instruction

z = z.mult(z).plus(c);

se transformerait en

float t = zRe * zRe - zIm * zIm + cRe;
zIm = 2 * zRe * zIm + cIm;
zRe = t;

Tout cela n'est pas très joli, mais dans certains rares cas où on cherche la performance à tout prix, il faut se passer du confort de l'abstraction et travailler à un niveau plus bas. Nous avons déjà fait cela en renonçant d'utiliser la fonction point() et en mettant les mains dans le cambouis (ou plutôt dans le tableau pixels pas aussi facile à manipuler). Nous allons continuer à contrecœur dans cette lancée et réécrire la classe Mandelbrot sans utiliser la classe Complexe. Cela consiste essentiellement à remplacer chaque variable Complexe par une paire de variables float et à effectuer à la main les opérations arithmétiques. Cette réécriture n'est pas très difficile si on part du code existant, mais développer la classe from skratch sans nombres complexes aurait été une tâche très laborieuse.

Exemple 20.8. Mandelbrot sans nombres complexes. Le code original est laissé en commentaire pour comparaison.

class Mandelbrot {
  // rectangle du plan complexe qu'on veut afficher sur la fenêtre
  //Complexe zMin;
  float zMinRe, zMinIm;
  //Complexe zMax;
  float zMaxRe, zMaxIm;

  int nMax; // nombre maximal d'itérations

  boolean julia; // dessinr l'ensemble de Julia ou de Mandelbrot ?
  //Complexe juliaC; // la constante c associé à l'ensemble de Julia
  float juliaCRe, juliaCIm;

  Mandelbrot(float reMin, float imMin, float largeur, int nMax) {
    float reMax = reMin + largeur;
    // Calculé en fonction de la taille de la fenêtre
    // pour garder la même échelle sur les deux axes
    float imMax = imMin + largeur * height / width;

    //zMin = new Complexe(reMin, imMin);
    zMinRe = reMin;
    zMinIm = imMin;
    //zMax = new Complexe(reMax, imMax);
    zMaxRe = reMax;
    zMaxIm = imMax;

    this.nMax = nMax;

    julia = false;
    //juliaC = new Complexe();
    juliaCRe = 0;
    juliaCIm = 0;
  }

  void dessiner() {
    colorMode(HSB, nMax, 1, 1);
    //Complexe zero = new Complexe();
    loadPixels();
    // On calcule z plus rapidement que avec map()
    float dRe = (zMaxRe - zMinRe) / width;
    float dIm = (zMaxIm - zMinIm) / height;
    // Pour chaque pixel (x, y) de la fenêtre
    int i = 0;
    float zIm = zMinIm;
    for (int y = 0; y < height; y++) {
      float zRe = zMinRe;
      for (int x = 0; x < width; x++) {
        // Calculer le nombre complexe correspondant à (x, y)
        //Complexe z = zPlan(x, y);
        // remplacé par zRe += dRe et zIm += dIm;
        // Itérer la suite
        int n;
        if (julia) {
          //n = itererSuite(z, juliaC);
          n = itererSuite(zRe, zIm, juliaCRe, juliaCIm);
        } else {
          //n = itererSuite(zero, z);
          n = itererSuite(0, 0, zRe, zIm);
        }
        // Si la suite n'a pas divergé, colorer le pixel en noir
        // sinon, colorer en fonction du nombre d'itérations
        if (n == nMax) {
          pixels[i] = color(0);
        } else {
          pixels[i] = color(n, 1, 1);
        }
        i++;
        zRe += dRe;
      }
      zIm += dIm;
    }
    updatePixels();
  }

  // Transforme un pixel de la fenêtre en point du plan complexe
  //Complexe zPlan(int x, int y) {
  //  float re = map(x, 0, width, zMin.re(), zMax.re());
  //  float im = map(y, 0, height, zMin.im(), zMax.im());
  //  return new Complexe(re, im);
  //}
  float zPlanRe(int x) {
    return map(x, 0, width, zMinRe, zMaxRe);
  }

  float zPlanIm(int y) {
    return map(y, 0, height, zMinIm, zMaxIm);
  }

  // Itère la suite z(0) = z0, z(n+1) = f(z(n)) pendant nMax itérations
  // ou jusqu'à ce qu'elle sorte du cercle de rayon 2
  // Renvoie le nombre d'itérations effectuées
  //int itererSuite(Complexe z0, Complexe c) {
  int itererSuite(float z0Re, float z0Im, float cRe, float cIm) {
    //Complexe z = z0;
    float zRe = z0Re;
    float zIm = z0Im;
    int n = 0;
    //while (n < nMax && z.mod2() < 4) {
    while (n < nMax && zRe * zRe + zIm * zIm < 4) {
      //z = z.mult(z).plus(c);
      float t = zRe * zRe - zIm * zIm + cRe;
      zIm = 2 * zRe * zIm + cIm;
      zRe = t;

      n++;
    }
    return n;
  }


  void zoom(int x, int y, float factor) {
    //Complexe z = zPlan(x, y);
    float zRe = zPlanRe(x);
    float zIm = zPlanIm(y);
    //zMin = z.plus(zMin.moins(z).div(factor));
    zMinRe = zRe + (zMinRe - zRe) / factor;
    zMinIm = zIm + (zMinIm - zIm) / factor;
    //zMax = z.plus(zMax.moins(z).div(factor));
    zMaxRe = zRe + (zMaxRe - zRe) / factor;
    zMaxIm = zIm + (zMaxIm - zIm) / factor;
  }

  void pan(int px, int py, int x, int y) {
    //Complexe delta = zPlan(px, py).moins(zPlan(x, y));
    float deltaRe = zPlanRe(px) - zPlanRe(x);
    float deltaIm = zPlanIm(py) - zPlanIm(y);
    //zMin = zMin.plus(delta);
    zMinRe += deltaRe;
    zMinIm += deltaIm;
    //zMax = zMax.plus(delta);
    zMaxRe += deltaRe;
    zMaxIm += deltaIm;
  }

  void switchJulia(int x, int y) {
    if (!julia) {
      //juliaC = zPlan(x, y);
      juliaCRe = zPlanRe(x);
      juliaCIm = zPlanIm(y);
    }
    julia = !julia;
  }
}

On n'a pas besoin de changer le code qui utilise la classe Mandelbrot car nous avons changé son organisation interne mais pas son interface publique. Les méthodes ont gardé les mêmes noms et les mêmes listes d'arguments et donc tout notre code qui utilisait l'ancienne version peut passer à la nouvelle sans aucun changement. C'est encore un argument contre l'utilisation directe des attributs : les attributs ont changé et si notre code les utilisait directement, il aurait fallu le modifier.

Sans la classe Complexe le code est devenu beaucoup moins lisible et beaucoup plus difficile à maintenir, mais au moins est-ce que ça valait la peine ? En mesurant à nouveau la vitesse de rendu d'images nous avons constaté qu'elle est passée de 4 à 9 fps, soit plus que doublée.

Il paraît difficile d'aller beaucoup plus vite car nous avons atteint les limites du processeur (CPU). Même avec la résolution modeste de 840x480, nous avons plus de 400 000 pixels à calculer. Avec une moyenne de 25 itérations par pixel, ça fait plus de 10 millions d'itérations pour rendre une seule image. Et si on pouvait calculer les couleurs de tous les pixels en parallèle ? Les cartes graphiques (GPU) sont conçues pour faire exactement ça. Elles sont équipées de beaucoup d'unités de calcul pas très puissantes, mais leur force est leur grand nombre. Lorsque vous jouez à un jeu vidéo gourmand en ressources, ces petites unités travaillent dur en concert. Elles exécutent toutes le même code, mais avec des données différentes pour calculer la couleur de chaque pixel sur votre écran et vous offrir une expérience fluide plus vraie que nature.

On peut programmer ces unités en utilisant un langage spécial appelé GLSL. Les petits programmes de calcul des couleurs des pixels s'appellent shaders. Les shaders et GLSL dépassent largement la portée de ce cours. Nous dirigeons les lecteurs intéressés vers cet excellent site pour une introduction générale aux shaders et GLSL et ce tutoriel pour leur utilisation en Processing. Ici nous allons simplement montrer l'utilisation d'un shader pour accélérer le rendu des ensembles de Mandelbrot et Julia.

Exemple 20.9. Shader

Notre classe Mandelbrot a un nouvel attribut de type PShader. Dans le constructeur on charge le code du shader à partir d'un fichier. Ce code sera compilé à la volée et exécuté sur le GPU. On utilise la méthode set() pour communiquer au shader quelques paramètres nécessaires pour les calculs qu'il va effectuer : la taille de la fenêtre, le nombre maximal d'itérations, la zone du plan complexe et les paramètres Julia. On utilise cette technique parce que le shader s'exécute sur le GPU et ne peut pas accéder aux variables d'instance correspondantes de notre classe qui sont stockés dans la mémoire du CPU. Le code de la méthode dessiner() a presque disparu. Tout ce qu'on fait dedans est de dessiner un rectangle qui occupe toute la fenêtre et laisser au shader le soin de calculer la couleur de chaque pixel de ce rectangle.

class Mandelbrot {
  // rectangle du plan complexe qu'on veut afficher sur la fenêtre
  Complexe zMin;
  Complexe zMax;

  int nMax; // nombre maximal d'itérations

  boolean julia; // dessinr l'ensemble de Julia ou de Mandelbrot ?
  Complexe juliaC; // la constante c associé à l'ensemble de Julia

  PShader shader; // notre nouveau héro

  Mandelbrot(float reMin, float imMin, float largeur, int nMax) {
    float reMax = reMin + largeur;
    // Calculé en fonction de la taille de la fenêtre
    // pour garder la même échelle sur les deux axes
    float imMax = imMin + largeur * height / width;

    zMin = new Complexe(reMin, imMin);
    zMax = new Complexe(reMax, imMax);
    this.nMax = nMax;

    julia = false;
    juliaC = new Complexe();

    shader = loadShader("mandelbrot.glsl");
    shader.set("size", float(width), float(height));
    shader.set("nMax", nMax);
    setShaderZ();
    setShaderJulia();
  }

  void setShaderZ() {
    shader.set("zMin", zMin.re(), zMin.im());
    shader.set("zMax", zMax.re(), zMax.im());
  }

  void setShaderJulia() {
    shader.set("julia", julia);
    shader.set("juliaC", juliaC.re(), juliaC.im());
  }

  void dessiner() {
    shader(shader);
    noStroke();
    rect(0, 0, width, height);
  }

  // Transforme un pixel de la fenêtre en point du plan complexe
  Complexe zPlan(int x, int y) {
    float re = map(x, 0, width, zMin.re(), zMax.re());
    float im = map(y, 0, height, zMax.im(), zMin.im());
    return new Complexe(re, im);
  }

  void zoom(int x, int y, float factor) {
    Complexe z = zPlan(x, y);
    zMin = z.plus(zMin.moins(z).div(factor));
    zMax = z.plus(zMax.moins(z).div(factor));
    setShaderZ();
  }

  void pan(int px, int py, int x, int y) {
    Complexe delta = zPlan(px, py).moins(zPlan(x, y));
    zMin = zMin.plus(delta);
    zMax = zMax.plus(delta);
    setShaderZ();
  }

  void switchJulia(int x, int y) {
    if (!julia) {
      juliaC = zPlan(x, y);
    }
    julia = !julia;
    setShaderJulia();
  }
}

On retrouve le code disparu de la méthode dessiner() dans le shader. Il est légèrement modifié (après tout, nous avons changé le langage et le paradigme de programmation). Mais vous ne serez pas complètement perdus en le lisant (sauf la fonction hsv2rgb() qui est un sortilège populaire en patois de Shaderland). La double boucle qui parcourt les pixels a disparu parce que les unités de calcul du GPU vont exécuter le code dans main() pour tous les pixels. Le GPU sait faire des calculs rapides avec des vecteurs et matrices de dimension 2, 3 et 4 et on profite de ce support natif pour représenter les nombres complexes par vec2.

#ifdef GL_ES
precision mediump float;
precision mediump int;
#endif

// Variables communes pour tous les pixels
// C'est ici que les valeurs passées par shader.set(...) sont stockées
uniform vec2 size;
uniform int nMax;
uniform vec2 zMin;
uniform vec2 zMax;
uniform bool julia;
uniform vec2 juliaC;

// Conversion HSV -> RGB
// Sur le GPU on représente les couleurs uniquement en RGB
// et il faut convertir manuellement.
vec3 hsv2rgb(vec3 c) {
  vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
  vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
  return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}


int itererSuite(vec2 z0, vec2 c) {
  vec2 z = z0;
  int n = 0;
  while (n < nMax && dot(z, z) < 4.0) {
    z = vec2(z.x * z.x - z.y * z.y, 2 * z.x * z.y) + c;
    n++;
  }
  return n;
}

// Ce programme s'exécute par les unités du GPU
// pour chaque pixel de la fenêtre.
// gl_FragCoord contient les coordonnées du pixel.
// On écrit dans gl_FragColor la couleur calculée
void main() {
  // Calculer le nombre complexe correspondant à (x, y)
  vec2 z = zMin + (zMax - zMin) * gl_FragCoord.xy / size;
  int n;
  // Itérer la suite
  if (julia) {
    n = itererSuite(z, juliaC);
  } else {
    n = itererSuite(vec2(0.0), z);
  }
  // Si la suite n'a pas divergé, colorer le pixel en noir
  // sinon, colorer en fonction du nombre d'itérations
  // La couleur du pixel est un vecteur de 4 éléments (rgba)
  // entre 0 et 1
  if (n == nMax) {
    gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
  } else {
    gl_FragColor = vec4(hsv2rgb(vec3(float(n) / nMax, 1.0, 1.0)), 1.0);
  }
}

Pour tester les limites de notre shader, le programme principal devient beaucoup plus exigeant :

• Nous augmentons la résolution à 960x720, soit presque 700 000 pixels à calculer ;
• Pour un rendu plus précis, nous passons le nombre maximal d'itérations de 100 à 300 ;
• Nous utilisons mouseDragged() pour bouger l'image tout le temps avec la souris et non seulement lorsque l'utilisateur relâche le bouton.

Mandelbrot mandel;

void setup() {
  // On ajoute le paramètre P2D qui permet de profiter d'accélération matérielle
  // On passe à une résolution beaucoup plus élevée
  size(960, 720, P2D);
  // On utilise nMax = 300 pour une rendue plus précise
  // Soyons fous, c'est le GPU qui gère
  mandel = new Mandelbrot(-2.5, -1.5, 4, 300);
  mandel.dessiner();
}

void draw() {
  //mandel.zoom(423, 321, 1.01);
  //mandel.dessiner();
  //if (frameCount == 1000) {
  //  float fps = frameCount * 1000.0 / millis();
  //  println(fps, "fps");
  //  noLoop();
  //}
}

void mouseWheel(MouseEvent event) {
  mandel.zoom(mouseX, mouseY, 1 - event.getCount() / 10.0);
  mandel.dessiner();
}

void mouseDragged() {
  if (mouseButton == LEFT) {
    mandel.pan(pmouseX, pmouseY, mouseX, mouseY);
    mandel.dessiner();
  }
}

void mousePressed() {
  if (mouseButton == RIGHT) {
    mandel.switchJulia(mouseX, mouseY);
  }
  mandel.dessiner();
}

En exécutant ce code, vous pouvez constater que votre balade interactive dans l'ensemble de Mandelbrot devient beaucoup plus fluide et agréable malgré le changement des paramètres qui demande beaucoup plus de calculs. En dé-commentant le code à l'intérieur de draw() vous pouvez mesurer la vitesse du rendu.

On n'a pas besoin d'une carte graphique de harcore gamer ultra-puissante et indécemment chère pour profiter de l'accélération matérielle. La plupart des CPUs modernes ont des unités graphiques intégrées qui qui savent exécuter des shaders avec une vitesse tout à fait respectable. Même les appareils mobiles modernes possèdent des accélérateurs graphiques (d'ailleurs les lignes 1-4 de notre shader servent justement à le rendre compatible avec les appareils mobiles). Sur notre système nous avons mesuré une vitesse d'environ 50 fps en utilisant uniquement le GPU intégré.