(* some handy definitions  *)

Clear[lam];  sMat = {{1, lam}, {0,1}}; tMat = {{0, -1}, {1,0}}; sMin =  {{1,  -lam}, {0,1}};  matActs[mat_, xVal_] := Block[{a = mat[[1,1]], b = mat[[1,2]], c = mat[[2,1]], d = mat[[2,2]]},
(a xVal + b)/(c xVal + d)];

(*  CF code --- this is old and far from optimized, but works in this situation --- one could check our results with the Pari code in the appendix of \cite{Rosen-Towse} *)

newQLamCFQn[qVal_Integer,valr_Real,iterBnd_Integer] :=        
  
    Block[ {d = valr,i,j,k,q = qVal,
            perStart =-1,perEnd, accuracyCount,
            iBound =iterBnd,xsi},


   (**   Computes lambda_q continued fractions --   
 give q, real value, and an iteration bound, returns q_n **)
  
     internalAcc = 5000;
     accuracyCount = 0;
     minPoly = MinimalPolynomial[2 Cos[Pi/q], x];
     xsi = {};
     xsiAb = {};
     epR = {};
     qSuite = {1, 0};
     i = 3;
     AppendTo[xsi,d];
     AppendTo[xsiAb,d];
     epsilon = {0};
     lambda = N[2 Cos[Pi/q],internalAcc];
     AppendTo[epR,Round[N[First[xsi] / lambda,internalAcc]]];
     newQn = Last[qSuite]*Abs[ Last[epR] ]*λ  + Sign[Last[epR] ]*qSuite[[-2]];
           AppendTo[ qSuite, newQn]; 
     
     AppendTo[xsiAb, N[1/( First[xsi] -  First[epR] lambda),internalAcc] ];
     AppendTo[xsi,Abs[ Last[xsiAb]] ];
     If[  Last[xsiAb] == First[xsiAb],
            i = iBound; perStart = 0; perEnd = 0 ];                   
     AppendTo[ epR,Sign[ Last[xsiAb]] * Round[N[Last[ xsi]/ lambda,internalAcc]]];
     newQn = Last[qSuite]*Abs[ Last[epR] ]*λ  + Sign[Last[epR] ]*qSuite[[-2]];
           AppendTo[ qSuite, newQn];  
     
     AppendTo[xsiAb,
                  N[(1/(Last[xsi] - 
                      Abs[Last[epR]] lambda)),internalAcc]];
     AppendTo[xsi,Abs[Last[xsiAb]]];
     If[ Last[xsiAb] == ComplexInfinity,
            accuracyCount = i; 
            i = iBound 
           ];     
      Do[  If[ Last[xsiAb] == xsiAb[[j]],
              i = iBound; perStart = j - 1; 
              perEnd = 1 
             ],
         {j, Length[xsiAb] - 1} ];
      While[ (i < iBound),
         AppendTo[epR, Sign[Last[xsiAb]] Round[N[Last[xsi]/ lambda,internalAcc]]];
           newQn = Last[qSuite]*Abs[ Last[epR] ]*λ  + Sign[Last[epR] ]*qSuite[[-2]];
           newQn = PolynomialRemainder[Expand[newQn], minPoly, x];
           AppendTo[ qSuite, newQn]; 
         AppendTo[xsiAb, N[(1/(Last[xsi] - Abs[Last[epR]] lambda)), internalAcc]];
         AppendTo[xsi,Abs[Last[xsiAb]] ];
         iHold = i;
          Do[ If[ Last[xsiAb] == xsiAb[[j]],
               perEnd = i - 1; perStart = j - 1; iHold = iBound],
             {j, Length[xsiAb] - 1}
           ];
          i = iHold;    
          If[Last[xsiAb] == ComplexInfinity,
                  accuracyCount = i; i = iBound ];    
          i = i + 1;
         ];  (* end of while loop *)
      If[ perStart > -1,
          Print[" Periodic: period begins at level = ",
                     perStart, " and ends at ", perEnd] ];
      If[ accuracyCount > 0,
        Print[" Accuracy kept until ",
                    accuracyCount, "  iterations. "] ];
  cutBeginQsuite = Drop[qSuite, {1,2}];
(* λ = 2 Cos[Pi/q];
 minPolyQn = Table[ MinimalPolynomial[cutBeginQsuite[[j]], x],
                                    {j,1, Length[cutBeginQsuite]}];
   hauteLog = Table[ N[Log[Max[Abs[CoefficientList[    minPolyQn[[j]], x ] ] ]]] ,
                                    {j,1, Length[cutBeginQsuite]}];                               
   λ = N[2 Cos[Pi/q],internalAcc]; valeursQnLog =  Table[ Log[N[cutBeginQsuite[[j]]] ],
                                    {j,1, Length[cutBeginQsuite]}];  
   Clear[λ];                             
     rapport = Table[ hauteLog[[j]]/valeursQnLog[[j]], {j, 2, Length[cutBeginQsuite]}];                             
   Return[rapport]  
   *)
   
   
   Return[  epR]
   (*** now want to return list of d-vectors of conjugates of q_n   ***)
   
 ] (* end newQLamCFQn  *)

