function L(x5,n)
%
%  Compute the matrix L(x) and det(I-L) see Section 3.1 in the paper
%  S.Yu.Orevkov, Asymptotics of the number of lattice triangulations of rectangles of width 4 and 5,
%  arXiv:2412.17065v1
%
%  We solve two equations K.G_lambda=F_lambda, lambda=1,2, where G[lambda] is the vector whose
%  entries are the values of g_{lambda,1},g_{lambda,2},h_{lambda,3},h_{lambda,4},h_{lambda,5},j_{lambda,2}
%  at the n-th roots of unity (shifted by 1/4n; see Section 4.1) and F[lambda] is the rhs of the
%  system of equations (3.8), (3.20)-(3.24).
%  In the MATLAB implementation, we define an n\times 2 matrix F whose columns are F_1 and F_2, and
%  compute G by the command G=F\K.
%
%  Example:
%           L(0.02547036302369,26)
%
%  (the determinant of I-L should be close to 0)
%
  x = x5^(1/5);
  z = exp(2*pi*i/n);
  z4= exp(pi*i/n/2);
  nn = 2*n^2*(n+1) + n*(n+1)/2 + n;
  K = eye(nn); F = zeros(nn,2);
  for j=1:nn
    K(j,j)=1;
  end;
  b=1/3;

  for iu=1:n
    for iv=1:n
      zu=z4*z^iu; u = (zu+b)/(b*zu+1);  zv=z4*z^iv; v = (zv+b)/(b*zv+1);

      r1uv  = 1 - x/v - v^2*x^2/u; %% r1(u,v)

      for is=1:n
        zs=z4*z^is; s = (zs+b)/(b*zs+1);                                % eq. (3.20):  g1(s,u,v) = ...

        p1suv = 1 - x/v - v^2*x^2/u - u^2*x^2/s/v + u^2*x^3/s/v^2; %% p1(s,u,v)
        % qTsuv_X_T2 is the list of coefficients of q(T,s,u,v)T^2 (as a polynomial in T)
        % p1usT_X_T2 is the list of coefficients of p(u,s,T)T^2 (as a polynomial in T)
        qTsuv_X_T2 = [-x^2*p1suv/s, 0, p1suv, -x*p1suv - x^2*s^2*r1uv/u, x^3*s^2*r1uv/u];
        p1usT_X_T2 =  [-x^2/s,      0,    1,     - x - x^2*s^2/u,     x^3*s^2/u    ];

        Psi1suv = Psi1(x,s,u,v);

        iG1 = indG1(is,iu,iv,n);
        K(iG1,indH4(is,iv,n))=K(iG1,indH4(is,iv,n)) - x*u*(1-x/v)/p1suv;
        K(iG1,indH5(is,iu,n))=K(iG1,indH5(is,iu,n)) - x*v/p1suv;
        F(iG1,1) = 1/(x^2*Psi1suv*p1suv);                       %% the r.h.s. with delta_{\labmda,1}

        for it=1:n
          zt=z4*z^it; t = (zt+b)/(b*zt+1); dt = zt*(1-b^2)/(b*zt+1)^2/n;

          qtsuv = polyval(qTsuv_X_T2,t)/t^2;   %% q(t,s,u,v)
          p1ust  = polyval(p1usT_X_T2,t)/t^2;  %% p1(u,s,t)
          r1st  = 1 - x/t - (t*x)^2/s;         %% r1(s,t)

          K(iG1,indG1(iu,is,it,n))=K(iG1,indG1(iu,is,it,n)) - x^2*v/s        /Psi1suv*(  p1ust/qtsuv - 1/p1suv)*dt;
          K(iG1,indG2(iv,is,it,n))=K(iG1,indG2(iv,is,it,n)) - x^2*u*(1-x/v)/s/Psi1suv*(   r1st/qtsuv - 1/p1suv)*dt;
          K(iG1,indG2(it,iu,iv,n))=K(iG1,indG2(it,iu,iv,n)) - x^2*r1uv       /Psi1suv*((1-x/t)/qtsuv + 1/p1suv)*dt;

          K(iG1,indH4(iu,it,n))=K(iG1,indH4(iu,it,n)) + x^3*v      /Psi1suv*( (1-x/t)/qtsuv + 1/p1suv )*dt;
          K(iG1,indH3(it,iv,n))=K(iG1,indH3(it,iv,n)) - x^3*(1-x/v)/Psi1suv/p1suv*dt;

        end;
      end;

      for it=1:n
        zt=z4*z^it; t = (zt+b)/(b*zt+1);                                     % eq. (3.21):  g2(t,u,v) = .....

        p2tuv = (1-x/t)*r1uv; %% p2(t,u,v)
        qtSuv_X_S2 = [-x^2*p2tuv/t/u, 0, p2tuv, -x^2*(u^2/v*(1-x/t)*(1-x/v) + t^2*r1uv), x^4*t^2*u^2/v*(1-x/v)];
        % qtSuv_X_S2 is the list of coefficients of q(t,S,u,v)S^2 (as a polynomial in S)
        
        Psi2tuv = Psi2(x,t,u,v);
        iG2 = indG2(it,iu,iv,n);

        K(iG2,indH3(it,iv,n))=K(iG2,indH3(it,iv,n)) - x*(1-x/v)/r1uv/Psi2tuv;
        K(iG2,indH4(iu,it,n))=K(iG2,indH4(iu,it,n)) - x*v/r1uv;

        for is=1:n
           zs=z4*z^is; s = (zs+b)/(b*zs+1); ds = zs*(1-b^2)/(b*zs+1)^2/n;

           qtsuv = polyval(qtSuv_X_S2,s)/s^2;  %% q(t,s,u,v)
           r1st  = 1 - x/t - (t*x)^2/s;        %% r1(s,t)
           p1ust = r1st - (s*x)^2/u/t*(1-x/t); %% p1(u,s,t)
           p1suv = r1uv - (u*x)^2/s/v*(1-x/v); %% p1(s,u,v)
           p2vst = (1-x/v)*r1st;               %% p2(v,s,t)

           K(iG2,indG1(is,iu,iv,n))=K(iG2,indG1(is,iu,iv,n)) - x^2/u    /Psi2tuv*( p1suv/qtsuv + 1/(1-x/t))*ds;
           K(iG2,indG1(iu,is,it,n))=K(iG2,indG1(iu,is,it,n)) - x^2*v/t/u/Psi2tuv*( p1ust/qtsuv - 1/r1uv   )*ds;
           K(iG2,indG2(iv,is,it,n))=K(iG2,indG2(iv,is,it,n)) - x^2/t    /Psi2tuv*( p2vst/qtsuv            )*ds;

           K(iG2,indH4(is,iv,n))=K(iG2,indH4(is,iv,n)) + x^3*(1-x/v) /Psi2tuv*(1/qtsuv          )*ds;
           K(iG2,indH5(is,iu,n))=K(iG2,indH5(is,iu,n)) + x^3*v/u     /Psi2tuv*(1/qtsuv + 1/p2tuv)*ds;

        end;
      end;
    end;
  end;

  % ----------------------------------------------------------------------------------------------

  for is=1:n                                                         % eq. (3.22):  h4(s,v) = .....
     zs=z4*z^is; s = (zs+b)/(b*zs+1);

     for iv=1:n
        zv=z4*z^iv; v = (zv+b)/(b*zv+1);

        Psi13sv = Psi13(x,s,v,30*n);

        iH4 = indH4(is,iv,n);

        F(iH4,2) = x^2/s/(1-x/v)/Psi13sv;       %% the r.h.s. with delta_{\labmda,2}

        for iu=1:n
           zu=z4*z^iu; u = (zu+b)/(b*zu+1); du = zu*(1-b^2)/(b*zu+1)^2/n;

           r1uv  = 1 - x/v - v^2*x^2/u;        %% r1(u,v)
           p1suv = r1uv - u^2*x^2/s/v*(1-x/v);  %% p1(s,u,v)
           qTsuv_X_T2 = [-x^2*p1suv/s, 0, p1suv, -x*p1suv - x^2*s^2*r1uv/u, x^3*s^2*r1uv/u];

           p1usT_X_T2 =  [-x^2/s, 0,  1,    - x - s^2*x^2/u,    s^2*x^3/u];
        %  p1usT_X_T2 is the list of coefficients of p(u,s,T)T^2 (as a polynomial in T)

           Phi1suv = Phi1(x,s,u,v);
           Psi1suv = 1 - (x^2/s)*p1suv*Phi1suv; %% Psi1(u,s,t)

           K(iH4,indG1(is,iu,iv,n))=K(iH4,indG1(is,iu,iv,n)) - x/s/v /Psi13sv*(                   Psi1suv )*du;
           K(iH4,indH5(is,iu,n)   )=K(iH4,indH5(is,iu,n)   ) - x^2/s /Psi13sv*( 1/(1-x/v) + x^2/s*Phi1suv )*du;

           for it=1:n
             zt=z4*z^it; t = (zt+b)/(b*zt+1); dt = zt*(1-b^2)/(b*zt+1)^2/n;

             qtsuv = polyval(qTsuv_X_T2,t)/t^2; %% q(t,s,u,v)
             p1ust = polyval(p1usT_X_T2,t)/t^2; %% p1(u,s,t)
             p2tuv = (1-x/t)*r1uv;              %% p2(t,u,v)

             K(iH4,indG1(iu,is,it,n))=K(iH4,indG1(iu,is,it,n)) + x^3/s^2/Psi13sv*(1/(1-x/v) + p1ust/qtsuv)*du*dt;
             K(iH4,indG2(it,iu,iv,n))=K(iH4,indG2(it,iu,iv,n)) + x^3/s/v/Psi13sv*(            p2tuv/qtsuv)*du*dt;
             K(iH4,indH4(iu,it,n)   )=K(iH4,indH4(iu,it,n)   ) - x^4/s  /Psi13sv*(          (1-x/t)/qtsuv)*du*dt;

           end;
        end;

        for it=1:n
           zt=z4*z^it; t = (zt+b)/(b*zt+1); dt = zt*(1-b^2)/(b*zt+1)^2/n;

           Psi2vst = Psi2(x,v,s,t);

           K(iH4,indG2(iv,is,it,n))=K(iH4,indG2(iv,is,it,n)) - x/s /Psi13sv*( Psi2vst )*dt;

        end;
     end;
  end;

  % --------------------------------------------------------------------------------------

  for is=1:n                                                 % eq. (3.23):  h5(s,u) = .....
     zs=z4*z^is; s = (zs+b)/(b*zs+1);

     for iu=1:is
        zu=z4*z^iu; u = (zu+b)/(b*zu+1);

        Psi14su = Psi14(x,s,u,30*n);
        iH5 = indH5(is,iu,n);
        for it=1:n
           zt=z4*z^it; t = (zt+b)/(b*zt+1); dt = zt*(1-b^2)/(b*zt+1)^2/n;

           p1ust = 1 - x/t - t^2*x^2/s - s^2*x^2/u/t + s^2*x^3/u/t^2; %% p1(u,s,t)
           r1st  = 1 - x/t - t^2*x^2/s;                               %% r1(s,t)

           qtsuV_X_V2 = [-x^2*p1ust/u, 0, p1ust, -x*p1ust - x^2*u^2*r1st/s, x^3*u^2*r1st/s];
         % qtsuV_X_V2 is the list of coefficients of q(t,s,u,V)V^2 (as a polynomial in V)


           Phi1ust = Phi1(x,u,s,t);
           Psi1ust = 1 - (x^2/u)*p1ust*Phi1ust; %% Psi1(u,s,t);

           K(iH5,indG1(iu,is,it,n))=K(iH5,indG1(iu,is,it,n)) - x/s /Psi14su*(                   Psi1ust )*dt;
           K(iH5,indH4(iu,   it,n))=K(iH5,indH4(iu,   it,n)) - x^2 /Psi14su*( 1 + x^2/u*(1-x/t)*Phi1ust )*dt;

           for iv=1:n
              zv=z4*z^iv; v = (zv+b)/(b*zv+1); dv = zv*(1-b^2)/(b*zv+1)^2/n;

              qtsuv = polyval(qtsuV_X_V2,v)/v^2; %% q(t,s,u,v)
              r1uv  = 1 - x/v - v^2*x^2/u;       %% r1(u,v)
              p2tuv = (1-x/t)*r1uv;              %% p2(y,u,v)
              p2vst = (1-x/v)*r1st;              %% p2(v,s,t)

              K(iH5,indG2(it,iu,iv,n))=K(iH5,indG2(it,iu,iv,n)) + x^3/u /Psi14su*( 1 + p2tuv/qtsuv )*dt*dv;
              K(iH5,indG2(iv,is,it,n))=K(iH5,indG2(iv,is,it,n)) + x^3/s /Psi14su*( 1 + p2vst/qtsuv )*dt*dv;
           end;
        end;

        for iv=1:n
           zv=z4*z^iv; v = (zv+b)/(b*zv+1); dv = zv*(1-b^2)/(b*zv+1)^2/n;

           p1suv = 1 - x/v - v^2*x^2/u - u^2*x^2/s/v + u^2*x^3/s/v^2; %% p1(s,u,v)

           Phi1suv = Phi1(x,s,u,v);
           Psi1suv = 1 - (x^2/s)*p1suv*Phi1suv; %% Psi1(s,u,v)

           K(iH5,indG1(is,iu,iv,n))=K(iH5,indG1(is,iu,iv,n)) - x/u /Psi14su*(                   Psi1suv )*dv;
           K(iH5,indH4(is,   iv,n))=K(iH5,indH4(is,   iv,n)) - x^2 /Psi14su*( 1 + x^2/s*(1-x/v)*Phi1suv )*dv;
        end;
     end;
  end;

  % ----------------------------------------------------------------------------------------

  for it=1:n
     zt=z4*z^it; t = (zt+b)/(b*zt+1);

     for iv=1:n                                                  % eq. (3.8):  h3(t,v) = .....
        zv=z4*z^iv; v = (zv+b)/(b*zv+1);

        iH3 = indH3(it,iv,n);
        K(iH3,indJ2(it,n))=K(iH3,indJ2(it,n)) - x/(1-x/v);

        for is=1:n
           zs=z4*z^is; s = (zs+b)/(b*zs+1); ds = zs*(1-b^2)/(b*zs+1)^2/n;

           K(iH3,indG2(iv,is,it,n))=K(iH3,indG2(iv,is,it,n)) - x/t*ds;
        end;

     end;

     iJ2 = indJ2(it,n);

     for is=1:n                                                      % eq. (3.24):  j2(t) = .....
        zs=z4*z^is; s = (zs+b)/(b*zs+1); ds = zs*(1-b^2)/(b*zs+1)^2/n;

        r1st  = 1 - x/t - (t*x)^2/s;   %% r1(s,t)

        F(iJ2,2)=F(iJ2,2) + x^3/(s*r1st*Psi0(x,s))*ds;     %% the r.h.s. with delta_{\labmda,2}

        for iu=1:n
           zu=z4*z^iu; u = (zu+b)/(b*zu+1); du = zu*(1-b^2)/(b*zu+1)^2/n;

           K(iJ2,indG1(iu,is,it,n))=K(iJ2,indG1(iu,is,it,n)) - x^2/t/s*ds*du;
        end;
     end;
  end;

  nn = 2*n^2*(n+1) + n*(n+1)/2 + n;
  G = K\F;

  L11 = 0; L21 = 0; L22 = 0; L12 = 1;
  for it=1:n
     zt=z4*z^it; t = (zt+b)/(b*zt+1); dt = zt*(1-b^2)/(b*zt+1)^2/n;
     L11 = L11 + x*G(indJ2(it,n),1)*dt;
     L21 = L21 + x*G(indJ2(it,n),2)*dt;
     for iu=1:n
        zu=z4*z^iu; u = (zu+b)/(b*zu+1); du = zu*(1-b^2)/(b*zu+1)^2/n;
        L12 = L12 + x^2*( x^2*G(indH3(it,iu,n),1) + G(indH4(iu,it,n),1) )*dt*du;
        L22 = L22 + x^2*( x^2*G(indH3(it,iu,n),2) + G(indH4(iu,it,n),2) )*dt*du;
     end;
  end;
  L11 = L11*x5; L12 = L12*x5;

  disp((1-L11)*(1-L22)-L12*L21); 
  L = [L11,L12; L21,L22]
