Global Functions  for TMJ article

The section where the function was originally discussed is in parenthesis.

Disclaimer:  These functions were coded by Barry H. Dayton in Mathematica 11.3 as simple illustrations of algorithms in  A Wolfram Language Approach to Real Numerical Algebraic Plane Curves  only.  They are not intended for general use.

This version of Global Functions contains examples of each function.  You can check the integrity of these functions on your Mathematica version by running this sheet with Evaluate Notebook. There should be no error messages produced.  Otherwise you should just run Evaluate Initial Cells before using these functions.

Index to Functions

({{aCurve[pts,x,y], 2.5`}, {affineInflectionPoints[f,x,y], 7.4`}, {allInfectionPoints[f,x,y], 7.4`}, {asymptotes[f,p,x,y], 7.1`}, {axisPoints[f,x,y], 9.3`}, {closestPoint[f,q,x,y], 3.5`}, {coefficientSigns[f,X], 9.2.1}, {complexSingularPoints[f,u,v]
complexProjectiveInfinitePoints, 5.3.2
6.4.1}, {complexProjectiveSingularPoints[h,u,v,tol], 6.4.1}, {conicWithFoci[p1,p2,q,x,y,options], 7.2`}, {criticalPoints[f,x,y], 3.5`}, {cTransform[f,p,x,y], 7.3`}, {dAss[f,pt,n], A.4}, {dDiagram[A, options], 9.2.2, GF9A}, {ddVerts[d], GF9A}, {descartesDiagram[sa], 9.2.1, GF9A}, {descartesSketch[B,d], 9.2.1, GF9A}, {deTrace[f,p,x,y,s,r], 4.6`}, {difOp[f,pt,i,j], A.4}, {dvAssoc[A], 9.22, GF9A}, {endFinder[f,n,CP,x,y], 4.6`}, {fAssoc[f,x,y], A.3}, {flt[{x,y},A], 6.2`}, {FLT[f,A,x,y], 6.2`}, {flti[p,A], 6.2`}, {FLTh[h,A,x,y,z], 6.2`}, {flth[{x,y,z},A], 6.2`}, {flts[{x,y},A], 6.2`}, {form[h,k,x,y], 3.3`}, {gaussAssoc[f,rr,x,y], 9.5`}, {gaussCurve[p,x,y,z], 1.4`}, {gaussGraph[f,gp,x,y,tol], 4.6`}, {gaussGraphNS[f,CP,x,y,tol], 4.6`}, {gaussNewton2[f,g,p0,x,y], 2.3`}, {gaussNewton3[f,g,h,p0,x,y], A.2}, {gaussPlot[f,r,pts], 4.2`}, {ghsMap[p,r], 8.3`}, {gPoints[f,r,x,y], 4.3`}, {graphSD[G,H], 8.5`}, {gTvec[f,p,x,y], 4.3`}, {homog[f,x,y,z], 5.2`}, {homothety[r,s], 6.3.2}, {infiniteRealPoints[f,x,y], 3.3`}, {inspectInfPoint[f,ip,r,x,y,options], 9.3`}, {intersectionMultiplicity[f,g,p,option], 5.3.1, A.4}, {inverseGHSmap[{t1,t2},r], 8.3`}, {invSyl[row,d,x,y], A.3}, {ip2z[ip,x,y], 6.3.3}, {iTransform[l,x,y], 6.3.3}, {klRotation[k,l,x,y], 6.3.1}, {kReflection[k,x,y], 6.3.1}, {kShear[k,x,y], 6.3.2}, {label13[+,+,txt], 9.2.2}, {label24[-,+,txt], 9.2.2}, {labeli[a,b,txt] i=I,II,III,IV, 9.2.2}, {line[p,q,x,y], 1.1`}, {macaulayD1[f,p,m], A.4}, {macaulayMatrix[f,g,p,m], A.4}, {maxForm[h,x,y], 3.3`}, {moebiusB[t,s], 8.7`}, {moebiusPhi[p], 8.7`}, {moebiusPLMB (constant), 8.7`}, {moebiusRMPP (constant), 8.7`}, {nDivide[h,g,x,y,tol], A.3}, {nearPoints[f,r,x,y,eps], 4.6`}, {newtonHyperbola[n,x,y], 1.3, GF1.3}, {nGCD[f,h,x,y,tol], A.3}, {p2aEq[u,x,y], A.5}, {p2aMatrix[F,d,t,x,y], A.5}, {p2cTransform (constant), 7.3.2}, {pathFinder[f,p,q,s,x,y,option], 4.4`}, {pathFinderDE[f,p,q,r,x,y,eps], 4.6`}, {pathFinderT[f,p,q,s,x,y,option], GF4.4}, {pExps[d], 2.5`}, {pickCyan[a,b], 9.2.2}, {pickMagenta[a,b], 9.2.2}, {pickSingularity[a,b], 9.2.2}, {plotAssoc[A,c1,c2], GF9A}, {plotDiamonds[f,rr,x,y], 9.5`}, {plotPath[f,G,V,x,y,s], 4.5`}, {plotPathDE[f,G,V,x,y,s], 4.6`}, {plotPathDES[f,G,V,x,y,eps], 4.6`}, {plotPathsS[f,G,V,x,y,s], 4.6`}, {plotPathsT[f,G,V,x,y,s], 4.5`}, {pointMinForm[f,p,x,y], 3.4`}, {pqTranslation[p,q], 6.3.1}, {randomIntegerPolynomial[d,kk,x,y], 2.5`}, {randomRealPolynomial[d,x,y], 2.5`}, {realRegularPoints[f,x,y,n], 5.3.2}, {realSingularPoints[f,x,y], 3.2`}, {sDiagram[A,options], GF9A}, {shiftFA[FA,q,n], A.3}, {singularFactor[f,x,y,tol], A.6}, {singularFactorInf[f,x,y,tol], A.6}, {sqAssoc (constant), GF9A}, {sqFree[f,x,y,tol], A.3}, {sqrSeq[p,A], GF9A}, {sylD1[f,m,x,y], A.3}, {sylvesterMatrix[f,g,m,x,y], A.3}, {t*Assoc (constant, *=1..4), GF9A}, {tangentRealPoints[f,k,x,y], 3.1`}, {tDeg[f,x,y], 2.5`}, {tLine[f,p,x,y], 3.1`}, {viroAssociation[f,x,y], 9.2.2}, {viroDiagram[f,x,y], 9.2.2}, {weierstrassNormalForm[f,ifp,x,y], 7.5`}, {zeroA[d], GF9A}})

The functions listed by the section of the book where they are first described.

line[p_,q_,x_,y_]:=Module[
{a,b,c,M},
M={
Switch[Length[p],2,Append[p,1],3,p],
Switch[Length[q],2,Append[q,1],3,q],
RandomReal[{-1,1},3]
};
Normalize[{a,b,c}/.First@Solve[M.{a,b,c}=={0,0,1}]].{x,y,1}
]

line[{2,3},{3,7},x,y]

line[{2,3},{1,4,0},x,y]

Expand[%/Coefficient[%,y]]

newtonHyperbola[n_,x_,y_]:=Module[{v,c,k,n2dig,s},
If[(n>2^28) ||(n<0), Print["n out of range"];Return[Fail]];
v={{0,0},{1,0},{0,1},{2,0},{1,1},{0,2},{3,0},{2,1},{1,2},{0,3},{4,0},{3,1},{2,2},{1,3},{0,4},{5,0},{4,1},{3,2},{2,3},{1,4},{0,5},{6,0},{5,1},{4,2},{3,3},{2,4},{1,5},{0,6}};
c={1,1,1,.2,1,.2,.008,.2,.2,.008,.000064,
.008,.2,.008,.000064,.2^10,.2^6,.008,.008,.2^6,.2^10,.2^15,.2^10,.2^6,.2^3,.2^6,.2^10,.2^15};
k=Which[ n<8,3,n<64,6,n<1024 ,10,n<32768,15,n<2^21,21,True,28];
n2dig=IntegerDigits[n,2,k];
s=Reverse[n2dig];
FromCoefficientRules[Table[v[[i]]->(-1)^s[[i]]*c[[i]], {i,k}],{x,y}]]

newtonHyperbola[213,x,y]

gaussCurve[p_,x_,y_,z_]:=ComplexExpand[Re[p/.{z->x+I*y}]];

gaussCurve[z^2-2z+3,x,y,z]

newton1[f_,p0_,z_]:=Module[{p1,p2},
p1=N[p0];
Chop[Catch[Do[p2=p1-(f/D[f,z])/.{z->p1};
If[Abs[p1-p2]<1.*^-14,Throw[p2]];
p1=p2,{30}]]]
]

newton1[x^3-3 x-1, 4.712,x]

gaussNewton2[f_,g_,p0_,x_,y_]:=Module[{Q,Q0,p1,p2,J0,n},
n=20;
p1=N[p0];
Catch[Do[Q={f,g}/.Thread[{x,y}->p1];
J0=D[{f,g},{{x,y}}]/.Thread[{x,y}->p1];
p2=p1-Inverse[J0].Q;
If[Norm[p1-p2]<1.*^-13,Throw[p2]];
p1=p2,{n}]]]

gaussNewton2[x-y+x y, x^2-y^2-3,{2,3},x,y]

tDeg[f_,x_,y_]:=Max[Table[Total[t],{t,Keys[CoefficientRules[f,{x,y}]]}]];

tDeg[x+y^2+x y^3, x,y]

pExps[d_]:=Module[{pex},
pex[0]={{0,0}};
pex[n_]:=pex[n]=Join[pex[n-1],Table[{n-j,j},{j,0,n}]];
pex[d]]

pExps[3]

randomRealPolynomial[d_,x_,y_]:=FromCoefficientRules[Table[e->RandomReal[{-1,1}/(4^Total[e])],{e,pExps[d]}],{x,y}]

randomRealPolynomial[3,x,y]

