Triangle de Sierpinski

[Billy the Kid]

Bonjour, je débute avec xcas et je voudrais programmer la construction du triangle de Sierpinski.
Voilà le principe, je suis parti d'un triangle équilatéral ABC, je place un point M en son centre puis j'effectue un tirage aléatoire qui donne un entier entre 1 et 3. En fonction du résultat du tirage, je place le point N1 au milieu de [AM], ou au milieu de [BM] ou au milieu de [CM] et j'effectue un autre tirage aléatoire qui donne naissance au milieu de [AN1], ou au milieu de [BN1] ou au milieu de [CN1], et ainsi de suite...
J'ai galéré comme un débutant pour pouvoir afficher les points dans le triangle mais j'ai trouvé quelque chose quand même.
Mais bon, voilà mon problème : quand je calcule les coordonnées des points N avec 100 itérations, cela fonctionne mais quand je veux faire 1000 itérations, xcas ne peut effectuer les calculs. Je voudrais afficher au moins 10000 points pour obtenir la figure mais aurais-je déjà dépassé les capacités de calculs de xcas? Bizarre, pourquoi cela plante-t-il et quelqu'un pourrait-il m'expliquer comment obtenir ces 10000 points?
Merci à vous.
Voici mon code :

Code : Tout sélectionner

Sier():={
local A,B,C,T,M,i,N,x,y,Tirage;
A:=point(-1/2,-sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
B:=point(1/2,-sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
C:=point(0,sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
T:=triangle(A,B,C,affichage=rouge);
M:=point(0,0,affichage=bleu)
x:=0;
y:=0;
N:=NULL;
pour i de 1 jusque 100 faire
Tirage:=1+alea(3);
si Tirage==1 alors
  x:=(-1/2+x)/2;
  y:=(-sqrt(3)/4+y)/2;
sinon
  si Tirage==2 alors
    x:=(1/2+x)/2;
    y:=(-sqrt(3)/4+y)/2;
  sinon
    x:=x/2;
    y:=(sqrt(3)/4+y)/2;
  fsi;
fsi;
N:=N,point(x,y,affichage=point_point)
fpour
retourne(T,A,B,C,M,N);
}
:;

J'ai besoin de votre aide.

[parisse]

Il faut faire les calculs en approché. Pour vous en convaincre faites c:=coordonnees(Sier()), puis par exemple c[98]. Il suffit de modifier l'affectation de x et y avec un evalf, par exemple
x:=evalf((-1/2+x)/2);

[parisse]

On peut aussi accélérer en utilisant l'affectation en place (sinon le remplissage de N est quadratique en fonction du nombre de points).

Code : Tout sélectionner

Sier():={
local A,B,C,T,M,i,N,x,y,Tirage;
A:=point(-1/2,-sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
B:=point(1/2,-sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
C:=point(0,sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
T:=triangle(A,B,C,affichage=rouge);
M:=point(0,0,affichage=bleu)
x:=0;
y:=0;
N:=[0$10000];
pour i de 1 jusque 10000 faire
Tirage:=1+alea(3);
si Tirage==1 alors
  x:=(-1/2+x)/2.;
  y:=(-sqrt(3)/4+y)/2.;
sinon
  si Tirage==2 alors
    x:=(1/2+x)/2.;
    y:=(-sqrt(3)/4+y)/2.;
  sinon
    x:=x/2.;
    y:=(sqrt(3)/4+y)/2.;
  fsi;
fsi;
N[i]=<point(x,y,affichage=point_point)
fpour
retourne(T,A,B,C,M,N);
}
:;

[alb]

Pour info, une version courte mais moins rapide, inspiré d'un exemple de guillaume connan:

Code : Tout sélectionner

Sierpinski(n):={
local T,G,L,k;
T:=[[-0.5,0],[0.5,0],[0,0.5*sqrt(3)]];
G:=(1/3)*(T[0]+T[1]+T[2]);
L:=G;
pour k de 1 jusque n faire
  G:=0.5*(G+T[hasard(3)]);
  L:=L,G;
fpour;
affichage(point_point);
nuage_points([L]);
}:;

[Billy the Kid]

Merci de prêter attention à mes soucis.
Je viens de suivre les conseils de B. Parisse et j'arrive à faire fonctionner ce code avec 7500 itérations. Mais cela ne passe pas pour 8000 itérations.
Ai-je atteint les limites de xcas? Avez-vous les mêmes problème que moi en compilant ce code?

Code : Tout sélectionner

Sier():={
local A,B,C,T,M,i,N,x,y,Tirage;
A:=point(-1/2,-sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
B:=point(1/2,-sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
C:=point(0,sqrt(3)/4,affichage=epaisseur_point_2+point_point+rouge);
T:=triangle(A,B,C,affichage=rouge);
M:=point(0,0,affichage=bleu)
x:=0;
y:=0;
N:=[0$7500];
pour i de 1 jusque 7500 faire
Tirage:=1+alea(3);
si Tirage==1 alors
  x:=evalf((-1/2+x)/2);
  y:=evalf((-sqrt(3)/4+y)/2);
sinon
  si Tirage==2 alors
    x:=evalf((1/2+x)/2);
    y:=evalf((-sqrt(3)/4+y)/2);
  sinon
    x:=evalf(x/2);
    y:=evalf((sqrt(3)/4+y)/2);
  fsi;
fsi;
N[i]=<point(x,y,affichage=point_point+epaisseur_point_2);
fpour
retourne(T,A,B,C,M,N);
}
:;

MErci àvous

[alb]

Chez moi Sier() passe pour 17000 mais plante pour 20000.
le Sierpinski(20000) précédent aboutit en 11 secondes.
Le plus joli et le plus rapide est obtenu avec:

Code : Tout sélectionner

Sierpinski(n):={
local T,G,L,k;
T:=[[-0.5,0],[0.5,0],[0,0.5*sqrt(3)]];
G:=(1/3)*(T[0]+T[1]+T[2]);
L:=[0$n];
L[0]=<G;
pour k de 1 jusque n faire
  G:=0.5*(G+T[hasard(3)]);
  L[k]=<G;
fpour;
affichage(point_point);
nuage_points([L]);
}:;

On a Sierpinski(100000) en moins de 2 secondes.
et Sierpinski(1000000) en 20 secondes.

[parisse]

Il y a en effet une reservation memoire trop grande que je corrige et qui fait planter Sier() pour de grandes valeurs, je la corrige.

[Billy the Kid]

Il y a en effet une reservation memoire trop grande que je corrige et qui fait planter Sier() pour de grandes valeurs, je la corrige.

Merci d'avoir prêté attention à mon problème. J'attends de vos nouvelles.

Chez moi Sier() passe pour 17000 mais plante pour 20000.
le Sierpinski(20000) précédent aboutit en 11 secondes.

C'est curieux car pour moi Sierpinski(8000) fonctionne mais Sierpinski(9000) plante. Comment expliquer ces différences avec Sierpinki(20000) : j'en suis loin.....
Je suis avec xcas 0.8.5 (c) 2000-8 sous Windows7(64bits) Processeur Intel(R)Core(TM)i5 CPU 3.20 GHz. Comment cela se fait-il qu'il y ait une différence si grande avec vos essais ?

[parisse]

Ca vient surement de la réservation mémoire maximale autorisée sous win et sous unix qui sont dépendantes de l'implémentation de la classe vector de la STL.

[Billy the Kid]

...dépendantes de l'implémentation de la classe vector de la STL.

Je penchais bien pour une différence entre les deux OS mais là, cela dépasse mes compétences et je ne comprends malheureusement pas ce que cela veut dire.
Merci quand même pour la réponse.