\\ Repair of strip bijection for 2 overlaid square grids with common origin,
\\ second grid rotated by Pi/4. Terms of A307110 are partially replaced.
\\ Author: Hugo Pfoertner, 2019-05-22
\\ Included for convenience: 17992 terms of A305575 and A305576
\\ providing an enumeration of points of square grid with primary sorting
\\ by radius and secondary sorting by polar angle
atan2(y, x)=if(x>0, atan(y/x), if(x==0, if(y>0, Pi/2, -Pi/2), if(y>=0, atan(y/x)+Pi, atan(y/x)-Pi)));
angle(x, y)=(atan2(y, x)+2*Pi)%(2*Pi);
{a004018(n) = if( n<1, n==0, 4 * sumdiv( n, d, (d%4==1) - (d%4==3)))};
a305575=vector(17992);
a305576=vector(17992);
n=0;
for(s=1,5725, my(r=a004018(s)); if (r>0, my(v=matrix(r,3), w=vector(r), m=sqrtint(s), L=0);\
   for(i=-m, m, my(k=s-i^2, kk);\
   if(k==0,\
    v[L++,1]=i; v[L,2]=0; v[L,3]=angle(i,0)\
   ,\
   if(issquare(k), kk=sqrtint(k);\
     forstep(j=-kk, kk, kk+kk, v[L++,1]=i; v[L,2]=j; v[L,3]=angle(i,j)))));\
   my(p=vecsort(v[,3], ,1));\
   for(k=1, L, a305575[n++]=v[p[k],1];a305576[n]=v[p[k],2];);));
\\ Terms now on arrays a305575 and a305576
\\
\\ Function to determine position of a grid point
findinring(i,j)={my(s=i*i+j*j);
  if(s==0,
    return(0),
    forstep(k=floor(Pi*(s+1))+sqrtint(s), 1, -1,
            if(i==a305575[k]&&j==a305576[k], return(k))))};
\\
\\ also included: 11000 terms of A307110 (unmodified strip bijection)
\\ (11000 because repair will address some points beyond target range of 10000)
a307110=vector(11000);
\\ Klaus Nagel's strip bijection
biject(i,j)=
   {my(C=cos(Pi/8), S=sin(Pi/8), T=S/C,
       gx=i*C-j*S,
       gy=i*S+j*C,
       k, xm, ym, v=[0,0]);
   k=round(gy/C);
   ym=C*k;
   xm=gx+(gy-ym)*T;
   v[1]=round((xm-ym*T)*C);
   v[2]=round((ym+v[1]*S)/C); v
   };
\\ apply to 11000 points given by corresponding locations in rings
for(n=1, 11000, my(v=biject(a305575[n],a305576[n]));a307110[n]=findinring(v[1],v[2]));
\\ A307731 inherits most of its terms from A307110
 a307731=a307110;
\\ Functions for rotations by -Pi/4
xr(x,y)=(x+y)/sqrt(2);
yr(x,y)=(y-x)/sqrt(2);
\\ Function to determine orientation of repair within grid
repair(i, k)={my(v=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]);
/* Index offsets of the 6 possible orientations of the repair neighborhood
  relative to the position of the point in grid G,
  for which the strip bijection of A307110 yields a distance to a bijection
  partner in grid H exceeding S2=sqrt(5)*sin(Pi/8) */
 if(i==0&&k==-1, v=[0, 0, 0, -1, 1, 0, 0, 0, 1, -1, 1, 0, 0, -1, 1, -1]);
 if(i==0&&k==1, v=[0, 0, 0, 1, -1, 0, 0, 0, -1, 1, -1, 0, 0, 1, -1, 1]);
 if(i==1&&k==0, v=[0, 0, -1, 0, 0, -1, 0, 0, -1, -1, 0, -1, -1, 0, -1, -1]);
 if(i==1&&k==1, v=[0, 0, 0, -1, 0, 1, 0, 0, -1, 1, -1, 0, -1, 0, -1, -1]);
 if(i==-1&&k==0, v=[0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1]);
 if(i==-1&&k==-1, v=[0, 0, 0, 1, 0, -1, 0, 0, 1, -1, 1, 0, 1, 0, 1, 1]); v};
\\ Main loop
{for(j=1, #a307110,
 my(xg=a305575[j],
    yg=a305576[j],
    kh=a305575[a307110[j]],
    mh=a305576[a307110[j]],
    xh, yh, xm, ym, dgrid, xo, yo, w=[0, 0, 0, 0]);
\\ Rotated positions in grid H
  xh=xr(kh, mh); yh=yr(kh, mh);
\\ Check if distance between bijection partners (xg, yg)<->(xh, yh) exceeds S2.
  if(sqrt((xh-xg)^2+(yh-yg)^2)>sqrt(5)*sin(Pi/8),
\\ Midpoint of line connecting bijection partners
  xm=(xg+xh)/2; ym=(yg+yh)/2;
\\ Find two points, one on each side,
\\ on the perpendicular bisector of the line joining the bisection partners
\\ at distance S1=cos(Pi/8) from the line's midpoint (xm, ym).
  dgrid=10.0;
  forstep(is=-1, 1, 2, my(d, xc, yc, ss, S1=cos(Pi/8));
\\ Special case yg=yh should never occur
  if(yg==yh,
       xc=xm; yc=ym+is*S1
    ,
       ss=(xg-xh)/(yh-yg);
       xc=xm+is*S1/sqrt(1+ss^2);
       yc=ym+ss*(xc-xm)
    );
\\ For the point that is closer to a grid point,
\\ the distances, rounded to the nearest integer,
\\ to (xg, yg) and (xh, yh) are determined.
  d=abs(round(xc)-xc)+abs(round(yc)-yc);
  if(d<dgrid, dgrid=d; xo=xc; yo=yc));
\\ Only 6 combinations of +-1 and 0 can be found as rounded distances.
\\ Feed x-distances into repair function. The y-distances are redundant.
\\ Its output is a vector of length 16 providing the offsets
\\ in grid coordinates of the revised junctions
\\ between the 4 pairs of bijection partners.
  w=repair(round(xo)-xg, round(-yr(xo, yo))-round(-yr(xh, yh)));
\\ Rotate back by "misusing" functions xr and yr
  xo=round(-yr(xh, yh)); yo=round(xr(xh, yh));
\\ convert 2d grid point coordinates
\\ of the affected 4 pairs back to 1d enumeration
\\ given by A305575 and A305576
\\ and overwrite corresponding positions in array a307731 which
\\ had been initialized as copy of a307110.
  for(i=0, 3, a307731[findinring(xg+w[4*i+1],yg+w[4*i+2])]
                    = findinring(xo+w[4*i+3],yo+w[4*i+4]))));}
\\ vector a307731 with applied repairs is now available
\\ De-activate comments to write b-file
\\ write("b307731.txt", 0, " ", 0);
\\ for(j=1, 10000, write("b307731.txt", j, " ", a307731[j]));
\\ Test: Print affected positions with index < 1000
for (k=1,1000,if(a307731[k]!=a307110[k],print1(k,", ")));
\\ Expected result
\\ 124, 128, 140, 144, 172, 176, 180, 184, 235, 239, 253, 257, 261, 266, 269, 274,
\\ 279, 283, 309, 315, 317, 321, 325, 326, 329, 330, 334, 338, 342, 346, 366, 370,
\\ 387, 388, 395, 396, 408, 416, 422, 426,