end

% ============================================================================================
% ============================================================================================
% ============================================================================================

function ps13 = Psi13(x,s,v,n)
   ps13 = 1 - x^4*(1-x/v)/s^2/v * Phi13(x,s,v,n);
end

function ph13 = Phi13(x,s,v,n)
   w = exp(2*pi*i/n); b=0; S=0;
   for k=1:n
      z=w^k; u=(z+b)/(b*z+1); du=z*(1-b^2)/(b*z+1)^2; S=S+u*du*Phi1(x,s,u,v);
   end
   ph13 = S/n;
end

function ps14 = Psi14(x,s,u,n)
   ps14 = 1 - x^4/s/u * Phi14(x,s,u,n);
end

function ph14 = Phi14(x,s,u,n)
   w = exp(2*pi*i/n); b=0; S=0;
   for k=1:n
      z=w^k; v=(z+b)/(b*z+1); dv=z*(1-b^2)/(b*z+1)^2; S=S+v*dv*Phi1(x,s,u,v);
   end
   ph14 = S/n;
end

function ps0 = Psi0(x,u)
  ps0 = 1 - x^2/u * Phi0(x,u);
end

function ph0 = Phi0(x,u)
  r = [-x^2/u, 0, 1, -x];
  dr = [-3*x^2/u, 0, 1]; v = roots(r);
  if abs(v(3))>1
    disp("Error in Phi0");
  end;
  ph0 = v(3)^2/polyval(dr,v(3));