randomIntegerPolynomial[d_,kk_,x_,y_]:=Module[{k,S,b,j,m,t},
b=Binomial[d+2,2];
If[kk<1 || kk>b, k=b,k=kk];
S=RandomSample[pExps[d],k];
m=Max[Table[Total[t],{t,S}]];
If[m<d,j=RandomInteger[{0,d}]; S=Append[Take[S,k-1],{j,d-j}]];FromCoefficientRules[Table[e->RandomInteger[{-9,9}],{e,S}],{x,y}]]

randomIntegerPolynomial[3,4,x,y]

aCurve[pts_,x_,y_]:=Module[{d,P,M,B,n,c,pow},
pow[a_,n_]:=If[n==0,1,a^n];
d=Switch[Length[pts],2,1,5,2,9,3,14,4,20,5,27,6,_,Return["number of points must be 2,5,9,14,20,27"]];
P=pExps[d];
n=Length[P];
M=Table[If[Length[p]==2,pow[p[[1]],e[[1]]]*pow[p[[2]],e[[2]]],pow[p[[1]],e[[1]]]*pow[p[[2]],e[[2]]]*pow[p[[3]],d-e[[1]]-e[[2]]]],{p,pts},{e,P}];
AppendTo[M,RandomReal[{-1,1},n]];
B=Append[Table[0,{n-1}],1];
c=LinearSolve[M,B];
FromCoefficientRules[Table[P[[i]]->c[[i]],{i,n}],{x,y}]
];

f=aCurve[{{1,1},{2,2},{3,5},{4,10}, {5,17}},x,y];
Expand[f/-Coefficient[f,y]]

tLine[f_,p_,x,y]:=line[p,{D[f,y],-D[f,x],0}/.Thread[{x,y}->p],x,y];

g=tLine[f, {3,5},x,y];
Expand[g/-Coefficient[g,y]]

tangentRealPoints[f_,k_,x_,y_]:=Module[{J, p,q,sol}, 
J=Det[D[{f,k},{{x,y}}]];
sol=NSolve[{f,J}, {x,y},Reals];
If[Length[sol]>0,{x,y}/.sol,{}] ];

f
tangentRealPoints[f, 4x-y,x,y]

form, finds degree k form in bivariate polynomial

