Octonacci[k_]:=Module[{octo,i},
octo=Table[0,{i,1,k}];
octo[[1]]=1;
octo[[2]]=2;
Do[ octo[[i]]=octo[[i-2]]+2*octo[[i-1]],{i,3,k}];
octo];

Octonacci[10]

modulomaille[coord_]:=Module[{cmaille,i},
cmaille={0,0};
Do[
{s=coord[[i]];
If[ s ≥0, cmaille[[i]]=N[FractionalPart[s]], cmaille[[i]]=N[FractionalPart[s]+1]];
},{i,1,2}]
;cmaille];

(***************************************************)
(* recherche des voisins dans la zone d'accpetance *)
(***************************************************)
voisins[i_,j_,ij_,etoile_,tt_,cote1_]:=Module[{existe,iix,iiy},
existe=False;
iix=i+etoile[[1,ij]];
iiy=j+etoile[[2,ij]];
If [(iiy ≥ 1) ∧ (iiy ≤ cote1)∧ ( iix≥tt[[iiy,1]]) ∧( iix ≤tt[[iiy,2]]),existe=True];
{existe,iix,iiy}
];

(****************************)
(*** Programme principal ****)
(***************************)

k=8; (* ordre de l'approximant*)

(******************************)
traceperp=False ;(* decide si l'on trace la selection dans E_perp *)
tracepar=False  ;(* decide si l'on montre la maille étendue dans E_par *)
traceparmod=False  ;(* decide si l'on montre la maille réduite dans E_par *)
traceparmodcom=False; (* decide si l'on trace la maille réduite commensurable dans E_par *)

startcpu=TimeUsed[];
startreel=SessionTime[];

rayonperp={0.5,0.1,0.1,0.05,0.02,0.01};
rayonpar={0.5,0.1,0.1,0.05,0.02,001};
rayonparmod={0.3,0.02,0.01,0.005,0.001,0.005};
(************************************)
(***** quelques définitions *********)
(************************************)
kpair=EvenQ[k];
octo=Octonacci[k+3];

gama=(({{0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1}, {0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1}, {0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1}, {0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1}})); (* les 16 sommets de l'hypercube *)
M=(({{1, 0, 1, -1}, {0, 1, 1, 1}, {1, 1, 1, 0}, {-1, 1, 0, 1}})); (* la matrice d'inflation *)
Print["M  : ",MatrixForm[M]];

x=N[1/Sqrt[2]];
f1=1+Sqrt[2];
f1k=f1^k;
f2=1-Sqrt[2];
f2k=f2^k;
pipar4=(({{1, 0, x, -x}, {0, 1, x, x}, {x, x, 1, 0}, {-x, x, 0, 1}}))/2;
pipar2=(({{1, 0, x, -x}, {0, 1, x, x}})); (*projecteur réduit à une matrice 2X4 pour engendrer les deux coordonnées finales *)
piperp4=(({{1, 0, -x, x}, {0, 1, -x, -x}, {-x, -x, 1, 0}, {x, -x, 0, 1}}))/2;
piperp2=(({{-x, -x, 1, 0}, {x, -x, 0, 1}}));


Mk=MatrixPower[M,k];
Print[" M^",k, "  = ", MatrixForm[Mk]];
Mkinv=MatrixPower[M,-k];
Print[" M^",-k, "  = ", MatrixForm[Mkinv]];

e=IdentityMatrix[4];
generperp={piperp2.e[[3]]*f2^k,piperp2.e[[4]]*f2^k};
alpha:=-octo[[k]];
If [EvenQ[k],alpha=-octo[[k]],alpha:=octo[[k]]];
If [EvenQ[k],delta=octo[[k+1]],delta=-octo[[k+1]]];
beta=alpha+delta;
cote=Abs[beta-2*alpha]+1;
cote1=cote-1;
Print[" alpha : ",alpha,"   beta : ", beta, "   delat : ", delta, "  cote : ",cote];
etoile=(({{alpha, -alpha, alpha, -alpha, beta, -beta, 0, 0}, {-alpha, alpha, alpha, -alpha, 0, 0, beta, -beta}}));



(*****************************************************************************************************)
(******** construction du tableau contenant les points sélectionnés dans la zone d'acceptance ********)
(*****************************************************************************************************)

tt=Table[0,{i,1,cote},{j,1,3}];
ix=Abs[alpha];iy=Abs[alpha];
(*bas du tableau*)
Do[{tt[[j,1]]=1;tt[[j,2]]=cote-1},{j,1,cote}];
Do[{tt[[j,1]]=ix+2-j;tt[[j,2]]=cote-ix+j-2},{j,1,iy}];

(*haut du tableau*)
Do[{tt[[cote-j,1]]=ix+1-j;tt[[cote-j,2]]=cote-ix+j-1},{j,1,iy-1}];
tt[[1,3]]=1;
Do[tt[[j,3]]=tt[[j-1,3]]+tt[[j-1,2]]-tt[[j-1,1]]+1,{j,2,cote-1} ];nsite=tt[[cote-1,3]]+tt[[cote-1,2]]-tt[[cote-1,1]];
Print["nb de sites : ",nsite];
Print["tt : ",tt];

(*********************************************)
(****** le span du générateur sur E_perp *****)
cgener={};
grgener={};
igener=2*Mk[[1,1]]; (*definit le nombre de lignes et colonnes à tracer sur M^-k avant sélection *)
(*r=rayonperp[[Min[5,k]]];*)
Do[
Do[
{v=i*generperp[[1]]+j*generperp[[2]];
v={i,j};
cgener=AppendTo[cgener,v];
(*If[traceperp,grgener=AppendTo[grgener,{Black,Disk[v,r]}]]*)
},{j,-igener,igener}],
{i,-igener,igener}];

(****************************)
(*projection de l'hypercube *)
(****************************)

gamainv=Transpose[(Mkinv).gama];
Print["gamainv : ",MatrixForm[gamainv]];
gamaperp=gamainv[[;;,3;;4]];

(*If[traceperp,{
rg=rayonperp[[Min[5,k]]]*2;
grgama={};
Do[grgama=AppendTo[grgama,{Blue,Disk[gamaperp[[i]],rg]}],{i,1,16}];
contour={3,9,4,11,13,8,14,6,3};
bord={};
Do[
bord=AppendTo[bord,{Dashed,Line[{gamaperp[[contour[[i]]]],gamaperp[[contour[[i+1]]]]}]}]
,{i,1,8}]
}];*)

(********************************************************************)
(** construction des sites sur M^-k(E_perp),ainsi que des liens *****)
(********************************************************************)

cperp=Table[0,{i,1,nsite},{j,1,2}]; 
nn=0;
If [kpair,choix=3,choix=13];

av=Reap[
Do[Do[{nn=nn+1;
cperp[[nn]]=gamaperp[[choix]]+{(tt[[j,1]]-tt[[1,1]]+i-1),(j-1)};

Do[{
{existe,ii,jj}=voisins[i+tt[[j,1]]-1,j,ij,etoile,tt,cote1];
If[existe,{numvois=tt[[jj,3]]+ii-tt[[jj,1]];Sow[{nn,numvois}]}]
}
,{ij,1,8}]
},{i,1,tt[[j,2]]-tt[[j,1]]+1}],{j,1,Length[tt]-1}]
]; (* sortie reap*)

liens=av[[2,1]];(*liens contient les paires de sites voisins*)
nliens=Length[liens];

(************** dessin des sites sélectionnés sur E_perp *******************)
(*If[traceperp,
{selection={};
Do[selection=AppendTo[selection,{Red,Circle[cperp[[i]],0.3]}]
,{i,1,nsite}];
Print[" trace dans l'espace perpendiculaire : "];
Graphics[{grgener,grgama,bord,selection}]//Print;
}];*)

(*****************************************************)
(************* retour sur E_parallel *****************)
(*****************************************************)

cpar=Table[f1k*x/2*{cperp[[i,1]]-cperp[[i,2]],cperp[[i,1]]+cperp[[i,2]]},{i,1,nsite}];

(* maille dans E_par *)
maillepar=N[f1k/2*{{1,0},{0,1}}];

(* rescaling global *)
fact=Abs[maillepar[[1,1]]];
cpar=cpar/fact;
maillepar=maillepar/fact;

(*If[tracepar,{
grcpar={};
ray=rayonpar[[Min[5,k]]];
Do[grcpar=AppendTo[grcpar,{Black,Disk[cpar[[i]],ray]}],{i,1,nsite}];
grmaillepar={{Red,Dashed,Line[{{0,0},maillepar[[1]]}]},{Red,Dashed,Line[{maillepar[[1]],maillepar[[1]]+maillepar[[2]] }]},
{Red,Dashed,Line[{maillepar[[1]]+maillepar[[2]],maillepar[[2]]}]},  {Red,Dashed,Line[{maillepar[[2]],{0,0}}]} };*)

Print["************************************************************ "];
Print["Trace de l'approximant dans l'espace parallèle, sans modulo "];
Print["************************************************************ "];

(*Graphics[{grmaillepar,grcpar}]//Print
}];*)



(***********************************************)
(**** sites à l'intérieur de la maille**********)
(***********************************************)

cparmod=Table[{0,0},{i,1,nsite}];
Do[cparmod[[i]]=modulomaille[cpar[[i]]],{i,1,nsite}];

(** preparation du tracé des liens  **)
al=Reap[
Do[{If[(liens[[i,1]]<liens[[i,2]])∧(Norm[cparmod[[liens[[i,1]]]]-cparmod[[liens[[i,2]]]]]<0.5),Sow[{Line[{cparmod[[liens[[i,1]]]],cparmod[[liens[[i,2]]]]}]}]]}
,{i,1,nliens}]
];(*fin reap*)

liensmaille=al[[2,1]];

(*If[traceparmod,{
ray=rayonparmod[[Min[5,k]]];
grcparmod={};
Do[grcparmod=AppendTo[grcparmod,{Black,Disk[cparmod[[i]],ray]}],{i,1,nsite}];*)

Print["************************************************************ "];
Print["Trace de l'approximant dans l'espace parallèle, avec modulo  :"];
Print["************************************************************ "];
(*Graphics[{grmaillepar,grcparmod,liensmaille}]//Print
}];*)

(******************************)
(** version commensurable ***)
(*******************************)

y=0.5;
pipar2com=(({{1, 0, y, -y}, {0, 1, y, y}}))/2;
(* maille dans E_par *)
mailleparcom={{(Mk[[1,1]]+2*y*Mk[[3,1]])/2,0},{0,(Mk[[1,1]]+2*y*Mk[[3,1]])/2}};

(* rescaling global *)
factcom=Abs[mailleparcom[[1,1]]];
cparcom=Table[{0,0},{i,1,nsite}];

Do[{
cparcom[[i]]=pipar2com.Mk.{0,0,cperp[[i,1]],cperp[[i,2]]}},{i,1,nsite}];

cparcom=cparcom/factcom;

mailleparcom=mailleparcom/factcom;

(**************)
cparmodcom=Table[{0,0},{i,1,nsite}];
Do[cparmodcom[[i]]=modulomaille[cparcom[[i]]],{i,1,nsite}];
(**************)

(*If[traceparmodcom,{
(** preparation du tracé des liens  **)
al=Reap[
Do[{If[(liens[[i,1]]<liens[[i,2]])∧(Norm[cparmod[[liens[[i,1]]]]-cparmod[[liens[[i,2]]]]]<0.5),Sow[{Line[{cparmodcom[[liens[[i,1]]]],cparmodcom[[liens[[i,2]]]]}]}]]}
,{i,1,nliens}]
];(*fin reap*)

liensmaillecom=al[[2,1]];
grmailleparcom={{Red,Dashed,Line[{{0,0},mailleparcom[[1]]}]},{Red,Dashed,Line[{mailleparcom[[1]],mailleparcom[[1]]+mailleparcom[[2]] }]},
{Red,Dashed,Line[{mailleparcom[[1]]+mailleparcom[[2]],mailleparcom[[2]]}]},  {Red,Dashed,Line[{mailleparcom[[2]],{0,0}}]} };

ray=rayonparmod[[Min[5,k]]];
grcparmodcom={};
Do[grcparmodcom=AppendTo[grcparmodcom,{Black,Disk[cparmodcom[[i]],ray]}],{i,1,nsite}];*)

Print["************************************************************************************ "];
Print["Trace de l'approximant, version commensurable,  dans l'espace parallèle, avec modulo "];
Print["************************************************************************************ "];
(*Graphics[{grmailleparcom,grcparmodcom,liensmaillecom}]//Print
}];*)
Print[" Fin du calcul "];
Print["  durée CPU depuis le début : ",TimeUsed[]-startcpu, "  durée réelle : ",SessionTime[] -startreel]  ;

liens2

Length[liens]

Length[cparmodcom]

Export["liensPBC8v2.dat",liens];

Export["PtsPBC8v2.dat",cparmodcom];

Length[cparmodcom]

219182.725005`11.491321319858516/3600

cparmodcom[[1]]

liens[[1]]