(************  example for q = 18    ************)

(*  we start with an initial candidate,  of periodic expansion  *)

lam = 2 Cos[Pi/18]; newQLamCFQn[18,N[ (lam + lam^3)/(lam^2 + 3),400],100]

(* we take the matrix corresponding to the period  *)

Clear[lam]; hyperbol = Expand[sMat. tMat.MatrixPower[sMat,4]. tMat.sMin.tMat.MatrixPower[sMin,4].tMat]

(* find its fixed points

Solve[ matActs[hyperbol, x] == x, {x}]

delta = 1 + 4 λ^4

(*  the sqrt of delta will satisfy the polynomial formed by evaluation the min'l poly of delta at x^2;  therefore, the minimal polynomial of the square root divides this latter polynomial.   If this polynomial factors, then we have found a square root in a ``small'' field  --- we thank E. Lanneau for this trick  *)

lam = 2 Cos[Pi/18];Factor[MinimalPolynomial[ 1 +4  lam^4, x^2]]

(* one factor for each sign of square root, we take the easier one to work with  *)

Solve[ 0 == -37 - 33 b - 3 d +  b d, {b}]

(* beta is our square root *)

beta = (37 + 3 delta)/(delta - 33)

λ =lam ;FullSimplify[beta^2 - delta]  (* check that this is indeed a square root   *)

Clear[λ]; Together[beta]   (* now looking for a nice form for beta  *)

Clear[λ]; Together[((-1+2 (λ^2)+beta)/(4 λ))]   (* here is the fixed point  *)

(*  we reduce the degrees of the polynomials involved by using the minimal polynomial of λ as a modulus   *) 
PolynomialRemainder[9-8 (λ^2)+(λ^4)+(λ^6), MinimalPolynomial[ 2 Cos[Pi/18], λ], λ]

newQLamCFQn[18,N[ (12-17 (lam^2)+7 (lam^4))/(2 lam (-8+(lam^4))),400],100]  (* here's  a check that the fixed point has purely periodic Rosen {q=18}-expansion  *)

(************  exemples for q = 20    ************)

lam = 2 Cos[Pi/20]; newQLamCFQn[20,N[ (lam + lam^3)/(lam^2 + 4), 900],100]

(* thus {1,-1,-1,-1,4,-1,-1,-1}  should be the period  *)

Clear[lam];  hyperbol = Expand[MatrixPower[sMat.tMat,4]. MatrixPower[sMin,4].tMat.MatrixPower[sMin.tMat,3]]

Solve[ matActs[hyperbol, x] == x, {x}]

lam = 2 Cos[Pi/20]; newQLamCFQn[20,N[(-4 lam+10 (lam^3)-6 (lam^5)+(lam^7)+(Sqrt[4 (lam^2)-4 (lam^4)+17 (lam^6)-32 (lam^8)+24 (lam^10)-8 (lam^12)+(lam^14)]))/(2 (4 (lam^2)-4 (lam^4)+(lam^6))), 900],100]

Clear[lam]; Factor[4 (lam^2)-4 (lam^4)+17 (lam^6)-32 (lam^8)+24 (lam^10)-8 (lam^12)+(lam^14)]

lam = 2 Cos[Pi/20];Factor[MinimalPolynomial[1+4 (lam^4)-4 (lam^6)+(lam^8), x^2]]

(*  thus square root t is  (a^2 - 29 a + 41)/(4 a - 16)   *)

Clear[λ, δ];  beta = (δ^2 - 29 δ + 41)/(4 δ - 16)

λ = lam; δ = 1+4 (λ^4)-4 (λ^6)+(λ^8);FullSimplify[beta^2 - δ  ]

Clear[λ]; δ = 1+4 (λ^4)-4 (λ^6)+(λ^8); Together[beta]

Clear[λ];beta = Simplify[(PolynomialRemainder[Numerator[Together[beta]], MinimalPolynomial[lam,λ], λ])/(PolynomialRemainder[Denominator[Together[beta]], MinimalPolynomial[lam,λ], λ])]

λ = lam;FullSimplify[beta^2 - δ  ]

Clear[lam]

(Factor[-4 lam+10 (lam^3)-6 (lam^5)+(lam^7)]+(Sqrt[Factor[4 (lam^2)-4 (lam^4)+17 (lam^6)-32 (lam^8)+24 (lam^10)-8 (lam^12)+(lam^14)]]))/(2  Factor[4 (lam^2)-4 (lam^4)+(lam^6)])

Simplify[((lam (-2+(lam^2)) (2-4 (lam^2)+(lam^4))+lam (-2+(lam^2))(Sqrt[(1+4 (lam^4)-4 (lam^6)+(lam^8))]))/(2 (lam^2) (((-2+(lam^2)))^2)))]

Clear[λ]; point =Together[ -((2-4 (λ^2)+(λ^4)+beta)/(4 λ-2 (λ^3)))]

lam = 2 Cos[Pi/20];point = Simplify[(PolynomialRemainder[Numerator[Together[point]], MinimalPolynomial[lam,λ], λ])/(PolynomialRemainder[Denominator[Together[point]], MinimalPolynomial[lam,λ], λ])]

λ = lam;newQLamCFQn[20,N[point, 500],100]

(************  exemples pour q = 24    ************)

lam = 2 Cos[Pi/24]; newQLamCFQn[24,N[ ( 2 lam^2 +13)/(lam + lam^3), 5000],200]

preperiod = Take[%, 49]

factors= {}; Clear[lam];  (* here we decided to simply ``peel off'' the perperiod

AppendTo[factors, tMat.MatrixPower[sMin, 4]. tMat.sMin.tMat.sMin.tMat.MatrixPower[sMat, 3]. tMat.sMin];
AppendTo[factors,  tMat.MatrixPower[sMin.tMat, 5]. MatrixPower[sMin, 2] ] ;
 AppendTo[factors,  tMat.MatrixPower[sMin.tMat, 5].MatrixPower[sMin, 2]] ; 
 AppendTo[factors,  tMat. sMat.tMat. sMin.tMat.sMin] ; 
 AppendTo[factors, tMat. sMat.MatrixPower[tMat.sMin, 7].tMat. sMin] ; 
 AppendTo[factors,MatrixPower[tMat.sMin, 8] .tMat. MatrixPower[ sMin, 2]] ; 
AppendTo[factors, MatrixPower[tMat.sMin, 2] . tMat. MatrixPower[ sMin, 2] . MatrixPower[tMat.sMin, 3] .tMat. MatrixPower[ sMin, 14] ] ; 
AppendTo[factors,tMat. sMat.tMat.sMat.tMat.sMat. tMat.sMat] ;

prodMat = IdentityMatrix[2]; Do[  prodMat = factors[[j]]. prodMat, {j,1, Length[factors]}]; prodMat = Expand[prodMat]; 
a = PolynomialRemainder[prodMat[[1,1]], MinimalPolynomial[ 2 Cos[Pi/24], lam], lam];
b = PolynomialRemainder[prodMat[[1,2]], MinimalPolynomial[ 2 Cos[Pi/24], lam], lam];
c = PolynomialRemainder[prodMat[[2,1]], MinimalPolynomial[ 2 Cos[Pi/24], lam], lam];
d = PolynomialRemainder[prodMat[[2,2]], MinimalPolynomial[ 2 Cos[Pi/24], lam], lam];
prodMat = {{a, b}, {c, d}}

Clear[lam];periodicValue = Together[ -matActs[prodMat,   ( 2 lam^2 +13)/(lam + lam^3)]]; 
periodicValue = (PolynomialRemainder[Expand[Numerator[periodicValue]], MinimalPolynomial[2Cos[Pi/24], lam], lam])/ (PolynomialRemainder[Expand[Denominator[periodicValue]], MinimalPolynomial[2Cos[Pi/24], lam], lam])

lam = 2 Cos[Pi/24];newQLamCFQn[24,N[periodicValue , 5000],200]