function res = J(x4,n)
%
%  Compute the function J(x) using the Riemann sums over the n-th roots of unity;
%  see Section 2.1 of the paper
%  S.Yu.Orevkov, "Asymptotics of the number of lattice triangulations of rectangles of width 4 and 5",
%  arXiv:2412.17065v1
%
%  We solve the equation K.G=F, where G is the vector whose entries are the values of g_1,g_2,h_3
%  at the n-th roots of unity (shifted by 1/4; see Section 4.1) and F is the rhs of the system 
%  of equations (2.11)-(2.13). Since g_2 is symmetric, i.e., g_2(t,u) = g_2(u,t), we compute it
%  only for (t,u) such that t=exp(i*theta), u=exp(i*phi), theta >= phi. Thus the size of the
%  system is n^2 [size of g_1] + n(n+1)/2 [size of g_2] + n [size of h_3] = 3n(n+1)/2
%
%  Example:
%           J(0.054114323851836479442, 40)
%
%  (the result should be close to 1)
%
  x = x4^(1/4); b=1/3; z = exp(2*pi*i/n); z4 = exp(pi*i/n/2);
  K = eye(3*n*(n+1)/2); F = zeros(3*n*(n+1)/2,1);
  for is=1:n                                                        % eq. (2.11):  g1(s,u) = .......
    for iu=1:n
      zs=z4*z^is; zu=z4*z^iu; s = (zs+b)/(b*zs+1); u = (zu+b)/(b*zu+1);
      psu = 1 - x/u - x^2*u^2/s; %% p(s,u)
      rsu = (1 - x/s)*(1 - x/u); %% r(s,u)
      qTsu_X_T2 = [ -x^2*psu/s, 0,  psu, -x*psu + (x/u-1)*(s*x)^2/u, (u/x-1)*(s*x^2/u)^2 ];
    % qTsu_X_T2 is the list of coefficients of q(T,s,u)T^2 (as a polynomial in T)
      phi1su = Phi1(x,s,u);
      psi1su = 1 - x^2*psu*phi1su/s;
      F( indG1(is,iu,n) ) = s/psu/psi1su;
      K( indG1(is,iu,n), indH3(is,n) ) = K( indG1(is,iu,n), indH3(is,n) ) - x^3*u/psu;
      for k=1:n
        zt=z4*z^k; t = (zt+b)/(b*zt+1); dt = zt*(1-b^2)/(b*zt+1)^2/n;
        qtsu = polyval(qTsu_X_T2,t)/t^2; %% q(t,s,u)
        pst = 1 - x/t - x^2*t^2/s;       %% p(s,t)
%
        K(indG1(is,iu,n),indG1(is,k,n))=K(indG1(is,iu,n),indG1(is,k,n)) + x^2*u/s/psi1su*( 1/psu - pst/qtsu )*dt;
        K(indG1(is,iu,n),indG2(k,iu,n))=K(indG1(is,iu,n),indG2(k,iu,n)) - x^2*s*(1-x/u)/psi1su*( 1/psu + (1-x/t)/qtsu )*dt;
      end;
    end
  end;
%
  for it=1:n                                                        % eq. (2.12):  g2(t,u) = .......
    for iu=1:n
      zt=z4*z^it; zu=z4*z^iu; t = (zt+b)/(b*zt+1); u = (zu+b)/(b*zu+1);
      ptu = 1 - x/u - x^2*u^2/t; %% p(t,u)
      rtu = (1 - x/t)*(1 - x/u); %% r(t,u)
      qtSu_X_S2 = [-x^2*rtu/u/t, 0, rtu, x^2*( u^2*(x/t-1) + t^2*(x/u-1) ), (u*t*x^2)^2 ];
    % qtSu_X_S2 is the list of coefficients of q(t,S,u)S^2 (as a polynomial in S)
      phi2tu = Phi2(x,t,u);
      psi2tu = 1 - x^2*rtu*phi2tu/t/u;
      for k=1:n
        zs=z4*z^k; s = (zs+b)/(b*zs+1); ds = zs*(1-b^2)/(b*zs+1)^2/n;
        qtsu = polyval(qtSu_X_S2,s)/s^2; %% q(t,s,u)
        psu = 1 - x/u - x^2*u^2/s;       %% p(s,u)
        pst = 1 - x/t - x^2*t^2/s;       %% p(s,t)
