prologues:=2;

input geometriesyr16;
input TEX;

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 cercles((0,0,0),(0.1,0,0),(0,0,0),(1,0,0),(0,1,0)) dashed evenly;
      label(TEX(""&decimal(nb)&""),Projette((0,0,0)));
    else:
      label(TEX(""&decimal(nb)&""),Projette((0,0,0))) withcolor 0.9green;
    fi;
    color ptd;
    ptd=(0,0,0);
    coef:=0.25;
    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 cercles(ptd+coef*(0,j,0),ptd+coef*(0,j,0)+(0.1,0,0),(0,0,0),(1,0,0),(0,1,0)) dashed evenly;% withcolor 0.9green;
	    label(TEX(""&decimal(k)&""),Projette(ptd+coef*(0,j,0)));
	  else:
	    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 cercles(ptd+coef*(j,0,0),ptd+coef*(j,0,0)+(0.1,0,0),(0,0,0),(1,0,0),(0,1,0)) dashed evenly;%withcolor 0.9green;
	    label(TEX(""&decimal(k)&""),Projette(ptd+coef*(j,0,0)));
	  else:
	    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 cercles(ptd+coef*(0,-j,0),ptd+coef*(0,-j,0)+(0.1,0,0),(0,0,0),(1,0,0),(0,1,0)) dashed evenly;%withcolor 0.9green;
	    label(TEX(""&decimal(k)&""),Projette(ptd+coef*(0,-j,0)));
	  else:
	    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 cercles(ptd+coef*(-j,0,0),ptd+coef*(-j,0,0)+(0.1,0,0),(0,0,0),(1,0,0),(0,1,0)) dashed evenly;%withcolor 0.9green;
	    label(TEX(""&decimal(k)&""),Projette(ptd+coef*(-j,0,0)));
	  else:
	    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(-7u,-5u,7u,5u);
trace feuillet withcolor blanc;
Initialisation(5,30,20,800);
trace ulam(41,4,49);
finespace;

end