\\ The function Rsqlatt(R) finds the grid point [x,y]
\\ that leads to the best approximation of the given resistance distance R in ohms between [0,0] and [x,y].
\\ Author: Hugo Pfoertner, 2022-09-09
\\ Examples: Rsqlatt(2) -> [106, 8], Rsqlatt(Pi) -> [2853, 2567], Rsqlatt(7) -> [636875169, 303176603].
\\ Increased precision, e.g. \p100, recommended for R > 10.

\\ Auxiliary functions

\\ Is n a sum of 2 squares (after Charles R Greathouse IV in A001481)
issum2sq(n) = {
  my (f=factor(n)); for(i=1, #f[, 1], if(f[i, 2]%2 && f[i, 1]%4==3, return(0))); 1
};

\\ Decomposition of primes of the form 4*k+1 into 2 squares
decomp2sq(p) = {
  my (m=(p-1)/4, r, x, limit=ceil(sqrt(p)));
  if (p > 4 && denominator(m)==1,
  forprime (c=2, oo, if (!issquare(Mod(c,p)), r=c;break));
  x = lift(Mod(r,p)^m);
  until (p < limit, r=p%x; p=x; x=r);
  if (p^2+x^2==4*m+1, p+I*x,0) ,0)
};

\\ Representations of n as sums of 2 squares n = a^2 + b^2.
\\ If none exists, output [[0,0]],
\\ otherwise output of all d decompositions [[a_1,b_1],[a_2,b_2],...,[a_d,b_d]] with a_i>=b_i and a_{i+1}>a_i.
d2sq(n) = {
  if (n<1, [[0,0]], if(n==1, [[1,0]],
  my (L=List(), V=List(), F=factor(n), Q=1, P=1, npf=#F[,1], npf0=1, P0=1);
  if (F[1,1]==2, Q=2^(F[1,2]\2); npf0=2; if (F[1,2]%2, P0+=I));
  for (np=npf0, npf, if (F[np,1]%4==3, if (F[np,2]%2==0, Q*=F[np,1]^(F[np,2]/2), return([[0,0]]))));
  for (np=npf0, npf, if(F[np,1]%4==1, my (d=decomp2sq(F[np,1]));
     for (k=1, F[np,2], listput(V,d))));
  my(bplus=2^(#V+1));
  for (selbin=0, 2^#V-1, my(P=P0);
    for(k=1, #V, P*=if(digits(selbin+bplus,2)[k+1], V[k], conj(V[k])));
      listput (L, Q*vecsort([abs(real(P)),abs(imag(P))],,4)));
  Set(L)))
};

\\ Asymptotic behavior of R(r)
Rsqasy(r) = {
  (log(r) + Euler + log(8)/2)/Pi
};

\\ Angular correction of asymptotic R
delsq(x,y) = {
  my (r2=x^2+y^2, a=1/(12*Pi*r2), phi=atan(y/x)-Pi/8); a*sin(4*phi)
};

\\ main function, input is target resistance in ohms,
\\ result is vector of grid point coordinates [x,y]
Rsqlatt(R) = {
  my(r2,r2list=vector(5),nlist=0,deltaR=oo,Dbest=[oo,oo]);
\\ Special cases, i.e., very small R, handled explicitly
  if (R < 1/4, return([0,0]));
  if (R < (1/Pi+1/4), return([1,0]));
  if (R < 0.6816901138, return([1,1]));
\\ use asymptotic formula to determine closest squared distance
  r2 = round (exp(Pi*R-Euler-log(8)/2)^2);
\\ List of candidate squared radii
  if (issum2sq(r2), r2list[1]=r2;nlist=1);
  for (dr2=1, oo, my(r2m=r2-dr2,r2p=r2+dr2);
    if (issum2sq(r2m), nlist++; if(nlist>#r2list||r2m<2, break); r2list[nlist]=r2m);
    if (issum2sq(r2p), nlist++; if(nlist>#r2list, break); r2list[nlist]=r2p));
  if (nlist>#r2list, nlist--);
\\ Find best match
\\ outer loop over candidate radii
  for (nr=1, nlist, my(r=sqrt(r2list[nr]), D=d2sq(r2list[nr]), d=#D, Rras=Rsqasy(r));
\\ inner loop over grid points on circle with candidate radius
    for (nd=1, d, my (Rapprox=Rras+delsq(D[nd][1],D[nd][2]), dR=abs(R-Rapprox));
     if (dR < deltaR, Dbest=D[nd]; deltaR=dR)
  ));
  Dbest
};