编辑代码

(**************************************)
(*
vertex :
boolean : libre ou pas (gelé)
tableau [0..4] : coordonnées ds l'espace
tableau [0..9] : sommets voisins (indice 0 = "puits")
*)
(**************************************)

type vertex = bool * (int array) * (int array)

type tiling = vertex array

let copie (t:tiling) =
  let q = Array.create (Array.length t) (false,[|0;0;0;0;0|],[|0;0;0;0;0;0;0;0;0;0|]) in
    for i=0 to Array.length t -1 do
      let (f,[|c1;c2;c3;c4;c5|],[|s0;s1;s2;s3;s4;s5;s6;s7;s8;s9|]) = t.(i) in q.(i) <- (f,[|c1;c2;c3;c4;c5|],[|s0;s1;s2;s3;s4;s5;s6;s7;s8;s9|])
    done;
    q


open Graphics;;
let ecran = (600,600) and  scale = 11.;;
let dessine = true;;
if dessine then open_graph (" "^(string_of_int (fst ecran))^"x"^(string_of_int (snd ecran)));;
let borne = 6;; (* nb tuiles max pour rechercher l'alternance *)
let thin_color = rgb 255 100 100 and fat_color = white;;
let thin_color = rgb 200 0 0 and fat_color = rgb 255 255 0;;
let thin_color = rgb 103 203 255 and fat_color = rgb 253 128 2;;

let perviy (a,b,c) = a
let vtoroy (a,b,c) = b
let tretiy (a,b,c) = c


(**************************************)
(********** faces d'un pavage *********)
(**************************************)

(* renvoie les listes sans répétition des faces thin et fat (4 sommets) du pavage t *)
let faces (t:tiling) =
  let cmp [a;b;c;d] [e;f;g;h] =
    if a>e then 1 else if a<e then (-1)
    else if b>f then 1 else if b<f then (-1)
    else if c>g then 1 else if c<g then (-1)
    else if d>h then 1 else if d<h then (-1)
    else 0
  in let rec insert e = function
    |[] -> [e]
    |f::r -> match cmp e f with
	|1 -> f::(insert e r)
	|0 -> e::r
	|_ -> e::f::r
  in let thin = ref [] and fat = ref [] in
    for i=1 to Array.length t-1 do
      let (f,c,s) = t.(i) and (x,y,z) = (ref 0, ref 0, ref 0) in
	while s.(!x)=0 do incr x done;
	while !z<10 do
	  y:= (!x+1) mod 10; incr z;
	  while s.(!y)=0 do y:= (!y+1) mod 10; incr z done;
	  if ((tretiy t.(s.(!x))).(!y)<>0)&&((tretiy t.(s.(!y))).(!x)<>0) then
	    begin
	      let [a;b;c;d] = [i; s.(!x); (tretiy t.(s.(!x))).(!y); s.(!y)] in
	      let e = List.hd (List.sort cmp [[a;b;c;d];[b;c;d;a];[c;d;a;b];[d;a;b;c]]) in
		match ((!y+10- !x) mod 10) with
		  |1 -> thin:=insert e !thin
		  |4 -> thin:=insert e !thin
		  |2 -> fat:=insert e !fat
		  |3 -> fat:=insert e !fat
		  |_ -> ()
	    end;
	  x:= !y
	done
    done;
    (!thin,!fat)


(**************************************)
(******** Flips sur les pavages *******)
(**************************************)

(* le i-eme sommet est-il flippable ? *)
let flippable i (t:tiling) =
  let (f,c,s) = t.(i) in f&&(Array.fold_left (fun k v -> k+(if v=0 then 0 else 1)) 0 s=3)

let flippables (t:tiling) =
  let l = ref [] in
    for i=0 to Array.length t-1 do
      if flippable i t then l:=i:: !l
    done;
    !l

(* renvoie l'index (a,b,c) des trois arêtes partant d'un sommet flippable, b etant celle séparant deux tuiles symetriques *)
let ordo s =
  let (a,b,c) = (ref 0,ref 0,ref 0) in
    for k=0 to 9 do
      if s.(k)<>0 then
	if (s.((k+7) mod 10)<>0)&&(s.((k+3) mod 10)<>0) then
	  (a:=(k+3) mod 10; b:=k; c:=(k+7) mod 10)
	else if (s.((k+6) mod 10)<>0)&&(s.((k+4) mod 10)<>0) then 
	  (a:=(k+4) mod 10; b:=k; c:=(k+6) mod 10)	  
    done;
    (!a,!b,!c)

(* fait un flip sur le i-eme sommet *)
let flip i (t:tiling) =
  let (f,c,s) = t.(i) in
  let (e1,e2,e3) = ordo s in
    (* deplace le sommet flippé ds l'espace *)
  let aux e =
    if e mod 2=0 then
      c.(e/2) <- c.(e/2)+1
    else
      c.(((e+5) mod 10)/2) <- c.(((e+5) mod 10)/2)-1
  in aux e1; aux e2; aux e3;
    (* voisinages des trois sommets de l'hexagone reliés à i *)
    let s1 = tretiy (t.(s.(e1)))
    and s2 = tretiy (t.(s.(e2)))
    and s3 = tretiy (t.(s.(e3))) in
      (* création des liens entre i et ses trois nouveaux voisins :*)
      (* d'abord chez ces voisins *)
      (tretiy (t.(s1.(e2)))).(e3) <- i;
      (tretiy (t.(s2.(e3)))).(e1) <- i;
      (tretiy (t.(s3.(e1)))).(e2) <- i;
      (* puis chez i *)
      s.((5+e1) mod 10) <- s2.(e3);
      s.((5+e2) mod 10) <- s3.(e1);
      s.((5+e3) mod 10) <- s2.(e1);
      (* suppression des liens entre i et ses trois anciens voisins : *)
      (* d'abord chez ces voisins *)
      s1.((5+e1) mod 10) <- 0;
      s2.((5+e2) mod 10) <- 0;
      s3.((5+e3) mod 10) <- 0;
      (* puis chez i *)
      s.(e1) <- 0;
      s.(e2) <- 0;
      s.(e3) <- 0


(**************************************)
(*************** dessin ***************)
(**************************************)

let proj [|c1;c2;c3;c4;c5|] =
  let pi5 = 8.*.atan 1./.5. in
  let a = (float_of_int c1)+.(float_of_int c2)*.(cos pi5)+.(float_of_int c3)*.(cos (2.*.pi5))+.(float_of_int c4)*.(cos (3.*.pi5))+.(float_of_int c5)*.(cos (4.*.pi5))
  and b = (float_of_int c2)*.(sin pi5)+.(float_of_int c3)*.(sin (2.*.pi5))+.(float_of_int c4)*.(sin (3.*.pi5))+.(float_of_int c5)*.(sin (4.*.pi5)) in
    (fst ecran/2 + int_of_float (a*.scale), snd ecran/2 + int_of_float (b*.scale))



let draw (t:tiling) = 
  let (thin,fat) = faces t in
  let rec aux coul = function
    |[] -> ()
    |[a;b;c;d]::r -> aux coul r;
	set_color coul;
	fill_poly [|proj (vtoroy t.(a));proj (vtoroy t.(b));proj (vtoroy t.(c));proj (vtoroy t.(d))|];
	set_color black;
	draw_poly [|proj (vtoroy t.(a));proj (vtoroy t.(b));proj (vtoroy t.(c));proj (vtoroy t.(d))|]
  in aux thin_color thin; aux fat_color fat


(* redessine un flip fait *)
let draw_flip i (t:tiling) =
  let draw_vertex j =
    let (f,c,s) = t.(j) in
    let (a,b) = proj c
    and voiz = ref [] in
      for k=9 downto 0 do
	if s.(k)<>0 then
	  let (f',c',s') = t.(s.(k)) in
	    voiz:=(f',proj c',k):: !voiz
      done;
      let (f',(x,y),k) = List.hd !voiz in voiz:= !voiz@[(f',(x,y),k+10)];
	let rec aux = function
	  |(f1,(x1,y1),j1)::(f2,(x2,y2),j2)::r -> aux ((f2,(x2,y2),j2)::r);
	      if f||f1||f2 then (* au moins un des trois sommets pas sur le bord ? *)
		begin
		  if ((j2-j1)=4)||((j2-j1)=1) then set_color thin_color
		  else set_color fat_color;
		  fill_poly [|(a,b);(x1,y1);(x2,y2)|];
		  set_color black;
		  moveto x1 y1; lineto a b; lineto x2 y2
		end
	  |_ -> ()
	in aux (!voiz)
  and s = tretiy t.(i)
  in let (a,b,c) = ordo s
  in draw_vertex i; draw_vertex (s.(a)); draw_vertex (s.(b)); draw_vertex (s.(c))

(*
TODO : représentation "intelligente" :
chaque face appartient à exactement trois ombres 3->2
dans chaque ombre, on peut associer une hauteur à cette face
-> 3 hauteurs, qui peuvent se coder en RGB
*)

let sommets (t:tiling) =
  set_color black;
  for i=0 to Array.length t-1 do
    let (x,y) = proj (vtoroy t.(i)) in
      moveto x y;
      draw_string (string_of_int i)
  done

(**************************************)
(************* sauvegarde *************)
(**************************************)


let to_svg (t:tiling) name dx =
  let (thin,fat) = faces t in
  let o = open_out name in
    Printf.fprintf o "<?xml version=\0341.0\034 standalone=\034no\034?>\n";
    Printf.fprintf o "<!DOCTYPE svg PUBLIC \034-//W3C//DTD SVG 1.1//EN\034 \034http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\034>\n";
    Printf.fprintf o "<svg xmlns=\034http://www.w3.org/2000/svg\034 width=\034%d\034 height=\034%d\034>\n" (fst ecran) (snd ecran);
    let rec aux coul = function
      |[] -> ()
      |[a;b;c;d]::r -> aux coul r;
	  let (xa,ya) = proj (vtoroy t.(a))
	  and (xb,yb) = proj (vtoroy t.(b))
	  and (xc,yc) = proj (vtoroy t.(c))
	  and (xd,yd) = proj (vtoroy t.(d))
	  in Printf.fprintf o "<polygon points=\034%d,%d %d,%d %d,%d %d,%d \034 fill=\034#%s\034 stroke=\034black\034/>\n" (xa+dx) (snd ecran-ya) (xb+dx) (snd ecran-yb) (xc+dx) (snd ecran-yc) (xd+dx) (snd ecran-yd) coul
    in
      aux "6fc9fd" thin; aux "ff700a" fat;
      Printf.fprintf o "</svg>";
      close_out o;
      (List.length thin, List.length fat)


let write_tiling (t:tiling) name =
  let o = open_out name in
    for i=0 to Array.length t-1 do
      let (f,c,s) = t.(i) in
	Printf.fprintf o "(%d,[|" (if f then 1 else 0);
	for j=0 to 4 do
	  Printf.fprintf o "%d;" c.(j);
	done;
	Printf.fprintf o "\008|],[|";
	for j=0 to 9 do
	  Printf.fprintf o "%d;" s.(j);
	done;
	Printf.fprintf o "\008|]);\n";
    done;
    close_out o

let array_of_string s =
  let rec aux i j =
    if i=String.length s then []
    else if (int_of_char s.[j]=45)||((int_of_char s.[j]>=48)&&(int_of_char s.[j]<=57)) then
      aux i (j+1)
    else
      if j>i then
	(int_of_string (String.sub s i (j-i)))::(aux (j+1) (j+1))
      else
	aux (j+1) (j+1)
  in Array.of_list (aux 0 0)

let read_tiling name =
  let t = ref [] in
  let o = open_in name in
  let fin = ref false in
    while !fin=false do
      try
	let l = array_of_string (input_line o) in
	  t:=(l.(0)=1,[|l.(1);l.(2);l.(3);l.(4);l.(5)|],[|l.(6);l.(7);l.(8);l.(9);l.(10);l.(11);l.(12);l.(13);l.(14);l.(15)|]):: !t
      with
	|End_of_file -> fin:=true
    done;
    close_in o;
    Array.of_list (List.rev !t)


(**************************************)
(**** Coupe et projection 5 vers 2 ****)
(**************************************)

(* verifie s'il existe un hypercube pointée ds vect(u,v) qui contient x *)

let permut = [|(0,1,2);(0,1,3);(0,1,4);(0,2,3);(0,2,4);(0,3,4);(1,2,3);(1,2,4);(1,3,4);(2,3,4)|]

(* vecteurs normaux et épaisseur de coupe pour les 10 ombres *)
let normep (u,v) =
  let norm = Array.create 10 (0.,0.,0.)
  and ep = Array.create 10 (0.,0.) in
    for p=0 to 9 do
      let (a,b,c)=permut.(p) in
      let (u1,u2,u3) = (u.(a),u.(b),u.(c))
      and (v1,v2,v3) = (v.(a),v.(b),v.(c)) in
      let (p1,p2,p3) = (u2*.v3-.u3*.v2,u3*.v1-.u1*.v3,u1*.v2-.u2*.v1) in
      let (d1,d2) = (ref 0.,ref 0.) in
	for i=0 to 1 do
	  for j=0 to 1 do
	    for k=0 to 1 do
	      let d = (float_of_int i)*.p1+.(float_of_int j)*.p2+.(float_of_int k)*.p3 in
		if (i,j,k)=(0,0,0) then (d1:=d; d2:=d)
		else (d1:=min !d1 d; d2:=max !d2 d)
	    done
	  done
	done;
	norm.(p) <- (p1,p2,p3);
	ep.(p) <- (!d1,!d2)
    done;
    (norm,ep)

let in_stripe (norm,ep) r x =
  let bol = ref true in
    for i=0 to 9 do
      if !bol then
	let (a,b,c)=permut.(i) in
	let (x1,x2,x3) = (float_of_int x.(a)-.r.(a),float_of_int x.(b)-.r.(b),float_of_int x.(c)-.r.(c))
	and (p1,p2,p3) = norm.(i)
	and (d1,d2) = ep.(i) in
	let s = (x1*.p1+.x2*.p2+.x3*.p3)
	in bol:= (d1<=s)&&(s<d2)
    done;
    !bol
	
(* pts entiers de [-n,n]^5 dans la bande dirigee par (u,v) decalee de r *)
(* utilise la connexite : cherche pts de voisins en voisins *)

let lex [|a;b;c;d;e|] [|a';b';c';d';e'|] =
  if a>a' then 1 else if a<a' then (-1) else
    if b>b' then 1 else if b<b' then (-1) else
      if c>c' then 1 else if c<c' then (-1) else
	if d>d' then 1 else if d<d' then (-1) else
	  if e>e' then 1 else if e<e' then (-1) else 0

exception DejaVu

let rec insert s = function
  |[] -> [s]
  |a::b -> match lex s a with
      |0 -> raise DejaVu
      |1 -> a::(insert s b)
      |(-1) -> s::(a::b)

(* coupe dans l'hypercube de côté n *)
let cut2 (u,v,r) n =
  let np = normep (u,v) in
  (* commence par trouver un point du réseau *)
  let prems = ref false and c=[|0;0;0;0;0|] in
    while not !prems do
      for i=0 to 4 do c.(i) <- (Random.int (2*n+1))-n done;
      prems:= in_stripe np r c
    done;
  (* trouve le reste de proche en proche *)
    let rec aux l ([|a;b;c;d;e|] as s) = 
      if (abs a>n)||(abs b>n)||(abs c>n)||(abs d>n)||(abs e>n)||(not (in_stripe np r s)) then l
      else
	try
	  aux (aux (aux (aux (aux (aux (aux (aux (aux (aux (insert s l) ([|a+1;b;c;d;e|])) ([|a-1;b;c;d;e|])) ([|a;b+1;c;d;e|])) ([|a;b-1;c;d;e|])) ([|a;b;c+1;d;e|])) ([|a;b;c-1;d;e|])) ([|a;b;c;d+1;e|])) ([|a;b;c;d-1;e|])) ([|a;b;c;d;e+1|])) ([|a;b;c;d;e-1|])
	with
	  |DejaVu -> l
    in aux [] c

(* coupe dans l'hypercube de côté n : version itérative *)
let cut2b (u,v,r)  n =
  let np = normep (u,v) and l = ref [] in
    for i=(-n) to n do
      for j=(-n) to n do
	for k=(-n) to n do
	  for l=(-n) to n do
	    for m=(-n) to n do
	      if in_stripe np r [|i;j;k;l;m|] then l:=[|i;j;k;l;m|]:: !l
	    done
	  done
	done
      done
    done;
    !l

(* coupe qui se projette sur l'écran *)
let cut (u,v,r) =
  let np = normep (u,v) in
  (* commence par trouver un point du réseau dans la coupe *)
  let prems = ref false and c=[|0;0;0;0;0|] and k = ref 0 in
    while not !prems do
      for i=0 to 4 do
	c.(i) <- (if !k mod 2=0 then !k/2 else - !k/2)
      done;
      incr k;
      prems:= in_stripe np r c
    done;
  (* puis trouve le reste de proche en proche *)
    let rec aux q ([|a;b;c;d;e|] as s) =
      let (x,y) = proj s in
	if (x<0)||(x>=fst ecran)||(y<0)||(y>=snd ecran)||(not (in_stripe np r s)) then q
	else
	  try
	    aux (aux (aux (aux (aux (aux (aux (aux (aux (aux (insert s q) ([|a+1;b;c;d;e|])) ([|a-1;b;c;d;e|])) ([|a;b+1;c;d;e|])) ([|a;b-1;c;d;e|])) ([|a;b;c+1;d;e|])) ([|a;b;c-1;d;e|])) ([|a;b;c;d+1;e|])) ([|a;b;c;d-1;e|])) ([|a;b;c;d;e+1|])) ([|a;b;c;d;e-1|])
	  with
	    |DejaVu -> q
    in aux [] c

(* determine si deux pts sont voisins, et selon quelle arête 
renvoie k=-2 si non voisin, arête sinon *)
let voisins x y =
  let k = ref (-1) in
  for i=0 to 4 do
    let d = y.(i)-x.(i) in
      if d=1 then
	if !k= (-1) then k:=2*i else k:=(-2)
      else if d= (-1) then
	if !k=(-1) then k:=(2*i+5) mod 10 else k:=(-2)
      else if d<>0 then k:=(-2)
  done;
  !k

(* relie les pts selectionnés par cut *)
(* entree : l = sortie de cut, n = entree de cut, pour marquer pts du bord *)
let link c =
  let c = Array.of_list ([|0;0;0;0;0|]::c) in
  let (t:tiling) = Array.init (Array.length c) (function i -> (true,c.(i),[|0;0;0;0;0;0;0;0;0;0|])) in
    for i=1 to Array.length c-1 do
      for j=1 to Array.length c-1 do
	let e = voisins c.(i) c.(j) in
	  if e>=0 then
	    let s = tretiy t.(i) in
	      s.(e) <- j
      done;
    done;
    t

(* gèle tous les sommets des faces qui ne doivent pas bouger (AC non vérifiable) *)
let freeze (t:tiling) =
  let rec aux d x y j =
    if (tretiy t.(j)).(x)=0 then false (* sorti du pavage *)
    else
      let k = ref ((x+d) mod 10) in
	while (tretiy t.(j)).(!k)=0 do k:= (!k+d) mod 10 done;
	if (d*(10+ !k-x)) mod 10>4 then false
	else if (!k=y)||(!k=(10+2*x-y+5) mod 10) then ((tretiy t.((tretiy t.(j)).(!k))).(x)<>0) (* tuile similaire ou symetrique trouvee *)
	else aux d x y ((tretiy t.(j)).(!k))
  in
    for i=1 to Array.length t-1 do
      let (f,c,s) = t.(i) and (x,y,z,libre) = (ref 0, ref 0, ref 0, ref true) in
	while s.(!x)=0 do incr x done;
	while !z<10 do
	  y:= (!x+1) mod 10; incr z;
	  while s.(!y)=0 do y:= (!y+1) mod 10; incr z done;
	  libre:= !libre && ((10+ !y- !x) mod 10<=4)
	  && (aux 1 !x !y ((tretiy t.(i)).(!y)))
	  && (aux 9 !y !x ((tretiy t.(i)).(!x)))
	  && (aux 9 !x ((!y+5) mod 10) i)
	  && (aux 1 !y ((!x+5) mod 10) i);
	  x:= !y
	done;
	t.(i) <- (!libre,c,s);
    done;
    t

let tiling (u,v,r) = freeze (link (cut (u,v,r)))

let tiling2 (u,v,r) n = freeze (link (cut2 (u,v,r) n))

(* vecteurs de coupe de Penrose *)
let cos2k k = floor (100000.*.cos (8.*.(float_of_int k)*.(atan 1.)/.5.))
let sin2k k = floor (100000.*.sin (8.*.(float_of_int k)*.(atan 1.)/.5.))
(* let u = [|cos2k 0;cos2k 1;cos2k 2;cos2k 3;cos2k 4|]
and v = [|sin2k 0;sin2k 1;sin2k 2;sin2k 3;sin2k 4|]
and r = [|0.2;0.2;0.2;0.2;0.2|] *)
let u = [|100000.; 30901.; -80902.; -80902.; 30901.|]
and v = [|0.; 95105.; 58778.; -58778.; -95105.|]
and r = [|0.2;0.2;0.2;0.2;0.2|]


(**************************************)
(********* Flips probabilisés *********)
(**************************************)

(* mesure la "qualite" du flip en i *)
(* pour ça : verifie l'alternance des tuiles dans la bande ou ca peut changer *)
(* non local a priori...majorer la distance de verif de l'AC ? *)

exception Bord

let qualite i (t:tiling) =
  let (e1,j,e3) = ordo (tretiy t.(i)) in
  let rec aux d k b n = (* verifie alternance des tuiles dans la bande dirigee par j* *)
    let s = tretiy t.(k) and l = ref ((j+d) mod 10) and a = ref 1 in
      while s.(!l)=0 do l:=(!l+d) mod 10; incr a done;
      if (!a=b)||(n=0) then 1 (* tuile identique retrouvee ou rien de trouve assez pres *)
      else if !a>4 then raise Bord (* bord atteint *)
      else if b=0 then aux d s.(!l) (!a) (n-1) (* premiere tuile -> fixe la cible *)
      else let c=[|4;3;2;1|] in
	if !a=c.(b-1) then 0 (* tuile cible trouvee *)
	else aux d s.(!l) b (n-1) (* continue la recherche *)
  in (aux 1 i 0 borne)+(aux 9 i 0 borne)



let mask_init = [|1;1;1|]
let mask_cv = [|0;1;1|]

(* un flip au pif par etape *)
let pif_flips mask (t:tiling) n =
  let k = ref 0 and m = Array.length t-1 in
    while !k<n do
      let i = 1+Random.int m in
	if flippable i t then
	  if mask.(qualite i t)=1 then
	    begin
	      flip i t;
	      incr k;
	      (* if dessine then draw_flip i t *)
	     end
    done;
    t

let mask2_init = [|50;50;50|]
let mask2_cv = [|0;30;95|]

(* a chaque etape, chaque flip a une certaine proba de se faire *)
(* renvoie (a,b,c) : a=bool convergence; b et c = nb flips (neutres/bons) faits *)
let step mask (t:tiling) =
  let (k1,k2) = (ref 0,ref 0) and cv = ref true in
  let rec liste_flips = function
    |(-1) -> []
    |i -> let r = liste_flips (i-1) in
	 if flippable i t then
	   try
	     begin
	       let q = qualite i t in
		 cv:= (!cv)&&(mask.(q)=0); (* y'a-t-il encore moyen de faire des flips ? *)
		 if Random.int 100<mask.(q) then (i,q)::r else r
	     end
	   with Bord -> r
	 else r
  in let rec fait_flips = function
    |[] -> ()
    |(i,q)::r -> fait_flips r;
	if flippable i t then (* flip-test car possibles collisions *)
	  begin
	    flip i t;
	    if dessine then draw_flip i t;
	    if q=1 then incr k1 else incr k2
	  end;
  in fait_flips (liste_flips (Array.length t-1)); (!cv,!k1, !k2)



let n_flips mask (t:tiling) n =
  if dessine then (clear_graph (); draw t);
  let cv = ref false and k = ref 0 and s = ref 0 in
    while (not (!cv))&&((n<0)||(!k<n)) do
      let (a,b,c) = step mask t in
	cv:=a;
	k:= !k+b+c;
	incr s;
    done
(*    Printf.printf "%d étapes, %d flips\n" (!s) (!k)  *)





(*************************************)

let _ =
  let o = open_out "5vers2.log" in
    Random.self_init ();
    Printf.printf "Coupe et projection..."; flush stdout;
    let t = tiling (u,v,[|0.5;0.5;0.5;0.5;0.5|]) in
      Printf.printf "fini.\n";
      let (nthin,nfat) = faces t in
	for rep = 1 to 1 do
(*	open_graph (" "^(string_of_int (fst ecran))^"x"^(string_of_int (snd ecran))); draw t; *)
	  Printf.fprintf o "%d thin rhombs, %d fat rhombs\n" (List.length nthin) (List.length nfat);
	  Printf.fprintf o "Initialisation par 100 millions de flips au pif.\n";
	  Printf.printf "Initialisation : 100 millions de flips au pif..."; flush stdout;
	  n_flips [|50;50;50|] t (100*1000000);
	  Printf.printf "fini.\n";
	  Printf.printf "Convergence par flips probabilisés..."; flush stdout;
	  let cv = ref false and (k1,k2,s) = (ref 0, ref 0, ref 0) in
	    to_svg t (Printf.sprintf "step%04d.svg" !s);
	    (* draw t; *)
	    Printf.fprintf o "étape, flips réversibles, flips irréversibles\n";
	    while not (!cv) do
	      let (a,b,c) = step [|0;30;70|] t in
		cv:=a;
		Printf.fprintf o "%4d, %d, %d\n" !s b c;
		k1:= !k1+b;
		k2:= !k2+c;
		incr s;
		(*draw t;*)
		to_svg t (Printf.sprintf "step%04d.svg" !s)
	    done;
	    Printf.fprintf o "Total : %d étapes, %d flips réversibles, %d flips irréversibles.\n" !s !k1 !k2; flush o;
	    Printf.printf "fini.\n"
      done;
      close_out o;
(*      close_graph *)


(************************************)

let norme2 = Array.fold_left (fun s x -> s+.x*.x) 0.

let f_abs x = if x>=0. then x else -.x

let norminf x = Array.fold_left (fun s x -> max s (f_abs x)) 0. x

let scal x y =
  let s = ref 0. in
    for i=0 to Array.length x-1 do s:= !s+.x.(i)*.y.(i) done;
    !s

let add x y = 
  let z = Array.create (Array.length x) 0. in
    for i=0 to Array.length x-1 do z.(i) <- x.(i)+.y.(i) done;
    z

let mult l x = 
  let z = Array.create (Array.length x) 0. in
    for i=0 to Array.length x-1 do z.(i) <- l*.x.(i) done;
    z

let dist x (u,v,r) =
  let x' = (add (mult (-.1.) r)) (Array.map float_of_int x) in
  let pxu = add (mult (-.(scal x' u)/.(norme2 u)) u) x' in
  let px = add (mult (-.(scal x' v)/.(norme2 v)) v) pxu in
    sqrt (norme2 px)
  
let max_dist (t:tiling) (u,v,r) =
  let k = ref 0. in
    for i=0 to Array.length t-1 do
      if flippable i t then
	k:= max !k (dist (vtoroy t.(i)) (u,v,r))
      done;
      for i=0 to Array.length t-1 do
	if flippable i t then
	  begin
	    let d = dist (vtoroy t.(i)) (u,v,r) in
	    let (x,y) = proj (vtoroy t.(i)) in 
	      set_color blue;
	      if d= !k then fill_circle x y 5;
	      set_color black;
	      moveto x y;
	      draw_string (string_of_int (int_of_float (100.*.d)))
	  end
      done;
      !k