FRACTALS et les ensembles fractals

FRACTALS, comment ça marche ?

Plus généralement...

Retour à la page principale...

FRACTALS, comment ça marche ?

Cher utilisateur, tu as sans doute pu remarquer que FRACTALS, bien qu'il tourne sur le 68030 à 16 Mhz de ton FALCON ne semble pas particulièrement lent... Pourquoi ? Pourtant ce genre de calcul, ça bouffe du CPU : pensez donc, pour chaque point de l'image jusqu'à 32 itérations (de base), une itération consistant en gros à élever un nombre complexe au carré et à vérifier qu'il ne sort pas d'une certaine zone...

Sur que si vous faites ça en C sur des nombres flottants bestialement pour chaque point de l'écran, un Pentium II 266 MHz ramera quand même !

Alors où sont les ruses ? A défaut de vous les dévoiler toutes, je vais vous dévoiler les miennes. Elles sont au nombre de 4:

1. Tous les calculs sont effectués sur des nombres entiers
2. Le DSP est mis à contribution
3. Le 68030 travaille en parallèle avec le DSP
4. Tous les points de l'image ne sont pas bestialement calculés

Mais tout d'abord, les bases:

L'ensemble de MANDELBROT

On va prendre le sacro-saint exemple de l'ensemble de MANDELBROT, bien connu. Sa définition est la suivante:

C'est l'ensemble M(c) des points du plan complexe tels que la suite z_n+1 = z²_n + c reste bornée.