end

function ps1 = Psi1(x,s,u,v)
  p1suv = 1 - x/v - v^2*x^2/u - u^2*x^2/s/v + u^2*x^3/s/v^2;
  ps1 = 1 - x^2/s*p1suv * Phi1(x,s,u,v);
end

function ph1 = Phi1(x,s,u,v)
   p1suv = 1 - x/v - v^2*x^2/u - u^2*x^2/s/v + u^2*x^3/s/v^2; r1uv  = 1 - x/v - v^2*x^2/u;
   q = [-x^2*p1suv/s, 0, p1suv, -x*p1suv - x^2*s^2*r1uv/u, x^3*s^2*r1uv/u];
   dq = [ 4*q(1), 0, 2*q(3), q(4) ]; t = roots(q);
   if abs(t(3))>1
     disp("Error in Phi1");
   end;
   if abs(t(4))>1
     disp("Error in Phi1");
   end;
   ph1 = t(3)^3/polyval(dq,t(3)) + t(4)^3/polyval(dq,t(4));
end

function ps2 = Psi2(x,t,u,v)
   r1uv =  1 - x/v - v^2*x^2/u;  p2tuv = (1-x/t)*r1uv;
   ps2 = 1 - x^2/t/u*p2tuv * Phi2(x,t,u,v);
