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,pas)=
  save $;
  picture $;
  $=image(
    if prem(nb)=true:
      fill fullcircle scaled 7mm withcolor 0.9green;
      label(TEX(""&decimal(nb)&""),(0,0));
    else:
      label(TEX(""&decimal(nb)&""),(0,0)) withcolor 0.9green;
    fi;
    pair ptd;
    ptd=(0,0);
    u:=8mm;
    k:=nb;
    for l=1 upto pas:
      if (l mod 2)=1:
	for j=1 upto l:
	  k:=k+1;
	  if prem(k)=true:
	    fill fullcircle scaled 7mm shifted (ptd+u*(0,j)) withcolor 0.9white;
	    label(TEX(""&decimal(k)&""),ptd+u*(0,j));
	  else:
	    label(TEX(""&decimal(k)&"$_{"&decimal(DIV)&"}$"),ptd+u*(0,j));
	  fi;
	endfor;
	ptd:=ptd+u*(0,l);
	for j=1 upto l:
	  k:=k+1;
	  if prem(k)=true:
	    fill fullcircle scaled 7mm shifted (ptd+u*(j,0)) withcolor 0.9white;
	    label(TEX(""&decimal(k)&""),ptd+u*(j,0));
	  else:
	    label(TEX(""&decimal(k)&"$_{"&decimal(DIV)&"}$"),ptd+u*(j,0));
	  fi;
	endfor;
	ptd:=ptd+u*(l,0);
      fi;
      if (l mod 2)=0:
	for j=1 upto l:
	  k:=k+1;
	  if prem(k)=true:
	    fill fullcircle scaled 7mm shifted (ptd+u*(0,-j)) withcolor 0.9white;
	    label(TEX(""&decimal(k)&""),ptd+u*(0,-j));
	  else:
	    label(TEX(""&decimal(k)&"$_{"&decimal(DIV)&"}$"),ptd+u*(0,-j));
	  fi;
	endfor;
	ptd:=ptd+u*(0,-l);
	for j=1 upto l:
	  k:=k+1;
	  if prem(k)=true:
	    fill fullcircle scaled 7mm shifted (ptd+u*(-j,0)) withcolor 0.9white;
	    label(TEX(""&decimal(k)&""),ptd+u*(-j,0));
	  else:
	    label(TEX(""&decimal(k)&"$_{"&decimal(DIV)&"}$"),ptd+u*(-j,0));
	  fi;
	endfor;
	ptd:=ptd+u*(-l,0);
      fi;
    endfor;
    );
  $
enddef;

input TEX;


beginfig(1);
  draw ulam(1,6);
endfig;


beginfig(2);
  draw ulam(41,15);
endfig;
end