Le premier terme z₀ de la suite étant défini égal à 0 (une autre valeur donne une variante de l'ensemble de Mandelbrot).

On démontre que cette suite devient non bornée dès que l'on trouve n tel que | z_n | > 2, c'est à dire que le point complexe correspondant se trouve hors du disque de rayon 2 et centré sur l'origine. Retenez le bien, ça va nous servir sous la forme:

| z_n |² > 4

Le calcul à effectuer pour chaque itération

En posant z_n = x_n + i *y_n et c = cx + i * cy (décomposition en parties réelle et imaginaire), on obtient:

x_n+1 = x_n² - y_n² + cx et y_n = 2 * x_n * y_n + cy (soit 5 multiplications et 2 additions)

D'autre part, on a besoin de | z_n |² =x_n² + y_n² à comparer à 4 (soit 2 multiplications et 2 additions)

Bref, si on est bestial, pour chaque itération, on a besoin de 7 multiplications et 4 additions.

OK, je te l'accorde, les élévations au carré peuvent très bien être calculées une seule fois, donc 5 multiplications et 4 additions.

Le calcul sur les entiers

Pourquoi sur des entiers ? Parce que ça va bien plus vite ! Oui mais, moi j'ai besoin de nombres "réels" te dis-tu ! Je vais avoir l'air malin avec un "entier" de 3,1415926536 (tiens, il me dit quelque chose, celui-là !).

Bon on va s'arranger...

Si on disait que le nombre "réel" X est représenté par l'entier 10 000 * X par exemple ? Oui mais si X est "grand", mon entier il va se faire exploser ! Certes mais qu'avons nous remarqué dernièrement ? Ah ! Faut suivre : que dès que on obtient un nombre supérieur en valeur absolue à 2, on sait que la suite va alors diverger, on peut donc s'arrêter de coder les "réels" à 2.0.

Bon, soit mais les régles de calcul sont elles les mêmes qu'avec les réels ?

Soient Xe et Ye les représentations entières du réel X et Y, on a:

Xe = 10 000 * X et Ye = 10 000 *Y, d'où:

Xe + Ye = 10 000 * (X + Y) et Xe * Ye = 10 000 * (10 000 * X * Y)

or 10 000 * (X + Y) est la représentation entière du réel X + Y et 10 000 * X * Y celle du réel X * Y

donc il est important de retenir :

1. Pour l'addition, l'addition des représentations entières donne directement la représentation entière de la somme réelle
2. Pour la multiplication, il suffit de diviser le résultat par le coefficient multiplicateur (ici 10 000)

Bref, grosso modo, on s'en sort bien !

Remarques:

1. En fait, on a bêtement changé d'échelle de représentation [-2;2] --> [-20 000; 20 000]
2. Dans le calcul effectif, on a besoin de coder la plage "réelle" [-8.0; 8.0]
3. Les entiers sont codés sur 32 bits et grâce aux multiplications 32x32 --> 64 bits du 68030, on peut avoir la précision voulue
4. En réalité, on ne prend pas 10 000 mais une puissance de 2 (2 ^28 précisément) car la division du résultat (nécessaire pour la multiplication) se résume à un "simple" décalage vers la droite sur 64 bits (bien plus économe en cycles CPU qu'une division)
5. Bien entendu, la limitation sur la précision et la restriction à une échelle [-8;8] ne peut pas être compatible avec tous les algorithmes nécessaires, en particulier, l'algorithme d'estimation de distance ne peut fonctionner que sur des nombres à virgule flottante et donc calculés par le FPU 68882

Utilisation du DSP 56001

Le principe de calcul reste le même sur le DSP 56001 si ce n'est quelques légères différences dues au format fractionnaire et non pas entier de ses registres. La différence principale réside dans le nombre de bits de codage (24 au lieu de 32), donc quand on zoome, le DSP tire sa révérence un peu avant le CPU. Un des gros intérêts du DSP par contre, c'est qu'il peut effectuer un multiplication et une addition en parallèle (1 cycle pour le tout à comparer aux quelques 40 cycles nécessaires pour une multiplication au 68030 ...) d'ou une vitesse de calcul bien plus grande...

Pour minimiser les erreurs de calcul dues à la différence de 8 bits de précision avec le CPU, le DSP recoit du CPU un précalcul des points à évaluer sur une ligne (s'il se contentait d'ajouter un offset à partir d'une valeur initiale, les erreurs s'accumuleraient).

Que fait le 68030 ?

Ben, il calcule lui aussi, tiens !

En gros, ca se passe ainsi:

1. Précalculs d'echelle pour le DSP (effectués par le CPU), puis envoi au DSP avec les coordonnées du point de départ
2. Le CPU demande au DSP de commencer son calcul
3. Le CPU commence aussi son calcul et quand le DSP lui signale qu'il a terminé le sien, le CPU récupère ses résultats, les stocke dans un tableau et fournit au DSP les nouvelles données de calcul
4. Et ainsi de suite jusqu'à ce que le nombre de lignes voulu soit calculé...
5. Ensuite, le CPU convertit toutes ces données (orbites) en une image au format ATARI. Cette partie nécessite également une optimisation conséquente (la fameuse routine indices 1 octet -> format bitplan) qui à elle seule mériterait une page WEB

Ce qu'il est aussi important de noter, c'est que bien que le CPU calcule beaucoup moins vite, il travaille tout de même, ce qui aurait été bien utile si ATARI avait un jour sorti un FALCON 040...

On ne calcule pas TOUS les points

Enfin, la dernière ruse : on évite de calculer tous les points de l'image. On suppose pour cela qu'un point et ses voisins sont suceptibles de renvoyer le même résultat de calcul sous certaines conditions. Vous avez surement remarqué qu'il y a parfois de grands paquets de la même couleur dans les fractals ? Et bien on va prendre comme règle l'affirmation suivante:

Soient trois points A, B et C de l'image tels que A est voisin de B et B voisin de C (on entend par voisin le point immédiatement suivant). Si A et C donnent le même résultat, on dira que B aussi et on évitera de le calculer.

Si on exploite cette hypothèse sur 3 lignes de calcul (tiens ! c'est pour ça qu'une image de FRACTALS a toujours un nombre de lignes multiple de 3 !), cela donne:

On notera C(x, y) le résultat du calcul (orbite) du point A(x,y). Chaque ligne comporte N points

1. Si C( 1, n ) = C( 1, n+2 ), alors C( 1, n+1 ) = C( 1, n ) sinon on calcule C( 1, n+1 ) (Ligne 1)
2. Si C( 3, n ) = C( 3, n+2 ), alors C( 3, n+1 ) = C( 3, n ) sinon on calcule C( 3, n+1 ) (Ligne 3)
3. Si C( 1, n ) = C( 3, n ), alors C( 2, n ) = C( 1, n ) sinon on calcule C( 2, n ) (Colonne 1)
4. Si C( 1, n+1 ) = C( 3, n+1 ), alors C( 2, n+1 ) = C( 1, n+1 ) sinon on calcule C( 2, n+1) (Colonne 2)
5. Si C( 1, n+2 ) = C( 3, n+2 ), alors C( 2, n+2) = C( 1, n+2 ) sinon on calcule C( 2, n+2 ) (Colonne 3)

Et on continue avec Si C(1, n+2) = C(1, n+4) ...

Donc dans le meilleur des cas (toutes les égalités vérifiées pour les points des 3 lignes), on aura calculé seulement :

• N / 2 points sur la ligne n° 1
• AUCUN point sur la ligne n° 2
• N / 2 points sur la ligne n° 3

, soit seulement 1 point sur 3 (!), ce qui est très appréciable d'autant plus que c'est dans les zones "noires" et donc uniformes que le temps de calcul est le plus long.

PLUS GENERALEMENT ...

Tout d'abord qu'est-ce qu'un fractal ? Disons qu'il s'agit d'un terme un peu générique englobant les objets ayant au moins une des caractéristiques défines dans le tableau ci dessous:

UN PEU D'HISTOIRE ...

Depuis le XIX ème siècle, certains mathématiciens comme Cantor, Lorentz, Von Neumann, Barnsley ...et bien d'autres ont commencé à inventer des objets étranges que le monde des Mathématiques n'avait pas pensé à explorer jusque là. Ces objets avaient en commun une définition des plus simples et des propriétés des plus étranges. Contrairement aux belles courbes bien continues et dérivables qui semblaient faire partie la majorité des ensembles, on commençait à découvrir des ensembles complètement 'fous', par exemple :

• Courbe continue en tout point et non dérivable en tout point (fonction de Weirstrass)
• Attracteurs étranges (courbe en ailes de papillon de Lorentz, courbe de Henon...)

Jusque dans les années 70, ces ensembles sont un peu restés des curiosités intellectuelles mais les progrès de l'informatique comme outil de calcul apporta alors une bouffée d'air et on eut alors accès à la représentation graphique de tels ensembles.

UN PEU PLUS SUR L'ENSEMBLE DE MANDELBROT ...

Outre le fait d'être le plus connu des ensembles fractals, cet objet possède des propriétés des plus étonnantes, ce qui le classe encore un peu à part dans la catégorie des fractals.

Gaston JULIA s'interressa à la suite complexe définie par z_n+1 = z²_n + c avec cette fois c fixé (initialisateur) et z₀ = point du plan complexe à tester.Très semblable (à priori) à l'ensemble de MANDELBROT où c est le point du plan à tester et z₀ = 0. On définit ainsi une infinité d'ensembles de JULIA (à chaque fois que l'on choisit un nombre c du plan complexe), alors qu'il y a qu'un ensemble de MANDELBROT.

On appelera par la suite J_c l'ensemble de JULIA d'initialisateur c et M l'ensemble de MANDELBROT.

Pourtant cette apparente ressemblance est fausse :

• M est l'ensemble des points tel que que l'ensemble J_c soit connexe (ie en un seul morceau, toute portion est reliée à une autre sur ce même ensemble)
• M est lui même connexe (ce qui paraît étrange au vu de l'extrême complexité de son contour !)
• M contient TOUS les ensembles de JULIA, si vous zoomez J_c et M au point c, vous obtiendrez à partir d'un certain niveau la même image
• Quelque soit le niveau de zoom sur M, l'image résultante ne sera pas simplement M mais M accompagné d'une multitude de détails toujours différents, l'ensemble de MANDELBROT n'a pas la propriété d'invariance d'échelle mais il est bien plus complexe

Ainsi l'ensemble de MANDELBROT est un peu le père des ensembles de JULIA...

Bibliographie

The science of FRACTAL IMAGES aux éditions Springer-Verlag

En Anglais, mais bourré d'infos, d'algorithmes et de superbes images