form[h_,k_,x_,y_]:=FromCoefficientRules[Select[CoefficientRules[h,{x,y}],Total[#[[1]]]==k&],{x,y}];

f=(x-1)^3-(y-2)^2;
form[f,2,x,y]

realSingularPoints[f_,x_,y_]:=Module[{eps},
eps=1.*^-6;
Flatten[Take[Select[Tally[tangentRealPoints[f,RandomReal[{-1,1},2].{x,y},x,y],Norm[#1-#2]<eps&],#[[2]]>1&],All,1],1]];

realSingularPoints[f, x,y]

maxForm, finds the highest degree non-zero  form of bivariate polynomial

maxForm[h_,x_,y_]:=form[h,tDeg[h,x,y],x,y];

f=2-x+y-x^2-x y +x^3-3x^2y
maxForm[f,x,y]

infiniteRealPoints[f_,x_,y_]:=Module[{fm,rnl,sol},
fm=maxForm[f,x,y];
rnl=RandomReal[{-1,1},2].{x,y}-1;
sol=NSolve[{fm,rnl},{x,y},Reals];
If[Length[sol]==0,Return[{}],sol={x,y}/.sol];sol=DeleteDuplicates[sol,Norm[#1-#2]<.001&];Table[Append[sol[[i]],0],{i,Length[sol]}]
];

f
infiniteRealPoints[f,x,y]

pointMinForm, unlike maxForm, the minimum form depends on the point used as reference, this function finds the minimum form at point p

pointMinForm[f_,p_, x_,y_]:=Module[{A,h,j,tol},
tol=1.*^-9;
h=Chop[Expand[f/.Thread[{x,y}->{x,y}+p]],tol];
A=<| CoefficientRules[h,{x,y}] |>;
j=Min[Table[Total[k],{k,Keys[A]}]];
form[h,j,x,y]];

f= Expand[(x-2)^3 (y-3)^2]
pointMinForm[f,{2,5},x,y]

criticalPoints[f_,x_,y_]:=DeleteDuplicates[{x,y}/.NSolve[{f, y D[f,x]-x D[f,y]},{x,y},Reals],Norm[#1-#2]<1.*^-6&];

f=x^4-3 x y+2 y^4-5
criticalPoints[f,x,y]

closestPoint[f_,q_,x_,y_]:=Module[{g,p, cp},
g=Expand[f/.Thread[{x,y}->{x+q[[1]],y+q[[2]]}]];
cp= criticalPoints[g,x,y];
p=MinimalBy[cp,Norm[#]&][[1]];
p+q];

f
closestPoint[f,{1,1},x,y]

gaussPlot[f_,r_,pts_]:=Module[{plta,pltb,pltc},
plta:=RegionPlot[x^2+y^2<r^2,{x,-1.1r,1.1r},{y,-1.1r,1.1r}, PlotStyle->LightOrange, BoundaryStyle->Orange];
pltb:=RegionPlot[{f>0&& x^2+y^2<r^2 },{x,-r,r},{y,-r,r},BoundaryStyle->Blue];
pltc:=ListPlot[Table[Labeled[pts[[n]],n],{n,Length[pts]}],PlotStyle->Red, PlotMarkers->Automatic];
Show[plta,pltb,pltc]];

g=x^3-x+1-y^2;
bp={x,y}/.NSolve[{x^2+y^2-16,g},{x,y},Reals]
gaussPlot[g,4,bp]

gTvec[f_,p_,x_,y_]:=Cross[Normalize[Chop[{D[f,x],D[f,y]}/.Thread[{x,y}->p],1.*^-6]]]

gTvec[g,bp[[2]],x,y]

gPoints[f_,r_,x_,y_]:= Module[{bp,cp,sol},
sol=NSolve[{f,x^2+y^2-r^2},{x,y}, Reals];
If[Length[sol]>0,bp={x,y}/.sol,bp={}];
bp=Reverse[SortBy[bp,ArcTan@@#&]];
cp= DeleteDuplicates[ Select[{x,y}/.NSolve[{f,y D[f,x]-x D[f,y]},{x,y},Reals],Norm[#]<r&],Norm[#1-#2]<.1&];
Join[bp,cp]];

gPoints[g,4,x,y]

Options[pathFinder]={maxit->30};
pathFinder[f_,p_, q_,s_,x_,y_,OptionsPattern[]]:= Module[{L,p0,p1,u, mxt,k},
mxt=OptionValue[maxit];
p0=p;
L=Reap[Sow[p];
k=0;
While[Norm[q-p0]>2s && k<mxt,
u=p0+s*gTvec[f,p0,x,y];
p1 =closestPoint[f,u,x,y];
Sow[p1];
p0=p1;
k++];
Sow[q]][[2,1]];
If[k≥ mxt, Print[Style["Warning: iteration limit reached at",Orange],p1]];
L];

g=x^3-x+1-y^2;
pathFinder[g,{-1,-1},{-1,1},.4,x,y]

(4.5) PathFinderT, tracing using tangent vector given by gTvec.  Same arguments as path finder but will be faster but may work less well in some cases.

Options[pathFinderT]:={maxit->20};
pathFinderT[f_,p_, q_,s_,x_,y_,OptionsPattern[]]:= Module[{L,p0,p1,tv,u,l, maxit2,k},
maxit2=OptionValue[maxit];
p0=p;
L=Reap[Sow[p];
k=0;
While[Norm[q-p0]>2s && k<maxit2,
tv=gTvec[f,p0,x,y];
u=p0+s*tv;
l=line[u,{tv[[2]],-tv[[1]],0},x,y];
p1={x,y}/.FindRoot[{f,l},{x,u[[1]]},{y,u[[2]]}];
Sow[p1];
p0=p1;
k++];
Sow[q]][[2,1]];
If[k≥ maxit2, Print[Style["Warning: iteration limit reached at",Orange],p1]];
L];

g
pathFinder[g,{-1,-1},{-1,1},.1,x,y]
pathFinder[g,{-1,-1},{-1,1},.1,x,y,maxit->10]

(4.5) Plot Path, f curve, G graph, V actual points of curve used as vertices, s step size

plotPath[f_,G_,V_,x_,y_,s_]:=Module[{n,L},
n=Length[G];
L=Flatten[Table[pathFinder[f,V[[G[[i,1]]]],V[[G[[i,2]]]],s,x,y],{i,n}],1]];

g
G={1->2,2->3,3->4};
V={{0,-1},{-1,-1},{-1,1},{0,1}};
Timing[plotPath[g,G,V,x,y,.3]]

In many cases we can use the faster version

plotPathT[f_,G_,V_,x_,y_,s_]:=Module[{n,L},
n=Length[G];
L=Flatten[Table[pathFinderT[f,V[[G[[i,1]]]],V[[G[[i,2]]]],s,x,y],{i,n}],1]];

g
G={1->2,2->3,3->4};
V={{0,-1},{-1,-1},{-1,1},{0,1}};
Timing[plotPathT[g,G,V,x,y,.3]]

(4.6) Curve tracing using NDSolveValue.  f curve, p  initial point, s maximal distance from initial point, r maximum distance from origin.  Returns pair of interpolating functions, i.e x,y  parameters of curve.  Interpolating function approximates curve, does not necessarily give actual points on curve.

deTrace[f_,p_,x_,y_,s_,r_]:=Module[{X,Y,F,u,t,tv,ttv,XS,YS},
F=f/.Thread[{x,y}->{X[t],Y[t]}];
tv=gTvec[f,p,x,y];
u=Normalize[RandomReal[{.8,1.2},2]];
ttv=u.tv;
NDSolveValue[{u.{X'[t],Y'[t]}==ttv,D[F,t]==0, {X[0],Y[0]}==p,
WhenEvent[X[t]^2+Y[t]^2>r^2||u[[1]]*X'[t]+u[[2]]*Y'[t]==0||(X[t]-p[[1]])^2+(Y[t]-p[[2]])^2>s^2,"StopIntegration"]},{X,Y},{t,0,20}]
]

Note : NDSolveValue::ndsz is not an error message but an information message which may or may not appear here depending on random value used.

g
{X,Y}=deTrace[g,{-1,-1},x,y,2,3]

(4.6) End point of interpolating function pair from detrace.

deEnd[{X_,Y_}]:=Module[{t1},
t1=X[[1,1,2]];
{{X[t1],Y[t1]},t1}]

deEnd[{X,Y}]

pathFinderDE[f_,p_,q_,r_,x_,y_,eps_]:=Module[{i,e,X,Y,s},
Off[NDSolveValue::ndsz];
e=p;
s=Norm[p-q];
i=1;
Reap[Sow[p];
While[s>eps && Norm[e]<2*r&& i<60,
{X,Y}=deTrace[f,e,x,y,s*eps,r];
e=deEnd[{X,Y}][[1]];
Sow[e];
s=Norm[e-q];
i++];
Sow[q]][[2,1]]];

pathFinderDE[g,{-1,-1},{-1,1},3,x,y,.5]

plotPathDE[f_,G_,V_,x_,y_,s_]:=Module[{r,n,L,LL,PF},
n=Length[G];
r=Max[Norm/@V];
L=Flatten[Table[pathFinderDE[f,V[[G[[i,1]]]],V[[G[[i,2]]]],r,x,y,s],{i,n}],1];
Append[L,V[[G[[n,2]]]]]]

Timing[plotPathDE[g,G,V,x,y,.5]]

endFinder[f_,n_,CP_,x_,y_,eps_]:=Module[{r,i,k,e,X,Y,s,CPD},
Off[NDSolveValue::ndsz];
r=Max[Norm/@CP];
e=CP[[n]];
CPD=CP;
CPD[[n]]=Nothing;
s=Norm[e-First[Nearest[CPD,e]]];
i=1;
e=Catch[
While[s>eps && Norm[e]<1.01*r&& i<60,
{X,Y}=deTrace[f,e,x,y,s,r];
e=deEnd[{X,Y}][[1]];
If[Norm[e]> r, Throw[e]];
s=Norm[e-First[Nearest[CPD,e]]];
i++];
Throw[e]];
If[Norm[e]>r,Return[0]];
k=1;
While[k≤Length[CP] && Norm[e-CP[[k]]]>s,k++];
If[k<=Length[CP],k,Fail]];

endFinder[g,2,V,x,y,.4]

nearPoints[f_,r_,x_,y_,pt_,ep_]:=Module[{np,stnpts},
np={x,y}/.NSolve[{f,(x-pt[[1]])^2+(y-pt[[2]])^2-ep^2},{x,y},Reals];
np=Select[np,Norm[#-pt+ep/2*Normalize[gTvec[f,#,x,y]]]> ep&];
Select[np,Norm[#]<r&]];

f=(-1+x)^3-(-2+y)^2;
nearPoints[f,4,x,y,{1,2},.4]

gaussGraphNS[f_,CP_,x_,y_,eps_]:=Module[{n,G},
n=Length[CP];
G=Select[Table[k->endFinder[f,k,CP,x,y,eps],
{k,n}],#[[2]]>0&];
Graph[G,VertexLabels->"Name",ImageSize->Small]];

g
gaussGraphNS[g,V,x,y,.4]

gaussGraph[f_,gp_,x_,y_,eps_]:=Module[{n,k,e,G,np,CP,r,S},
CP=gp;
r=Max[Norm[#]&/@gp];
n=Length[CP];
G={};
k=1;
S=Reap[While[k≤n, 
If[Norm[gTvec[f,CP[[k]],x,y]]<eps/2, 
np=nearPoints[f,r,x,y,CP[[k]],1.25*eps];
Do[Sow[k->n+j],{j,Length[np]}];
CP=Join[CP,np];
n=Length[CP];
e=0,
e=endFinder[f,k,CP,x,y,eps]];
If[e >0,Sow[k->e]];
k++]];
G=S[[2,1]];
Echo[{Show[ContourPlot[f==0, {x,-1.2r,1.2r},{y,-1.2r,1.2r}],ListPlot[Table[Labeled[CP[[j]],j],{j,Length[CP]}]]],
Graph[G,VertexLabels->"Name"]}];
{CP,G}]

f
Vf={{2,1},{1,2},{2,3}};
gaussGraph[f,Vf,x,y,.4]

plotPathS[f_,G_,V_,x_,y_,s_]:=Module[{n, L,r},
n=Length[G];
r=Max[Norm/@V];
L=Flatten[Table[If[Norm[gTvec[f,V[[G[[i,1]]]],x,y]]<s/2, {V[[G[[i,1]]]],V[[G[[i,2]]]]},
pathFinder[f,V[[G[[i,1]]]],V[[G[[i,2]]]],s,x,y]],
{i,n}],1];
AppendTo[L,V[[G[[n,2]]]]];
L];

plotPathDES[f_,G_,V_,x_,y_,eps_]:=Module[{r,n, L},
n=Length[G];
r=Max[Norm[#]&/@V];
L=Flatten[Table[If[Norm[gTvec[f,V[[G[[i,1]]]],x,y]]<eps/2, {V[[G[[i,1]]]],V[[G[[i,2]]]]},
pathFinderDE[f,V[[G[[i,1]]]],V[[G[[i,2]]]],r,x,y,eps]],
{i,n}],1];
AppendTo[L,V[[G[[n,2]]]]];
L];

Vsf={{2,1},{1,2},{1.4196433776070834`,2.271844506346037`},{2,3}};
Gsf={1->2,2->3,3->4};
Timing[plotPathS[f,Gsf,Vsf,x,y,.3]]
Timing[plotPathDES[f,Gsf,Vsf,x,y,.3]]

pLine[p_,q_,x_,y_,z_]:=Module[{a,b,c,M},
M={p,q,RandomReal[{-1,1},3]};
({a,b,c}/.(Solve[M.{a,b,c}=={0,0,1}][[1]])).{x,y,z}
];

pl=pLine[{2,3,4},{-1,3,0}, x,y,z]
Expand[pl/Coefficient[pl,x]]

pLine[{2+I,3,1},{1,I,0},x,y,z]

homog[f_,x_,y_,z_]:=Expand[z^tDeg[f,x,y](f/.Thread[{x->x/z,y->y/z}])];

homog[x y^2-x^2+y^2-xy-x+1,x,y,z]

dTol:=1.*^-12;

f=x-y+x^3;
g=x-y+y^3;
intersectionMultiplicity[f,g,{0,0},dTol]

complexSingularPoints[f_,u_,v_]:=Module[{a,b,sol1,sol2,solz},
{a,b}=RandomComplex[{-1-I,1+I},2];
sol1={u,v}/.NSolve[{f,a D[f,u]+b D[f,v]}];
{a,b}=RandomComplex[{-1-I,1+I},2];
sol2={u,v}/.NSolve[{f,a D[f,u]+b D[f,v]}];
solz=Join[sol1,sol2];
Flatten[Take[Select[Tally[solz,Norm[#1-#2]<.001&],#[[2]]>1&],All,1],1]]

h=9-10 x^2+x^4+6 y^2+2 x^2 y^2+y^4;
complexSingularPoints[h,x,y]

realRegularPoints[f_,x_,y_,n_]:=Module[{a,b,c,d,k},
If[tDeg[f,x,y]<2, Echo["Degree must be > 1"]; Return[{}]];
k=0;
              Reap[While[k<n, 
{a,b}=RandomReal[{-10,10},2];
{c,d}=closestPoint[f,{a,b},x,y];
If[Abs[Det[D[{f,a x-b y},{{x,y}}]]/.Thread[{x,y}->{c,d}]]>.01,Sow[{c,d}];k++]]][[2,1]]];

f=x^4-4 x y+y^4;
realRegularPoints[f,x,y,4]

flt[{x_,y_},A_]:=Module[{u,v,w},
{u,v,w}= A.{x,y,1};
If[Abs[w]<1.*^-6,{},{u/w,v/w}]];

flti[p_,A_] :=Module[{u,v,w},
If[Length[p]==2, {u,v,w}= A.Append[p,1], {u,v,w}= A.p];
If[Abs[w]<1.*^-6, {u,v,0}, {u/w,v/w}]];

flts[{x_,y_},A_]:=Module[{u,v,w},
{u,v,w}= A.{x,y,1};
{u/w,v/w}];

FLT[f_,A_,x_,y_]:=Module[{fh,gh, AI,B,d,z},
d=tDeg[f,x,y];
fh=Expand[z^d(f/.{x->x/z,y->y/z})];
AI=Inverse[A];
B=AI.{x,y,z};
gh=Expand[fh/.Thread[{x,y,z}->B]];
Chop[gh/.{z->1},1*^-10]
];

flth[{x_,y_,z_},A]:=A.{x,y,z};

 FLTh[h_,A_,x_,y_,z_]:=Expand[h/.Thread[{x,y,z}->Inverse[A].{x,y,z}]];

A={{-2,2,7},{-1,-6,0},{3,3,8}};
flt[{2,1},A]
flt[{-3,1/3},A]
flti[{-3,1/3},A]
flts[{x,y},A]
FLT[y-x^2,N[A],x,y]
flth[{-9,1,3},A]
FLTh[y z-x^2,A,x,y,z]

pqTranslation[p_,q_]:={{1,0,q[[1]]-p[[1]]},{0,1,q[[2]]-p[[2]]},{0,0,1}};

flt[{0,1},pqTranslation[{0,1},{2,4}]]

klRotation[k_,l_,x_,y_]:= Module[{p,a1,a2,b1,b2,T,R, A,B,m},
a1=Coefficient[k,x];
b1=Coefficient[k,y];
a2=Coefficient[l,x];
b2=Coefficient[l,y];
If[Abs[a1*b2-a2*b1]< 1.*^-9,Print["No affine center"];
Return[Fail]];
p={x,y}/.NSolve[{k,l}][[1]];
T={{1,0,-p[[1]]},{0,1,-p[[2]]},{0,0,1}};
R=N[RotationMatrix[{{a1,b1},{a2,b2}}]];
R=Join[Append[R,{0,0}],{{0},{0},{1}},2];
Inverse[T].R.T
]

flt[{1,0},klRotation[x,y,x,y]]

kReflection[k_,x_,y_]:=Module[{a,b,c,T,RX,RO},
RX={{1,0,0},{0,-1,0},{0,0,1}};
a=Coefficient[k,x];
b=Coefficient[k,y];
c=k/.Thread[{x,y}->{0,0}];
If[a==0,
T=pqTranslation[{0,0},{0,-c/b}];
Return[T.RX.Inverse[T]],
RO=klRotation[y,k,x,y];
Return[RO.RX.Inverse[RO]];
]
];

flt[{1,0},kReflection[x,x,y]]

kShear[k_,x_,y_]:=Module[{a,b,c,d,e},
c=Coefficient[k,x];
d=Coefficient[k,y];
e=k/.{x->0,y->0};
a=-e/d;
b=-c/d;
{{1,0,0},{-b,1,-a},{0,0,1}}
];

FLT[x-y,kShear[x-y,x,y],x,y]

homothety[r_, s_] := {{r, 0, 0}, {0, s, 0},  {0, 0, 1}};

flt[{1,1},homothety[2,-3]]

ip2z[ip_,x_,y_]:=Module[ {p,ral,A,B},
If[Length[ip] ≠ 3,Print["not infinite point"];Abort[],If[Abs[ip[[3]]]>1.*^-6,Print["not infinite point"];Abort[],p=Chop[ip,1.*^-6]]];
B={{1.414213562373095`,0.`,0.`},{0.`,0.`,1.414213562373095`},{0.`,-0.7071067811865475`,-0.7071067811865475`}};
If[Abs[p[[1]]]<1.*^-6,Return[B]];
ral=line[p,{0,0},x,y];
B.klRotation[ral,x,x,y]]

flti[{0,1,0},ip2z[{0,1,0},x,y]]

iTransform takes line l to the infinite line moving the original infinite line to the x-axis.  This uses a fixed “random” vector r1.  Essentially iTransform deletes the line l from the projective plane forming a new affine plane.  The function flti will transform points from the original affine plane to the new one.

iTransform[l_,x_,y_]:= Module[{ph,r1},
ph=Normalize[N[{Coefficient[l,x],Coefficient[l,y], l/.Thread[{x,y}->{0,0}]}]];
r1={-0.7812879011987071`,-0.3168196688992504`,0.9191639531555853`};
{r1,{0,0,1},ph}]

flti[{5,7},iTransform[x-y,x,y]]

flti[{3.2,3.2},iTransform[x-y,x,y]]

flti[{2,5,0}, iTransform[3x-2y+5,x,y]]

complexProjectiveSingularPoints[h_,u_,v_,eps_]:=Module[{f,w,A,a,b,sol1,sol2,solz,sing,sing1,sing2},A={RandomInteger[{-100,100},3],{0,0,RandomInteger[{100,200}]},RandomInteger[{-1,1},3]};
f=FLT[h,A,u,v];
{a,b}=RandomInteger[{-100,100},{2,2}].{1,I};
sol1={u,v}/.NSolve[{f,a D[f,u]+b D[f,v]}];
{a,b}=RandomInteger[{-100,100},{2,2}].{1,I};
sol2={u,v}/.NSolve[{f,a D[f,u]+b D[f,v]}];
solz=Join[sol1,sol2];
sing=Flatten[Take[Select[Tally[solz,Norm[#1-#2]<eps&],#[[2]]>1&],All,1],1];
sing1=Table[{sing[[j,1]],If[Abs[sing[[j,2]]]>eps,sing[[j,2]],0]},{j,Length[sing]}];
sing2=flti[#,Inverse[A]]&/@sing1;
sing2=Table[Which[Length[p]==2,p,Abs[p[[1]]]>eps,p/p[[1]],Abs[p[[2]]]>eps,p/p[[2]],True,{0,0,p[[3]]}],{p,sing2}];
Chop[sing2,eps]];

h=x^2+x^4+6 y^2+2 x^2 y^2+y^4;
complexProjectiveSingularPoints[h,x,y,1.*^-7]

complexProjectiveIntersectionPoints[f_,g_,x_,y_]:=Module[{f1,g1,A,sol, sol2},
A={RandomReal[{-1,1},3],{0,0,RandomReal[{1,2}]},RandomReal[{-1,1},3]};
f1=FLT[f,A,x,y];
g1=FLT[g,A,x,y]; 
sol={x,y}/.NSolve[{f1,g1}];
flti[#,Inverse[A]]&/@sol]

complexProjectiveIntersectionPoints[x,x y-1,x,y]

asymptotes[f_,p_,x_,y_]:=Module[{g,gm,A,sol,rnd,k},
A=ip2z[p,x,y];
g=FLT[f,A,x,y];
gm=pointMinForm[g,{0,0},x,y];
rnd=RandomReal[{-1,1},2].{x,y}-1;
sol=DeleteDuplicates[{x,y}/.NSolve[{gm,rnd},{x,y},Reals],Norm[#1-#2]<.001&];
k=Length[sol];
If[k==0,Return[{}]];
Table[FLT[sol[[i,2]]x-sol[[i,1]]y,Inverse[A],x,y],{i,k}]
];

fa=y(x^2-1)-x^3-2x^2;
ipfa=infiniteRealPoints[fa,x,y]
asymptotes[fa,{0,1,0},x,y]
asymptotes[fa,{1,1,0},x,y]

Options[conicWithFoci]={conic->"ell"};
conicWithFoci[p1_,p2_,q_,x_,y_,OptionsPattern[]]:=Module[{c,p0,l,qc,R,T,M,qq,a,b,aa,bb,cc,g,pp0},
p0=N[(p1+p2)/2];
c=Norm[p1-p0];
l=line[p1,p2,x,y];
qc=closestPoint[l,q,x,y];
If[Norm[qc-q]>c ||Abs[ Norm[p1-qc]+Norm[p2-qc]-2c]>1.*^-10, Echo["Bad Point q "];Abort[]];
If[p1[[1]]==0&& p2[[1]]==0,R=IdentityMatrix[3],R=klRotation[l,x,x,y]];
pp0=flt[p0,R];
T=pqTranslation[pp0,{0,0}];
M=T.R;
qq=flt[q,M];
Switch[StringTake[OptionValue[conic],3],"ell",
sol=NSolve[{c^2-b^2+a^2,qq[[1]]^2*b^2+qq[[2]]^2*a^2-a^2b^2},{a,b},Reals],"hyp",sol=NSolve[{c^2-b^2-a^2,-qq[[1]]^2*b^2+qq[[2]]^2*a^2-a^2b^2},{a,b},Reals],_,Echo[Null,"options ell, hyp only"];Abort[]];
If[Length[sol]==0, Echo[Null,"no solution"]; Abort[]];
{aa,bb}={a,b}/.sol[[1]];
Switch[StringTake[OptionValue[conic],3],"ell",g=x^2/aa^2+y^2/bb^2-1,"hyp",g=-x^2/aa^2+y^2/bb^2-1];
FLT[g,Inverse[M],x,y]
]

p1={1,4};
p2={3,1};
q={2,1};
conicWithFoci[p1,p2,q,x,y, conic->"hyperbola"]

cTransform[f_,p_,x_,y_]:=Module[{fh,ph, nh, t,cs,A,d},
d=tDeg[f,x,y];
fh=Expand[t^d*(f/.Thread[{x,y}-> {x/t,y/t}])];
ph=If[Length[p]==2,N[Append[p,1]],N[p]];
ph=ph/Norm[ph];
nh={D[fh,x],D[fh,y], D[fh,t]}/.Thread[{x,y,t}->ph];
nh=nh/Norm[nh];
cs=Cross[nh,ph];
cs= cs/Norm[cs];
A={cs,ph,nh}];

circ1=x^2+y^2-1;
FLT[circ1,cTransform[circ1,{0,1},x,y],x,y]

p2cTransform:= pqTranslation[{0,-1},{0,0}].ip2z[{0,1,0},x,y];

FLT[y-x^2,p2cTransform,x,y]

affineInflectionPoints[f_,x_,y_]:=Module[{H,fh,d,z,sol},
d=tDeg[f,x,y];
fh=homog[f,x,y,z];
H=Det[D[fh,{{x,y,z},2}]];
sol=NSolve[{fh,H,z-1},{x,y,z},Reals];
If[Length[sol]==0,Return[{}],{x,y}/.sol]
];

nh210=newtonHyperbola[210,x,y];
aip=affineInflectionPoints[nh210,x,y]

allInflectionPoints[f_,x_,y_]:=Module[{H,fh,d,z,l,ips,sol},
fh=homog[f,x,y,z];
H=Det[D[fh,{{x,y,z},2}]];
l=RandomReal[{-10,10},3].{x,y,z};
sol=NSolve[{fh,H,l-1},{x,y,z},Reals];
If[Length[sol]==0,Return[{}]];
ips={x,y,z}/.sol;
ips=DeleteDuplicates[ips,Norm[#1-#2]<1.*^-10&];
Which[Abs[#[[3]]]>1.*^-9,Take[#/#[[3]],2],Abs[#[[2]]]>1.*^-9,#/#[[2]],True,#/#[[1]]]&/@ips
];

f
allInflectionPoints[f,x,y]

weierstrassNormalForm[f_,ip_,x_,y_]:=Module[{B1,B2,B3,B4,B,f1,f2,f3,f4,j,k,c4,c6,cy2},
B1=cTransform[f,ip,x,y];
f1=Chop[FLT[f,B1,x,y], 1.*^-9];
cy2=Coefficient[f1,y^2];
If[Length[cy2]>0,Print["Sorry, this example requires special handling"];Abort[]];
k=Expand[Coefficient[f1,y]/cy2/2+y];
B2={{1,0,0},{Coefficient[k,x],1,k[[1]]},{0,0,1}};
f2=FLT[f1,B2,x,y];
B3=homothety[CubeRoot[-(Coefficient[f2,x^3]/Coefficient[f2,y^2])],1];
f3 = FLT[f2,B3,x,y];
f3=-Expand[f3/Coefficient[f3,y^2]];
B4={{1,0,Coefficient[f3,x^2]/3},{0,1,0},{0,0,1}};
f4=FLT[f3,B4,x,y];
B=B4.B3.B2.B1;
c4=Coefficient[f4,x]/27;
c6=f4[[1]]/54;
j=1728c4^3/(c4^3-c6^2);
{f4,B,j}];

weierstrassNormalForm[nh210,aip[[1]],x,y]

ghsMap[p_,r_]:=If[Length[p]==2,N[6*r/Pi*ArcTan[Norm[p]/6/r*Pi]]*Normalize[p],3*r*N[Normalize[Take[p,2]]]];

ghsMap[{1,1,0},3]

inverseGHSmap[q_,r_]:=Tan[Norm[q]/6/r*Pi]*6*r/Pi*Normalize[q];

inverseGHSmap[{x,y},5]

graphSD[G_,H_]:=With[{GH=GraphIntersection[G,H]},FindCycle[EdgeList[GraphUnion[GraphDifference[G,GH],GraphDifference[H,GH]]]]][[1]];

G={14\[UndirectedEdge]11,11\[UndirectedEdge]13,13\[UndirectedEdge]12,12\[UndirectedEdge]15,15\[UndirectedEdge]9,9\[UndirectedEdge]14};H={18\[UndirectedEdge]11,11\[UndirectedEdge]13,13\[UndirectedEdge]12,12\[UndirectedEdge]19,19\[UndirectedEdge]7,7\[UndirectedEdge]18};
graphSD[G,H]

moebiusB[t_,s_]:=Module[{a=N[Sqrt[2]-1],mb,MB,tt},
mb[tt_]:={Cos[2tt+Pi]/(1+a Cos[tt]), Sin[2tt+Pi]/(1+a Cos[tt]),1/(1+a Cos[tt+Pi/2])};
s*mb[t]+(1-s)*mb[t+Pi]];

moebiusPLMB:=Show[ParametricPlot3D[moebiusB[t,s], {t,-Pi,Pi},{s,0,.5}, Mesh->None,PlotStyle->Opacity[0.8]],ParametricPlot3D[moebiusB[t,.5],{t,0,Pi},PlotStyle->{Thick,Blue,Dashed}],Boxed->False,Axes->False];

Show[moebiusPLMB]

moebiusRMPP:=ListLinePlot[{{{-Pi,0},{Pi,0}},{{-Pi,.5},{Pi,.5}},{{Pi,0},{Pi,.5}},{{-Pi,0},{-Pi,.5}},{{0,0},{0,.5}}},PlotStyle->{Brown,Orange,Green,Green,Cyan},Axes->False,Frame->True,AspectRatio->1/3];

moebiusRMPP

moebiusPhi[p_]:=
If[Length[p]==2&&(p[[1]]>0||p[[2]]≠0 ),Return[{ArcTan@@N[p],N[ArcTan[Norm[p]]/Pi]}],Echo[p,"ambiguous point "]; Return[{}]];

L=moebiusPhi/@{{1,1},{1.4,1},{1.2,1.2},{1,1}}
Show[moebiusPLMB,Graphics3D[{Black,Thick,Line[moebiusB@@#&/@L]}]]

For Chapter 9 we list functions with syntax used in the first sections.  Much actual code is given in section 9.5 and to keep this notebook readable we put the code at the end of this notebook since it does not illustrate curve theory or algorithms.

coefficientSigns[f_, X_]:=Module[{A},
A=Association[CoefficientRules[f,X]];
Association[Table[p->-(Sign[A[p]]-1)/2, {p,Keys[A]}]]];

A=coefficientSigns[x^3+x^2-x+1,{x}]

descartesDiagram[A]

B=A;
AssociateTo[B,{{-3}->1,{-2}->0,{-1}->0}];
descartesSketch[B,3]

x^2-x y+y^2-x+y-1;
A2=coefficientSigns[x^2-x y+y^2-x+y-1,{x,y}]
dvAssoc[A2]

viroAssociation[f_,x_,y_]:=dvAssoc[coefficientSigns[f,{x,y}]];
viroDiagram[f_,x_,y_]:=dDiagram[viroAssociation[f,x,y],pcolr->Green,ncolr->Red];

f=1-x+y+x^2-x y+y^2;
Af=viroAssociation[f,x,y]
viroDiagram[f,x,y]

DDf=dDiagram[Af, pcolr->Orange, ncolr->Brown]

pickCyan[a_,b_]:=Graphics[{White,Rectangle[{a+.1,b+.1},{a+.9,b+.9}],Thick,Blue,Circle[{a+1,b+1},.5,{Pi,3Pi/2}], Circle[{a,b},.5,{0,Pi/2}]}];

Show[DDf,pickCyan[0,-1]]

pickMagenta[a_,b_]:=Graphics[{White,Rectangle[{a+.1,b+.1},{a+.9,b+.9}],Thick,Blue,Circle[{a+1,b},.5,{Pi/2,Pi}],Circle[{a,b+1},.5,{3Pi/2,2Pi}]}];

Show[DDf,pickMagenta[0,-1]]

pickSingularity[a_,b_]:=Graphics[{White,Rectangle[{a+.1,b+.1},{a+.9,b+.9}],Thick,Blue,Line[{{a,b+.5},{a+1,b+.5}}],Line[{{a+.5,b},{a+.5,b+1}}]}]

Show[DDf,pickSingularity[0,-1]]

labelI[a_,b_,txt_]:=ListPlot[{{a+.5,b+.5}},PlotStyle->Directive[Transparent,Thin],PlotLabels->Placed[txt,After], Axes->None]

Show[DDf,labelI[-.2,1.2,"A"]]

labelII[a_,b_,txt_]:=ListPlot[{{a-.5,b+.5}},PlotStyle->Directive[Transparent,Thin],PlotLabels->Placed[txt,Above], Axes->None]

labelIII[a_,b_,txt_]:=ListPlot[{{a-.5,b-.5}},PlotStyle->Directive[Transparent,Thin],PlotLabels->Placed[txt,Below], Axes->None]

Show[DDf,labelIII[-1,.1,"A"]]

labelIV[a_,b_,txt_]:=ListPlot[{{a+.5,b-.5}},PlotStyle->Directive[Transparent,Thin],PlotLabels->Placed[txt,After], Axes->None]

label13[a_,b_,txt_]:=If[a≥0&&b≥0,ListPlot[{{{a+.5,b+.5}},{{-a-.5,-b-.5}}},PlotStyle->Directive[Transparent,Thin],PlotLabels->{Placed[txt,After],Placed[txt,Below]}, Axes->None],Echo["label13 error:must have a,b positive"]];

Show[DDf,label13[0,1,"A"]]

label24[a_,b_,txt_]:=If[a≤0&&b≥0,ListPlot[{{{a-.5,b+.5}},{{-a+.5,-b-.5}}},PlotStyle->Directive[Transparent,Thin],PlotLabels->{Placed[txt,Above],Placed[txt,After]}, Axes->None],Echo["label24error: must have a≤0,b≥0"]];

toggle[A_,K_]:=Module[{B},
B=A;
Do[If[B[p]==0,B[p]=1,B[p]=0],{p,K}];
B]

Af
Af1=toggle[Af,{{0,0},{1,0},{-1,1}}]
{DDf,dDiagram[Af1]}

axisPoints[f_,x_,y_]:=
Quiet[Check[SortBy[{x,y}/.NSolve[{f,x y},{x,y},Reals],Abs[#[[2]]]&],{}]]

axisPoints[1-x+y-x y,x,y]

Options[inspectInfPoint]={output->All};
inspectInfPoint[f_,ip_,r_,x_,y_,OptionsPattern[]]:=Module[{A,g,bp,ap,sol,GP},
A=ip2z[ip,x,y];
g=FLT[f,A,x,y];
sol=NSolve[{g,x^2+y^2-r^2},{x,y},Reals];
If[Length[sol]>0,bp={x,y}/.sol;
ap=flti[#,Inverse[A]]&/@bp,
bp={}; ap={}];
GP:=Show[gaussPlot[g,r,bp], ContourPlot[y==0,{x,-r,r},{y,-r,r},ContourStyle->Directive[Orange,Dashed]], ImageSize->Small];
Switch[OptionValue[output],"plot",Return[GP],"data",Return[{A,g,bp,ap}],Default,Continue];
Echo[A,"Transformation"];
Echo[g," Transformed equation"];
Echo[bp,"boundary points coordinates relative to infinite point"];
Echo[ap, "affine coordinates boundry points"];
GP]

inspectInfPoint[y(x^2-1)-1,{0,1,0},1,x,y]

inspectInfPoint[y(x^2-1)-1,{0,1,0},1,x,y,output->"plot"]

{A,g,bp,ap}=inspectInfPoint[y(x^2-1)-1,{0,1,0},1,x,y,output->"data"]

plotDiamonds[f_,rr_, x_,y_]:=Module[{m},
m=1.1*Max[rr];
If[Length[rr]≠tDeg[f,x,y], Echo[" wrong number of diamonds","WARNING"]];
Show[RegionPlot[f>0,{x,-m,m},{y,-m,m}],ContourPlot[Abs[x]+Abs[y],{x,-m,m},{y,-m,m},Contours->rr,ContourShading->None]]];

fd=y^4-y^2-x^4+2x^2;
rrfd={.2,.4,1.3,3};
plotDiamonds[fd,rrfd,x,y]

gaussAssoc[f_,rr_,x_,y]:=Module[{A,r,V,d,sgn},
sgn=Switch[Sign[#],1,0,-1,1,0,0]&;
d=tDeg[f,x,y];
If[d≠ Length[rr], Echo[d,"rr size"];Abort[]];
r=Prepend[rr,0];
V=ddVerts[d];
<|Table[v->sgn[f/.Thread[{x,y}->r[[Total[Abs[v]]+1]]/(Total[Abs[v]]+.001)*v]],{v,V}]|>
];

dDiagram[gaussAssoc[fd,rrfd,x,y]]

n2s[n_]:=Module[{k,n2dig,s},
k=Which[ n<8,3,n<64,6,n<1024 ,10,n<32768,15,n<2097152,21,n<2^28,28, n<2^36,36,n<2^45,45,True,Fail];
n2dig=IntegerDigits[n,2,k];
s=Reverse[n2dig]
];

n2s[38]

s2n[s_]:=Module[{a},
a=Reverse[s];
FromDigits[a,2]
];

s2n[{0,1,1,0,0,1,0,0}]

Functions From Appendix on Numerical Linear Algebra

gaussNewton3[f_,g_,h_,p0_,x_,y_]:=Module[{Q,Q0,p1,p2,J0,n},
n=30;
p1=N[p0];
Catch[Do[Q={f,g,h}/.Thread[{x,y}->p1];
J0=D[{f,g,h},{{x,y}}]/.Thread[{x,y}->p1];
p2=p1-PseudoInverse[J0].Q;
If[Norm[p1-p2]<1.*^-7,Throw[p2]];
p1=p2,{n}]];
p2]

f=10+15 x+9 (x^2)-9 y;
g=13+9 x-21 y+9 (y^2);
h=-26-42 x+30 y+36 x y;
NumberForm[{x,y}/.FindRoot[{f,g},{x,-.3},{y,.7}],10]
NumberForm[gaussNewton3[f,g,h,{-.3,.7},x,y],10]

fAssoc[f_,x_,y_]:=Module[{A,B,d,FA},
A=Association[CoefficientRules[f,{x,y}]];
d=Max[Table[Total[p],{p,Keys[A]}]];
FA=Association[Table[If[MissingQ[A[p]],p->0,p->A[p]],{p,pExps[d]}]];
{FA,d}];

f
FA=fAssoc[f,x,y]

shiftFA[FA_,q_,n_]:=Module[{S,K,d},
K=Keys[FA];
S=Association[Table[p+q->FA[p],{p,K}]];
Association[Table[If[MissingQ[S[p]],p->0,p->S[p]],{p,pExps[n]}]]];

f
FR=shiftFA[FA[[1]],{1,2},FA[[2]]+3]
FromCoefficientRules[Normal[FR],{x,y}]
Expand[f*x y^2]

sylD1[f_,n_,x_,y_]:=Module[{FA,d},
{FA,d}=fAssoc[f,x,y];
If[d>n,Print["Must have n≥deg(f)"];Abort[]];Table[Values[shiftFA[FA,q,n]],{q,pExps[n-d]}]]

sylD1[f,3,x,y]//MatrixForm

sylvesterMatrix[f_,g_,n_,x_,y_]:=Join[sylD1[f,n,x,y],sylD1[g,n,x,y]];

g=2+3x+4y;
S3=sylvesterMatrix[f,g,2,x,y];
S3//MatrixForm

invSyl[row_,d_,x_,y_]:=FromCoefficientRules[Table[pExps[d][[i]]->row[[i]],{i,Length[row]}],{x,y}];

invSyl[S3[[1]],2,x,y]
invSyl[S3[[3]],2,x,y]

nDivide[h_,g_,x_,y_,tol_]:=Module[{d1,d2,P,S,f},
d1=tDeg[g,x,y];
d2=tDeg[h,x,y];
If [d1>d2,Print["Does Not Divide"];Return[Fail]];
P=PseudoInverse[N[sylD1[g,d2,x,y]],Tolerance->tol];S=Chop[sylD1[h,d2,x,y].P];f=invSyl[S[[1]],d2-d1,x,y];
If[Norm[Flatten[sylD1[f*g-h,d2,x,y]]]>d2*tol,Print["Does not divide at this tolerance"];Return[Fail]];
f];

h=Expand[f*g]
nDivide[h,f,x,y,1.*^-12]

nGCD[ff_,hh_,x_,y_,tol_]:=Module[{S,Sk,U,V,u,v,su,d1,d2,dv,f,h,n,m,k,W,lnull},
lnull[M_]:=Dimensions[M][[1]]-Length[SingularValueList[N[M],Tolerance->tol]];
d1=tDeg[ff,x,y];
d2=tDeg[hh,x,y];
If[d1≥d2,f=ff; h=hh;m=d1;n=d2,
f=hh;h=ff;m=d2;n=d1];
k=m-1;
Sk={{0}};
While[k<m+n && lnull[Sk]==0,
k++;                                                   (* Main loop  *)
Sk=sylvesterMatrix[f,h,k,x,y]];
{U,S,V}=SingularValueDecomposition[N[Sk],Tolerance->tol];
W=Flatten[Take[U,All,-1]];
dv=k-m;
If[dv≤0, Return[h]];
su=Chop[Take[W,Binomial[dv+2,2]],tol];
u=invSyl[su,dv,x,y];
nDivide[h,u,x,y,tol]]

ff=Expand[(2+3x+4y)(4x^2-2y+1)]
h
nGCD[ff,h,x,y, 1.*^-12]

sqFree[f_,x_,y_,tol_]:=Module[{S,Sk,U,V,v,su,dv,h,n,m,k,W},
h=Expand[{D[f,x],D[f,y]}.RandomReal[{-1,1},2]];
lnull[M_]:=Dimensions[M][[1]]-Length[SingularValueList[N[M],Tolerance->tol]];
m=tDeg[f,x,y];
n=tDeg[h,x,y];
k=m-1;
Sk={{0}};
While[k<m+n && lnull[Sk]==0,
k++;
Sk=sylvesterMatrix[f,h,k,x,y]];
{U,S,V}=SingularValueDecomposition[N[Sk],Tolerance->tol];
W=Flatten[Take[U,All,-1]];
dv=k-n;
su=Chop[Take[W,-Binomial[dv+2,2]],tol];
invSyl[su,dv,x,y]]

f3=Expand[y^3(x+y)^2]
sqFree[f3,x,y,1.*^-12]

difOp[f_,pt_,i_,j_]:=N[(1/Factorial[i]/Factorial[j]D[f,{x,i},{y,j}])/.Thread[{x,y}->pt]];

f
difOp[f,{-1/3,2/3},1,1]
difOp[x*y*f,{-1/3,2/3},3,1]

dAss[f_,pt_,0]=<|{0,0}->0|>;
dAss[f_,pt_,n_]:=dAss[f,pt,n]=Append[dAss[f,pt,n-1],Table[{i,n-i}->difOp[f,pt,i,n-i],{i,0,n}]];

f
dAss[f,{-1/3,2/3},0]
dAss[f,{-1/3,2/3},2]

macaulayD1[f_,p_,n_]:=Module[{M},
Table[Values[shiftFA[dAss[f,p,n-Total[q]],q,n]],{q,pExps[n-1]}]];

f
macaulayD1[f,{-1/3,2/3},2]//MatrixForm

macaulayMatrix[f_,g_,p_,n_]:=Join[macaulayD1[f,p,n],macaulayD1[g,p,n]];

f
g=13+9 x-21 y+9 (y^2)
macaulayMatrix[f,g,{-1/3,2/3},2]//MatrixForm

Options[intersectionMultiplicity]={HF->Off};
intersectionMultiplicity[f_,g_,p_,tol_,OptionsPattern[]]:=Module[{M,c,r,m,n,del, maxm,i,L,x,y},
maxm=20;
If[Abs[f/.Thread[{x,y}->p]]>tol,Print["p not zero of f"];Abort[]];
If[Abs[g/.Thread[{x,y}->p]]>tol,Print["p not zero of g"];Abort[]];
del=1;
m=1;
L={1};
For[i=1,i<maxm && del>0,i++,
M=macaulayMatrix[f,g,p,i];
c=Dimensions[M][[2]];
r=Length[Select[SingularValueList[N[M]],#>tol&]];
n=c-r;
del=n-m;
m=n;
AppendTo[L,m]];
If[m==maxm,Print["Multiplicity Greater or equal to ",maxm]];If[OptionValue[HF]==On,Print["depth = ",i-2,"  HF =  ",L]];
m]

f=10+15 x+9 (x^2)-9 y;
g =3+9 x-21 y+9 (y^2);
intersectionMultiplicity[f,g,{-1/3,2/3},1.*^-12]
intersectionMultiplicity[f,g,{-1/3,2/3}, 1.*^-12,HF->On]

p2aMatrix[F_,d_,t_,x_,y_]:=Module[{m,n,M,c,ca,ra},
m=Max[Exponent[F[[1]],t],Exponent[F[[2]],t]];
n=m*d;
M={};
Do[c=CoefficientRules[Expand[(x^p[[1]] y^p[[2]])/.Thread[{x,y}->F]],t];
ca=Association[c];
ra=Association[Table[If[MissingQ[ca[{k}]],{k}->0,{k}->ca[{k}]],{k,0,n}]];
AppendTo[M,Values[ra]],{p,pExps[d]}];
M];

PA=p2aMatrix[{t^2,t^3},3,t,x,y]

p2aEq[u_,x_,y_]:=Module[{P,d},
If[Length[Dimensions[u]]>1,Print["u must be a vector"];Abort[]];
d=Switch[Length[u],3,1,6,2,10,3,15,4,21,5,28,6,36,7,45,8,55,9];
If[Length[d]>1,Abort[]];
P=pExps[d];
FromCoefficientRules[Table[P[[k]]->u[[k]],{k,(d+2)(d+1)/2}],{x,y}]];

NS=NullSpace[Transpose[PA]]
p2aEq[NS[[1]],x,y]

singularFactor[f_,x_,y_,tol_]:=Module[{h,csp,d,sp,M,M1,cr,c0,rhs,rr,XX,h1},
csp=Chop[complexSingularPoints[f,x,y]];
d=tDeg[f,x,y];
If[d<3,Print["degree must be at least 3"];Abort[]];
If[Length[csp]<((d-1)(d-2)/2+1),Print["Not enough affine singular points"];Abort[]];
csp=RandomSample[csp,(d-1)(d-2)/2+1];
sp=If[d==3,csp,Join[csp,realRegularPoints[f,x,y,d-3]]];
XX=pExps[d-2];
M=Table[Which[e=={0,0},1,e[[1]]==0,p[[2]]^e[[2]],e[[2]]==0,p[[1]]^e[[1]],True,p[[1]]^e[[1]]*p[[2]]^e[[2]]],{p,sp},{e,XX}];
rr=Table[0,{Length[XX]}];
rr[[-1]]=1;
M1=Append[M,RandomReal[{-1,1},Length[XX]]];
rhs=Table[0,{Length[XX]}];
rhs[[-1]]=1;
cr=LeastSquares[M1,rhs];
c0=Min[Select[Abs[cr],#>.1&]];
h=FromCoefficientRules[Table[XX[[i]]->cr[[i]],{i,Length[XX]}],{x,y}];
h=ComplexExpand[Re[h]];
h1=nGCD[f,h,x,y,tol];
If[tDeg[h1,x,y]==0,Print["GCD not found, Returning h"];Return[{h,1}]];
Chop[{h1,nDivide[f,h1,x,y,tol]},10*tol]];

sf=Expand[(x^2+y^2-1)(x^2-y)]
singularFactor[sf,x,y,1.*^-7]

singularFactorInf[f_,x_,y_,tol_]:=Module[{d,fh,f1,z,A,h,k},
d=tDeg[f,x,y];
A=iTransform[RandomReal[{-1,1},3].{x,y,1},x,y];
f1=FLT[f,A,x,y];
{h,k}=singularFactor[f1,x,y,tol];
Chop[{FLT[h,Inverse[A],x,y],FLT[k,Inverse[A],x,y]},10*tol]
];

singularFactorInf[x^3-x y,x,y,1.*^-7]

(9A)   Appendix:  Subroutines for Chapter 9

Note :  If the code for Chapter 9 works for all examples then the code below should also be OK.  So we are not giving examples since it is unlikely that you will be using any of these directly.  But if there is a problem with the code above for Chapter 9 the problem could be here.

ddVerts[d_]:=Module[{V},
V[0]={{0,0}};
V[n_]:=V[n]=Join[V[n-1],Reverse[Table[{j,n-j},{j,0,n}]],Table[{-j,n-j},{j,1,n}],Table[{-j,j-n},{j,0,n-1}],Table[{j,j-n},{j,1,n-1}]];
V[d]];

zeroA[d_]:=Association[Table[v->0,{v,ddVerts[d]}]]

plotAssoc[A_,c1_,c2_]:=Module[{K,C1,C2},
K=Keys[A];
C1=Select[K,A[#]==0&];
C2=Select[K,A[#]==1&];
If[Length[C1]==0, Return[ListPlot[K,Axes->False, PlotStyle->c2]]];
If[Length[C2]==0, Return[ListPlot[K,Axes->False, PlotStyle->c1]]];
ListPlot[{C1,C2},PlotMarkers->{●,●},PlotStyle->{c1,c2},Axes->False]];

sqrSeq[p_,A_]:={A[p],A[p+{1,0}],A[p+{0,1}],A[p+{1,1}]};

ps0[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a+1,b+1},{a,b+1},{a,b}}]}];

ps1[a_,b_]:=Graphics[{Thick,Blue,Circle[{a,b},.5,{0,Pi/2}], Thin,Gray,Line[{{a,b},{a+1,b},{a+1,b+1},{a,b+1},{a,b}}]}];

ps2[a_,b_]:=Graphics[{Thick,Blue,Circle[{a+1,b},.5,{Pi/2,Pi}], Thin,Gray,Line[{{a,b},{a+1,b},{a+1,b+1},{a,b+1},{a,b}}]}];

ps3[a_,b_]:=Graphics[{Thick,Blue,Line[{{a,b+.5},{a+1,b+.5}}], Thin,Gray,Line[{{a,b},{a+1,b},{a+1,b+1},{a,b+1},{a,b}}]}];

ps4[a_,b_]:=Graphics[{Thick,Blue,Circle[{a,b+1},.5,{3Pi/2,2Pi}], Thin,Gray,Line[{{a,b},{a+1,b},{a+1,b+1},{a,b+1},{a,b}}]}];

ps5[a_,b_]:=Graphics[{Thick,Blue,Line[{{a+.5,b},{a+.5,b+1}}], Thin,Gray,Line[{{a,b},{a+1,b},{a+1,b+1},{a,b+1},{a,b}}]}];

ps6[a_,b_]:=Graphics[{Thick,Blue,Circle[{a+1,b+1},.5,{Pi,3Pi/2}], Thin,Gray,Line[{{a,b},{a+1,b},{a+1,b+1},{a,b+1},{a,b}}]}];

ps7[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a+1,b+1},{a,b+1},{a,b}}],Thick,Cyan,Dashed,Circle[{a+1,b+1},.5,{Pi,3Pi/2}], Circle[{a,b},.5,{0,Pi/2}],Magenta,Circle[{a+1,b},.5,{Pi/2,Pi}],Circle[{a,b+1},.5,{3Pi/2,2Pi}]}];

sqAssoc=<| {0,0,0,0}->ps0,{1,1,1,1}->ps0,{1,0,0,0}->ps1,{0,1,1,1}->ps1, {0,1,0,0}->ps2,{1,0,1,1}->ps2,{1,1,0,0}->ps3,{0,0,1,1}->ps3,{0,0,1,0}->ps4,{1,1,0,1}->ps4,{1,0,1,0}->ps5,{0,1,0,1}->ps5,{1,1,1,0}->ps6,{0,0,0,1}->ps6,{0,1,1,0}->ps7,{1,0,0,1}->ps7 |>;

triseq1[p_,A_]:={A[p],A[p+{1,0}],A[p+{0,1}]};
triseq2[p_,A_]:={A[p],A[p+{-1,0}],A[p+{0,1}]};
triseq3[p_,A_]:={A[p],A[p+{-1,0}],A[p+{0,-1}]};
triseq4[p_,A_]:={A[p],A[p+{1,0}],A[p+{0,-1}]};

pt10[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a,b+1},{a,b}}]}];

pt11[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a,b+1},{a,b}}],
Thick,Blue,Circle[{a,b},.5,{0,Pi/2}]}];

pt12[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a,b+1},{a,b}}],Thick,Blue,Line[{{a,b+.5},{a+.5,b+.5}}]}];

pt13[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a,b+1},{a,b}}],Thick,Blue,Line[{{a+.5,b},{a+.5,b+.5}}]}];

t1Assoc=<|{0,0,0}->pt10,{1,1,1}->pt10,{1,0,0}->pt11,{0,1,1}->pt11,
{1,1,0}->pt12,{0,0,1}->pt12, {0,1,0}->pt13, {1,0,1}->pt13 |>;

pt20[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a-1,b},{a,b+1},{a,b}}]}];
pt21[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a-1,b},{a,b+1},{a,b}}],
Thick,Blue,Circle[{a,b},.5,{Pi/2,Pi}]}];
pt22[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a-1,b},{a,b+1},{a,b}}],Thick,Blue,Line[{{a,b+.5},{a-.5,b+.5}}]}];
pt23[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a-1,b},{a,b+1},{a,b}}],Thick,Blue,Line[{{a-.5,b},{a-.5,b+.5}}]}];
t2Assoc=<|{0,0,0}->pt20,{1,1,1}->pt20,{1,0,0}->pt21,{0,1,1}->pt21,
{1,1,0}->pt22,{0,0,1}->pt22, {0,1,0}->pt23, {1,0,1}->pt23 |>;
pt30[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a-1,b},{a,b-1},{a,b}}]}];
pt31[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a-1,b},{a,b-1},{a,b}}],
Thick,Blue,Circle[{a,b},.5,{Pi,3Pi/2}]}];
pt32[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a-1,b},{a,b-1},{a,b}}],Thick,Blue,Line[{{a,b-.5},{a-.5,b-.5}}]}];
pt33[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a-1,b},{a,b-1},{a,b}}],Thick,Blue,Line[{{a-.5,b},{a-.5,b-.5}}]}];
t3Assoc=<|{0,0,0}->pt30,{1,1,1}->pt30,{1,0,0}->pt31,{0,1,1}->pt31,
{1,1,0}->pt32,{0,0,1}->pt32, {0,1,0}->pt33, {1,0,1}->pt33 |>;
pt40[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a,b-1},{a,b}}]}];
pt41[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a,b-1},{a,b}}],
Thick,Blue,Circle[{a,b},.5,{3Pi/2,2Pi}]}];
pt42[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a,b-1},{a,b}}],Thick,Blue,Line[{{a,b-.5},{a+.5,b-.5}}]}];
pt43[a_,b_]:=Graphics[{Thin,Gray,Line[{{a,b},{a+1,b},{a,b-1},{a,b}}],Thick,Blue,Line[{{a+.5,b},{a+.5,b-.5}}]}];
trAssoc=<|{0,0,0}->pt40,{1,1,1}->pt40,{1,0,0}->pt41,{0,1,1}->pt41,
{1,1,0}->pt42,{0,0,1}->pt42, {0,1,0}->pt43, {1,0,1}->pt43 |>;

Options[dDiagram]={pcolr->Orange,ncolr->Brown};
dDiagram[A_,OptionsPattern[]]:=Module[{d,sqrs,t1,t2,t3,t4,PS,PT,PB,PA},
d=Max[Total[#]&/@Keys[A]];
sqrs=Select[Keys[A],KeyExistsQ[A,#+{1,0}]&&KeyExistsQ[A,#+{1,1}]&&KeyExistsQ[A,#+{0,1}]&];
PS:=Show[Table[sqAssoc[sqrSeq[p,A]]@@p,{p,sqrs}]];
t1=Select[Keys[A],#[[1]]+#[[2]]==d-1&];
t2=Select[Keys[A],-#[[1]]+#[[2]]==d-1&];
t3=Select[Keys[A],-#[[1]]-#[[2]]==d-1&];
t4=Select[Keys[A],#[[1]]-#[[2]]==d-1&];
PT:=Show[Table[t1Assoc[triseq1[p,A]]@@p,{p,t1}],Table[t2Assoc[triseq2[p,A]]@@p,{p,t2}],Table[t3Assoc[triseq3[p,A]]@@p,{p,t3}],Table[trAssoc[triseq4[p,A]]@@p,{p,t4}]];
PB:=ListLinePlot[{{{-d,0},{d,0}},{{0,-d},{0,d},{d,0},{0,-d},{-d,0},{0,d}}},PlotStyle->Black,Axes->False];
PA:=plotAssoc[A,OptionValue[pcolr],OptionValue[ncolr]];
Show[PS,PT,PB,PA, PlotRange->d+1]
];

Options[sDiagram]={pcolr->Orange,ncolr->Brown};
sDiagram[A_,OptionsPattern[]]:=Module[{d,sqrs,PS,PA},
sqrs=Select[Keys[A],KeyExistsQ[A,#+{1,0}]&&KeyExistsQ[A,#+{1,1}]&&KeyExistsQ[A,#+{0,1}]&];
PS:=Show[Table[sqAssoc[sqrSeq[p,A]]@@p,{p,sqrs}]];
PA:=plotAssoc[A,OptionValue[pcolr],OptionValue[ncolr]];
Show[PS,PA,Axes->True]
];

descartesDiagram[sa_]:=Module[{d,sqrs,PS,PA,ssa,K,A},
d=Max[Flatten[Keys[sa]]];
ssa=sa;
AssociateTo[ssa,Table[{-j}->If[Mod[j,2]==0,sa[{j}],1-sa[{j}]],{j,d}]];
descartesSketch[ssa,d,pcolr->Green,ncolr->Red]
]

Options[descartesSketch]={pcolr->Orange,ncolr->Brown};
descartesSketch[A_,d_,OptionsPattern[]]:=Module[{SA,T,B,AA},
T=Association[Table[{x,1}->1,{x,-d,d}]];
B=Association[Table[{x,-1}->0,{x,-d,d}]];
AA=Association[Table[{x,0}->A[{x}],{x,-d,d}]];
SA=Union[B,AA,T];
sDiagram[SA,pcolr->OptionValue[pcolr],ncolr->OptionValue[ncolr]]
]

s2A[s_]:=Module[{n,d,A,K},
n=Length[s];
d=Which[ n==3,1,n==6,2,n==10 ,3,n==15,4,n==21,5,n==28,6, n==36,7,n==45,8,n==55,9,True,Fail];
K=Select[ddVerts[d],(#[[1]]≥0 && #[[2]]≥0)&];
A=Association[Table[K[[i]]->s[[i]],{i,n}]]];

dvAssoc[A_]:=Module[{d,K,V,v,pver},
K=Keys[A];
d=Max[Table[Total[p],{p,K}]];
v=ddVerts[d];
V=<| |>;
Do[Which[p[[1]]≥0 && p[[2]]≥0,AppendTo[V,p->A[p]],
p[[1]]<0&&p[[2]]>=0&&Mod[-p[[1]],2]==0, AppendTo[V,p->A[Abs[p]]],
p[[1]]<0&&p[[2]]>=0&&Mod[-p[[1]],2]==1, AppendTo[V,p->(1-A[Abs[p]])],
p[[2]]<0&&p[[1]]>=0&& Mod[-p[[2]],2]==0, AppendTo[V,p->A[Abs[p]]],
p[[2]]<0&&p[[1]]>=0&& Mod[-p[[2]],2]==1, AppendTo[V,p->(1-A[Abs[p]])],
Sign[(-1)^p[[1]]]*Sign[(-1)^p[[2]]]==1, AppendTo[V,p->A[Abs[p]]],
Sign[(-1)^p[[1]]]*Sign[(-1)^p[[2]]]==-1, AppendTo[V,p->(1-A[Abs[p]])]],{p,v}];
V];

A2=coefficientSigns[x^2-x y+y^2-x+y-1,{x,y}]
dvAssoc[A2]

Examples from TMJ article

Plotting the hyperbola on the Möbius band.

f=x^2-y^2-4;
{A1,g1,bp1,ap1}=inspectInfPoint[f,{1,1,0},1,x,y,output->"data"];
{A2,g2,bp2,ap2}=inspectInfPoint[f,{1,-1,0},1,x,y,output->"data"];
{RegionPlot[f>0, {x,-3,3},{y,-3,3},Epilog->{{Red, PointSize[Medium],Point[ap1]},{Purple,PointSize[Medium],Point[ap2]}},ImageSize->Small],inspectInfPoint[f,{1,1,0},1,x,y,output->"plot"],inspectInfPoint[f,{1,-1,0},1,x,y,output->"plot"]}

The code below depends on the order of points in ap1, ap2, bp1 and bp2.   Unfortunately the code above may return different orders,  so we sort to get a standard ordering.

sp1=SortBy[ap1,#[[1]]&]
sp2=SortBy[ap2, #[[1]]&]

cp1=SortBy[bp1, First[flti[#,Inverse[A1]]]&]
cp2=SortBy[bp2, First[flti[#,Inverse[A2]]]&]

Note the following gives sp1 and sp2

flti[#,Inverse[A1]]&/@cp1
flti[#,Inverse[A2]]&/@cp2

The following diagram shows the 5 segments that must be traced with their boundaries.

RegionPlot[x^2-y^2>4,{x,-5,5},{y,-5,5},Epilog->{{Text["L3p", {-4.5,3.5}]},{Text["L1p",{-3,1.5}]},{Text["L2",{2.5,0}]},{Text["L4p",{4.5,3.5}]},{Text["L1n",{-3,-1.5}]},{Text["L3n",{-4.5,-3.5}]},{Text["L4n",{4.5,-3.5}]},{Text["sp1[[1]]",{-3,-3}]},{Text["{-2,0}",{-1.4,0}]},{Text["sp2[[1]]",{-3,3}]},{Text["{-1,1,0}",{-4.5,4.8}]},{Text["{1,1,0}",{4.5,4.8}]},{Text["{1,-1,0}",{4.5,-4.8}]},{Text["{-1,-1,0}",{-4.5,-4.8}]},{Text["sp1[[2]]",{2.8,3}]},{Text["sp2[[2]]",{2.8,-3}]}}]

L1p=pathFinderT[-f,{-2,0},sp2[[1]],.05,x,y]

L1n=pathFinderT[f,{-2,0},sp1[[1]],.05,x,y]

L2=pathFinderT[f,sp2[[2]],sp1[[2]],.2,x,y,maxit->50]

L3p=flti[#,Inverse[A2]]&/@pathFinderT[g2,cp2[[1]],{0,0},.075,x,y]

L3n=flti[#,Inverse[A1]]&/@pathFinderT[-g1,cp1[[1]],{0,0},.075,x,y]

L4p=flti[#,Inverse[A1]]&/@pathFinderT[g1,cp1[[2]],{0,0},.075,x,y]

L4n=flti[#,Inverse[A2]]&/@pathFinderT[-g2,cp2[[2]],{0,0},.075,x,y]

Check Work

Show[ContourPlot[x^2-y^2==4,{x,-5,5},{y,-5,5},ContourStyle->LightGray],ListLinePlot[{L1p,L1n,L2,Drop[L3p,-1],Drop[L3n,-1],Drop[L4p,-1],Drop[L4n,-1]}]]

M1p=moebiusPhi[#]&/@L1p;
M1p[[1]]=N[{Pi,ArcTan[2]/Pi}];
M1p

M1n=moebiusPhi[#]&/@L1n;
M1n[[1]]=N[{-Pi,ArcTan[2]/Pi}];
M1n

M2=moebiusPhi[#]&/@L2

Note that M3p and M3n both have the same ambiguous point (1,1,0) but need to be replaced by different points.

M3p=moebiusPhi[#]&/@L3p;
M3p[[-1]]={N[ArcTan[-1,1]],.5};
M3p

M3n=moebiusPhi[#]&/@L3n;
M3n[[-1]]={N[ArcTan[-1,-1]],.5};
M3n

M4n=moebiusPhi[#]&/@L4n;
M4n[[-1]]={N[ArcTan[1,-1]],.5};
M4n

M4p=moebiusPhi[#]&/@L4p;
M4p[[-1]]={N[ArcTan[1,1]],.5};
M4p

Check work

Show[moebiusRMPP, ListLinePlot[{M1p,M1n,M2,M3p,M3n,M4p,M4n},PlotStyle->Directive[Thick,Blue]]]

MM1p=moebiusB@@#&/@M1p;
MM1n=moebiusB@@#&/@M1n;
MM2=moebiusB@@#&/@M2;
MM3n=moebiusB@@#&/@M3n;
MM3p=moebiusB@@#&/@M3p;
MM4n=moebiusB@@#&/@M4n;
MM4p=moebiusB@@#&/@M4p;

Show[moebiusPLMB, Graphics3D[{{Thick,Black,Line[MM1p]},{Thick,Black,Line[MM1n]},{Thick,Black,Line[MM2]},{Thick,Black,Line[MM3p]},{Thick,Black,Line[MM3n]},{Thick,Black,Line[MM4p]},{Thick,Black,Line[MM4n]}}],ViewPoint->{0.4,3.13,-1.23},ImageSize->Medium]