input geometriesyr16;

%input marith;% à utiliser pour de grandes valeurs des côtés !!!

vardef pgcd(expr a,b)=
  save $;
  numeric $;
  boolean reste;
  reste=true;
  if a>b:
    dividende0=a;
    diviseur0=b;
  else:
    dividende0=b;
    diviseur0=a;
  fi;
  re0=diviseur0;
  j=0;
  forever: exitunless reste;
    if (dividende[j] mod diviseur[j])=0:
      reste:=false;
    else:
      j:=j+1;
      re[j]=(dividende[j-1] mod diviseur[j-1]);
      dividende[j]=diviseur[j-1];
      diviseur[j]=re[j];
    fi;
  endfor;
  $=re[j];
  $
enddef;

vardef ppcm(expr a,b)=
  save $;
  numeric $;
  $=a*b/pgcd(a,b);
  $
enddef;

prologues:=2;

vardef roulades(expr num,deno,total)=
  save $;
  picture $;
  path cc,poly;
  %définition des constantes de bases
  cote:=3;
  totalcen:=0;
  totaltour:=num-1;
%définition des points de bases
  pair O,A[],B[],C,D;
  O=(0,0);
  B2=u*(-5,-5);
  B0=u*(5,-5);
  B1=rotation(B0,B2,60);
  B3=B0;
  A0=B0;
  A1=(num/deno)[B0,B1];
  A2=rotation(A1,A0,-60);
  A3=A0;
  %définition des centres de rotation
  marque_p:="plein";
  pair cen[];
  cen[0]=B0;
  for k=0 upto 2:
    for j=1 upto deno:
      cen[deno*k+j]:=(j/deno)[B[k],B[k+1]];
      pointe(cen[deno*k+j]);
      totalcen:=totalcen+1;
    endfor;
  endfor;
  totalcentre:=totalcen;
  for p=1 upto totaltour:
    for k=0 upto totalcen:
      cen[p*totalcen+k]=cen[k];
      totalcentre:=totalcentre+1;
    endfor;
  endfor;
%définition des polygones
  poly=A0--A1--A2--cycle;
  cc=B0--B1--B2--cycle;
%traçage
  trace cc withcolor bleu;
  trace poly withcolor violet;
  pair I[],J[];
  I0=A0;
  J0=iso(A0,A1,A2);
  w:=0;
  for k=num step num until total:
    if 7*((k div 7)+1)-k<num:
      if (k mod 3)=0:
	for j=1 upto 120:
	  I[120*w+j]=rotation(I[120*w+j-1],cen[7*((k div 7)+1)],1);
	endfor;
	w:=w+1;
      else:
	for j=1 upto 120:
	  I[120*w+j]=rotation(I[120*w+j-1],cen[k],1);
	endfor;
	for j=1 upto 120:
	  I[120*w+120+j]=rotation(I[120*w+120+j-1],cen[7*((k div 7)+1)],1);
	endfor;
	w:=w+2
      fi;
    elseif (k mod 7)=0:
      for j=1 upto 240:
	I[120*w+j]=rotation(I[120*w+j-1],cen[k],1);
      endfor;
      w:=w+2;
    else:
      for j=1 upto 120:
	I[120*w+j]=rotation(I[120*w+j-1],cen[k],1);
      endfor;
      w:=w+1;
    fi;
  endfor;
  path cd;
  cd=I0
  for j=2 step 1 until (120*w):
    ..I[j]
  endfor;
  $=image(
    trace cd withcolor orange;
    );
  $
enddef;
figure(-12u,-12u,12u,12u);
trace roulades(2,7,4);
fin;
figure(-12u,-12u,12u,12u);
trace roulades(2,7,6);
fin;
figure(-12u,-12u,12u,12u);
trace roulades(2,7,8);
fin;
end