prologues:=2;

input geometriesyr16;
input TEX;

vardef DessineObjetColoreUlam(expr col)=
  OrdrepourColore;
  for l=1 upto NF:
    color cc,dd;
    dd=Vision(Gc[l*100+1]);
    cc=Normal(Sommet[Gc[l*100+1]],Sommet[Gc[l*100+2]],Sommet[Gc[l*100+3]]);
    if (ProduitScalaire(dd,cc)<0):
      if pointilles="oui":
	drawoptions(dashed dashpattern(on3bp off9bp));
	trace for k=1 upto Gc[100*l]:
	  Projette(Sommet[Gc[100*l+k]])--
	endfor
	cycle;
      fi;
    else:
      fill for k=1 upto Gc[100*l]:
	Projette(Sommet[Gc[100*l+k]])--
      endfor
      cycle withcolor (ProduitScalaire(dd,cc)/(Module(dd)*Module(cc)))*2*col;%Co[l mod 7];
      trace for k=1 upto Gc[100*l]:
	Projette(Sommet[Gc[100*l+k]])--
      endfor
      cycle withpen pencircle scaled0.5bp withcolor gris;
    fi;
    drawoptions();
  endfor;
enddef;

vardef Icosaedreulam(expr cent)=
  picture ico;
  ico=image(
%%Sommets
    NbS:=12;
    Sommet1:=cent+0.08*(0.8944271,0,0.4472137);
    Sommet2:=cent+0.08*(0.2763932,0.8506507,0.4472137);
    Sommet3:=cent+0.08*(-0.7236067,0.5257311,0.4472137);
    Sommet4:=cent+0.08*(-0.7236067,-0.5257311,0.4472137);
    Sommet5:=cent+0.08*(0.2763932,-0.8506507,0.4472137);
    Sommet6:=cent+0.08*(0,0,1);
    Sommet7:=cent+0.08*(0,0,-1);
    Sommet8:=cent+0.08*(-0.8944271,0,-0.4472137);
    Sommet9:=cent+0.08*(-0.2763932,-0.8506507,-0.4472137);
    Sommet10:=cent+0.08*(0.7236067,-0.5257311,-0.4472137);
    Sommet11:=cent+0.08*(0.7236067,0.5257311,-0.4472137);
    Sommet12:=cent+0.08*(-0.2763932,0.8506507,-0.4472137);
%%Faces
    NF:=20;
    Fc[100]:=3;Fc[101]:=1;Fc[102]:=2;Fc[103]:=6;
    Fc[200]:=3;Fc[201]:=2;Fc[202]:=3;Fc[203]:=6;
    Fc[300]:=3;Fc[301]:=3;Fc[302]:=4;Fc[303]:=6;
    Fc[400]:=3;Fc[401]:=4;Fc[402]:=5;Fc[403]:=6;
    Fc[500]:=3;Fc[501]:=5;Fc[502]:=1;Fc[503]:=6;
    Fc[600]:=3;Fc[601]:=10;Fc[602]:=1;Fc[603]:=5;
    Fc[700]:=3;Fc[701]:=1;Fc[702]:=10;Fc[703]:=11;
    Fc[800]:=3;Fc[801]:=11;Fc[802]:=2;Fc[803]:=1;
    Fc[900]:=3;Fc[901]:=2;Fc[902]:=11;Fc[903]:=12;
    Fc[1000]:=3;Fc[1001]:=12;Fc[1002]:=3;Fc[1003]:=2;
    Fc[1100]:=3;Fc[1101]:=3;Fc[1102]:=12;Fc[1103]:=8;
    Fc[1200]:=3;Fc[1201]:=3;Fc[1202]:=8;Fc[1203]:=4;
    Fc[1300]:=3;Fc[1301]:=4;Fc[1302]:=8;Fc[1303]:=9;
    Fc[1400]:=3;Fc[1401]:=4;Fc[1402]:=9;Fc[1403]:=5;
    Fc[1500]:=3;Fc[1501]:=5;Fc[1502]:=9;Fc[1503]:=10;
    Fc[1600]:=3;Fc[1601]:=7;Fc[1602]:=8;Fc[1603]:=12;
    Fc[1700]:=3;Fc[1701]:=7;Fc[1702]:=9;Fc[1703]:=8;
    Fc[1800]:=3;Fc[1801]:=7;Fc[1802]:=10;Fc[1803]:=9;
    Fc[1900]:=3;Fc[1901]:=7;Fc[1902]:=11;Fc[1903]:=10;
    Fc[2000]:=3;Fc[2001]:=7;Fc[2002]:=12;Fc[2003]:=11;
    DessineObjetColoreUlam(vert);
    );
  ico
enddef;