end

function ph2 = Phi2(x,t,u,v)
   r1uv =  1 - x/v - v^2*x^2/u;  p2tuv = (1-x/t)*r1uv;
   q = [-x^2*p2tuv/t/u, 0, p2tuv, -x^2*(u^2/v*(1-x/t)*(1-x/v) + t^2*r1uv), x^4*t^2*u^2/v*(1-x/v)];
   dq = [ 4*q(1), 0, 2*q(3), q(4) ]; s = roots(q);
   if abs(s(3))>1
     disp("Error in Phi2");
   end;
   if abs(s(4))>1
     disp("Error in Phi2");
   end;
   ph2 = s(3)^3/polyval(dq,s(3)) + s(4)^3/polyval(dq,s(4));
end

function res=indG1(i,j,k,n)
         res=n^2*(i-1)+n*(j-1)+k;
end

function res=indG2(i,j,k,n)
         res=n^2*(i-1)+n*(j-1)+k + n^3;
end

function res=indH3(i,j,n)
         res=n*(i-1)+j + 2*n^3;
end

function res=indH4(i,j,n)
         res=n*(i-1)+j + 2*n^3 + n^2;
end

function res=indH5(i,j,n)
    if j<i 
        ij = 2*n^2*(n+1)+i*(i-1)/2+j;
    else
        ij = 2*n^2*(n+1)+j*(j-1)/2+i;
    end;
    res=ij;
end

function res=indJ2(i,n)
         res=2*n^2*(n+1) + n*(n+1)/2 + i;
end