%
        if it >= iu
          K(indG2(it,iu,n),indG1(k,iu,n))=K(indG2(it,iu,n),indG1(k,iu,n)) - x^2/u/s/psi2tu*( 1/(1-x/t) + psu/qtsu )*ds;
          K(indG2(it,iu,n),indG1(k,it,n))=K(indG2(it,iu,n),indG1(k,it,n)) - x^2/t/s/psi2tu*( 1/(1-x/u) + pst/qtsu )*ds;
          K(indG2(it,iu,n),indH3(k,n))   =K(indG2(it,iu,n),indH3(k,n))    + x^5/s/qtsu/psi2tu*ds;
        end;
      end;
    end
  end;
%
  for is=1:n                                              %  eq. (2.13):  h3(s) = ......
    zs = z4*z^is; s = (zs+b)/(b*zs+1);
    psi3s = Psi3(x,s,6*n);
    for iu=1:n
      zu = z4*z^iu; u = (zu+b)/(b*zu+1); du = zu*(1-b^2)/(b*zu+1)^2/n;
      psu = 1 - x/u - x^2*u^2/s; %% p(s,u)
      K(indH3(is,n),indG1(is,iu,n))=K(indH3(is,n),indG1(is,iu,n)) - 2*Psi1(x,s,u)/x/s/psi3s*du;
      qTsu_X_T2 = [ -x^2*psu/s, 0,  psu, -x*psu + (x/u-1)*(s*x)^2/u, (u/x-1)*(s*x^2/u)^2 ];
    % qTsu_X_T2 is the list of coefficients of q(T,s,u)T^2 (as a polynomial in T)
      for it=1:n
        zt = z4*z^it; t = (zt+b)/(b*zt+1); dt = zt*(1-b^2)/(b*zt+1)^2/n; qust=polyval(qTsu_X_T2,t)/t^2;
        K(indH3(is,n),indG2(iu,it,n))=K(indH3(is,n),indG2(iu,it,n)) + x*(1-x/u)*(1-x/t)/qust/psi3s*du*dt;
      end;
    end;
  end;
%
  G = K\F; int=0;
  for is=1:n
    for iu=1:n
      zs=z4*z^is; zu=z4*z^iu; 
      int=int+G( indG1(is,iu,n) )*zs*(1-b^2)/(b*zs+1)/(zs+b) * zu*(1-b^2)/(b*zu+1)/(zu+b);
    end;
  end;
  res = x4*int/n^2;
end

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

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

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

function ps3 = Psi3(x,s,n)
  ps3 = 1 - x^4*Phi3(x,s,n)/s^2;
end

function ph1 = Phi1(x,s,u)
   q = [ x^2*(x/u + u^2*x^2/s - 1)/s, 0, 1 - x/u - u^2*x^2/s,... 
         x*(x/u - 1 - s^2*x/u + u^2*x^2/s + s^2*x^2/u^2), s^2*x^3*(1 - x/u)/u ];
   dq = [ 4*q(1), 0, 2*q(3), q(4) ];
   t = roots(q);
   ph1 = t(3)^3/polyval(dq,t(3)) + t(4)^3/polyval(dq,t(4));
end

function ph2 = Phi2(x,t,u)
   q = [ -x^2*(1 - x/t)*(1 - x/u)/u/t, 0, (1 - x/t)*(1 - x/u),...
          x^2*( u^2*(x/t-1) + t^2*(x/u-1) ), u^2*t^2*x^4 ];
   dq = [ 4*q(1), 0, 2*q(3), q(4) ];
   s = roots(q);
   ph2 = s(3)^3/polyval(dq,s(3)) + s(4)^3/polyval(dq,s(4));
end

function ph3 = Phi3(x,s,n)
   w = exp(2*pi*i/n); b=1/3; int=0;
   for k=1:n
      z=w^k; u=(z+b)/(b*z+1); du=z*(1-b^2)/(b*z+1)^2; int=int+u*du*Phi1(x,s,u);
   end
   ph3 = int/n;
end

function res = indG1(k,l,n)
         res = n*(k-1)+l;
end

function res = indG2(k,l,n)
    if l<k
         kl=n^2+k*(k-1)/2+l;
    else
         kl=n^2+l*(l-1)/2+k;
    end; res=kl;
end

function res = indH3(k,n)
         res = n*(3*n+1)/2+k;
end

% xxx0 = 0.0541143238518364794421