vardef Cubeulam(expr cent)=
  picture cub;
  cub=image(
    NbS:=8;
    Sommet1:=cent+(-0.18/4,-0.18/4,-0.18/4)+(0.18/2,0,0);
    Sommet2:=cent+(-0.18/4,-0.18/4,-0.18/4)+(0.18/2,0.18/2,0);
    Sommet3:=cent+(-0.18/4,-0.18/4,-0.18/4)+(0,0.18/2,0);
    Sommet4:=cent+(-0.18/4,-0.18/4,-0.18/4)+(0,0,0);
    Sommet5:=cent+(-0.18/4,-0.18/4,-0.18/4)+(0,0,0.18/2);
    Sommet6:=cent+(-0.18/4,-0.18/4,-0.18/4)+(0.18/2,0,0.18/2);
    Sommet7:=cent+(-0.18/4,-0.18/4,-0.18/4)+(0.18/2,0.18/2,0.18/2);
    Sommet8:=cent+(-0.18/4,-0.18/4,-0.18/4)+(0,0.18/2,0.18/2);
%%Faces
    NF:=6;
    Fc[100]:=4;Fc[101]:=1;Fc[102]:=4;Fc[103]:=3;Fc[104]:=2;
    Fc[200]:=4;Fc[201]:=4;Fc[202]:=5;Fc[203]:=8;Fc[204]:=3;
    Fc[300]:=4;Fc[301]:=1;Fc[302]:=6;Fc[303]:=5;Fc[304]:=4;
    Fc[400]:=4;Fc[401]:=5;Fc[402]:=6;Fc[403]:=7;Fc[404]:=8;
    Fc[500]:=4;Fc[501]:=2;Fc[502]:=3;Fc[503]:=8;Fc[504]:=7;
    Fc[600]:=4;Fc[601]:=1;Fc[602]:=2;Fc[603]:=7;Fc[604]:=6;
    DessineObjetColoreUlam(jaune);
    );
  cub
enddef;

vardef pgcd(expr A,B) = 
    save a,b,r;
    numeric a,b,r;
    a := A; b := B;
    forever: 
        r := a mod b;
        a := b; b := r;
        exitunless r > 0;
    endfor;
    a
enddef;

vardef prem(expr P)=
  boolean reponse;
  reponse=false;
  numeric Div,DIV;%DIV pour compter le nombre de diviseurs
  Div=0;DIV=0;
  if P=1:
    Div:=1;
  else:
    for g=2 upto P-1:
      if (pgcd(P,g)<>1):
	Div:=Div+1;
	if (P mod g)=0:
	  DIV:=DIV+1;
	fi;
      fi;
    endfor;
  fi;
  if Div=0:
    reponse:=true;
  fi;
  reponse
enddef;

vardef ulam(expr nb,nbdiag,arret)=
  save spiulam;
  picture spiulam;
  spiulam=image(
    if prem(nb)=true:
      trace Cubeulam((0,0,0));
      label(TEX(""&decimal(nb)&""),Projette((0,0,0)));
    else:
      trace Icosaedreulam((0,0,0));
      label(TEX(""&decimal(nb)&""),Projette((0,0,0))) withcolor 0.9green;
    fi;
    color ptd;
    ptd=(0,0,0);
    coef:=0.35;
    k:=nb;
    for l=1 upto (nbdiag*2+1):
      if (l mod 2)=1:
	for j=1 upto l:
	  if ((k+1)<(nb+arret)):
	  k:=k+1;
	  if prem(k)=true:
	    trace Cubeulam(ptd+coef*(0,j,0));
	    label(TEX(""&decimal(k)&""),Projette(ptd+coef*(0,j,0)));
	  else:
	    trace Icosaedreulam(ptd+coef*(0,j,0));
	    label(TEX(""&decimal(k)&"$_{"&decimal(DIV)&"}$"),Projette(ptd+coef*(0,j,0)));
	  fi;
	  fi;
	endfor;
	ptd:=ptd+coef*(0,l,0);
	for j=1 upto l:
	  if ((k+1)<(nb+arret)):
	  k:=k+1;
	  if prem(k)=true:
	    trace Cubeulam(ptd+coef*(j,0,0));
	    label(TEX(""&decimal(k)&""),Projette(ptd+coef*(j,0,0)));
	  else:
	    trace Icosaedreulam(ptd+coef*(j,0,0));
	    label(TEX(""&decimal(k)&"$_{"&decimal(DIV)&"}$"),Projette(ptd+coef*(j,0,0)));
	  fi;
	  fi;
	endfor;
	ptd:=ptd+coef*(l,0,0);
      fi;
      if (l mod 2)=0:
	for j=1 upto l:
	  if((k+1)<(nb+arret)):
	  k:=k+1;
	  if prem(k)=true:
	    trace Cubeulam(ptd+coef*(0,-j,0));
	    label(TEX(""&decimal(k)&""),Projette(ptd+coef*(0,-j,0)));
	  else:
	    trace Icosaedreulam(ptd+coef*(0,-j,0));
	    label(TEX(""&decimal(k)&"$_{"&decimal(DIV)&"}$"),Projette(ptd+coef*(0,-j,0)));
	  fi;
	  fi;
	endfor;
	ptd:=ptd+coef*(0,-l,0);
	for j=1 upto l:
	  if((k+1)<(nb+arret)):
	  k:=k+1;
	  if prem(k)=true:
	    trace Cubeulam(ptd+coef*(-j,0,0));
	    label(TEX(""&decimal(k)&""),Projette(ptd+coef*(-j,0,0)));
	  else:
	    trace Icosaedreulam(ptd+coef*(-j,0,0));
	    label(TEX(""&decimal(k)&"$_{"&decimal(DIV)&"}$"),Projette(ptd+coef*(-j,0,0)));
	  fi;
	  fi;
	endfor;
	ptd:=ptd+coef*(-l,0,0);
      fi;
    endfor;
    );
  spiulam
enddef;

figureespace(-7.5u,-5u,6.5u,5u);
trace feuillet withcolor blanc;
Initialisation(5,30,20,650);
pointilles:="non";
drawoptions(withcolor orange);
for t=-3 upto 3:
  trace segment((t*0.35,-3*0.35,0),(t*0.35,3*0.35,0));
  trace segment((-3*0.35,t*0.35,0),(3*0.35,t*0.35,0));
endfor;
drawoptions();
trace ulam(41,4,49);
finespace;
end