Mathematica 12, 12.1       5/21/21

These functions were coded by me in Mathematica 11.3-12.0 as simple illustrations of algorithms in A Numerical Approach to Real Algebraic Curves with the Wolfram Language only.

This appendix, also periodically called the global functions notebook, contains examples of most of the functions I have written for this book. 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.

Chapter 1

(1.1) line: Gives a line through two points, finite or infinite. You can enter points as either pairs or triples. Note that there is randomization; rerunning the same data in a session will give the same answer, but it will vary by a constant multiple in different sessions.

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

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

(1.3) newtonHyperbolas: 1≤n≤7 degree 1, 8≤n≤63 degree 2, 64≤n≤1023 degree 3, 1024≤n≤32767 degree 4, 32768≤n≤2097151 degree 5, 2097152≤n<(2^28) degree 6.

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]

(1.4) gaussCurve: p is a polynomial in z with possibly complex coefficients. ComplexExpand assumes that variables x, y are real.

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

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

(1.5) sqFree: sqFree[g,x,y,tol] gives the square free decomposition of polynomial g(x,y) to tolerance tol. See Appendix (A.3) below for code and example.

Chapter 2

(2.2) newton1: One-dimensional Newton method.

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]

(2.2) gaussNewton2: Two-dimensional Newton method for two functions. In practice, use the built-in FindRoot instead.

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]

(2.4) tDeg: The total degree of a bivariate polynomial.

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

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

(2.4) pExps: A list of exponent pairs in a bivariate polynomial of degree d can be constructed recursively.

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]

pMons: Monomials of given degree and variables

pMons[d_,x_,y_]:=Table[x^p[[1]] y^p[[2]],{p,pExps[d]}]

(2.4) randomRealPolynomial: A procedure to construct random numerical bivariate polynomials of degree d.

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

randomRealPolynomial[3,x,y]

(2.4) randomIntegerPolynomial: The function to produce a random integer polynomial has a new variable kk, which is the maximum number of nonzero terms which must be a positive integer. At least one term will have total degree d.

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]

(2.4) aCurve: Gives a curve of degree d through (({{d+2}, {2}}))-1 points provided that these points are sufficiently random.

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]]

Chapter 3

(3.1) tLine: Tangent line at a point.

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

f= 2 -2 x+ (x^2)- y;
g=tLine[f, {3,5},x,y];
Expand[g/-Coefficient[g,y]]

(3.1) tangentRealPoints: Points of f=0 with tangents parallel to line k.

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= 2 -2 x+ (x^2)- y;
tangentRealPoints[f, 4x-y,x,y]

(3.2) realSingularPoints: Real affine singular points. Note that we have a similar function for complex projective singular points in (6.4).

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]];

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

(3.3) form: Homogeneous form consisting of all monomials of degree k in polynomial h.

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]

(3.3) maxForm: Maximal form of polynomial h.

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

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

(3.3) infiniteRealPoints: Note that we have a similar function for complex projective singular points in Chapter 6.

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=(x-1)^3-(y-2)^2;
infiniteRealPoints[f,x,y]

(3.4) pointMinForm: Minimum form at a point not necessarily 0. The tolerance tol can be changed.

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=(x-1)^3-(y-2)^2;
pointMinForm[f,{1,2},x,y]

(3.5) criticalPoints:

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-1)^3-(y-2)^2;
cp=criticalPoints[f,x,y]

(3.5) closetPoint: Closest point on a curve to q. If several are equidistant, only one is given.

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=(x-1)^3-(y-2)^2;
closestPoint[f,{5,3},x,y]

Chapter 4

(4.2) gaussPlot: Utility function to plot a curve inside a disk, showing critical and/or boundary points, if given. f is the curve, r is the radius of the disk, and pts is the points, possibly an empty set.

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]

(4.3) gTvec: Directed tangent vector clockwise about components of f>0.

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

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

(4.3) gPoints: gPoints of the curve f inside a disk of radius r.

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]];

g=x^3-x+1-y^2;
gPoints[g,4,x,y]

(4.4) pathFinder: pathFinder uses naive tracing using the closest point. f is the curve, p is the initial point, q is the final point, and s is the step size.

pathFinder[f_,p_, q_,s_,x_,y_]:= Module[{L,p0,p1,u, maxit,k},
maxit=30;
p0=p;
L=Reap[Sow[p];
k=0;
While[Norm[q-p0]>2s && k<maxit,
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≥ maxit, Print["Warning: iteration limit reached at",p1]];
L];

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

(4.4) pathFinderT: pathFinderT traces using the tangent vector given by gTvec. It has the same arguments as pathFinder. It will be faster, especially in higher degrees, but may not work as well in some cases. An option can be used to change the maximum number of iterations.

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["Warning: iteration limit reached at",p1]];
L];

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

(4.5) plotPath: Curve plotting using several points. f is the curve, G is the graph, V is the actual points of the curve used as vertices, and s is the 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=x^3-x+1-y^2;
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=x^3-x+1-y^2;
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) deTrace: Curve tracing using NDSolveValue. f is the curve, p is the initial point, s is the maximal distance from the initial point, and r is the maximum distance from the origin. deTrace returns a pair of interpolating functions, i.e. x, y parameters of the curve. Interpolating function approximates the curve but does not necessarily give actual points on the 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 that may or may not appear here depending on the random value used.

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

(4.6) deEnd: deEnd uses end point and maximum value for the domain of the interpolating function pair from deTrace.

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

(* using  output {X,Y} from example above *)
deEnd[{X,Y}]

(4.6) pathFinderDE: Here, r is the maximum distance from the origin and eps is the tolerance. When the distance from a calculated point to q is determined, the algorithm stops.

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]]];

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

(4.6) plotPathDE: Has the same syntax as plotPath.

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]]]]]]

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

(4.6) endFinder: Routine to attempt to find the next point of a set of points CP on a curve reached by tracing from the n(^th) point of CP. The goal is met when the trace is within eps of a point in CP.

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]];

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

(4.6) nearPoints: Are points nearby a singularity pt of f but with gTvec at those points pointing away from the singularity. r is the radius of a bounding circle and ep is the distance of the nearPoints from the singularity. These are necessary because one can trace into a singularity but cannot trace out of a singularity.

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=(((x-1))^3)-(((y-2))^2);
nearPoints[f,4,x,y,{1,2},.4]

(4.6) gaussGraphNS: An attempt at a black box function that traces a curve f inside a disk of gPoints CP and returns a plot of the curve with gPoints labeled and a graph of the curve inside the disk. eps is used in nearPoints and endFinder. If run without, it will return the gPoints and graph rules. gaussGraphNS is simpler and works in the nonsingular case. gaussGraph works best for Gauss curves, possibly even with singularities.

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=x^3-x+1-y^2;
V={{0,-1},{-1,-1},{-1,1},{0,1}};
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=(((x-1))^3)-(((y-2))^2);
Vf={{2,1},{1,2},{2,3}};
gaussGraph[f,Vf,x,y,.4]

(4.6) plotPathS and plotPathDES: Two routines to plot singular curves from graphs that include user near points. The first is without the DE path tracing and the second is with.

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];

f=(((x-1))^3)-(((y-2))^2);
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]]

Chapter 5

(5.1) pLine: Projective line through two points (triples) of the real or complex projective plane.

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

pLine[{1,2,1},{0,1,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^2-x-y^3+2,x,y,z]

(5.3) dTol: Default tolerance (constant).

dTol:=1.*^-12;

(5.3) intersectionMultiplicity[f,g,p]: Gives the intersection multiplicity of curves f,g at intersection point p. See (A.4).

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

(5.3) complexSingularPoints: Finds complex affine singular points.

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]

(5.3) realRegularPoints: Finds n random real regular points for curve f of degree >1.

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]

Chapter 6

(6.2) Fractional linear transformations: FLT

flt - Transforms affine points—always returns a pair.

flti - Transforms possibly infinite points—may return a triple.

flts - Transforms symbolic points—returns a pair of possibly rational symbolic expressions.

flth - Transforms projective point to projective point.

FLT - Transforms curve equations—gives an affine equation.

FLTh - Transforms homogeneous curve equation results in a homogeneous equation.

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]
(* answer {} expected *)

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]

(6.3) pqTranslation: Translation taking point p to point q. An example of a Euclidean transformation as FLT.

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}]]

(6.3) klRotation: Rotation taking line k to line l. An example of a Euclidean transformation as FLT.

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]]

(6.3) kRelfection: Reflection in line k. An example of a Euclidean transformation as FLT.

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]]

(6.3) kShear: A standard shear moving line k to the x-axis. An example of affine transformations as FLT.

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]

(6.3) homothety: Homothety where r,s are positive reals. If one of r,s is negative, this adds in a reflection about the given axis. Both negative adds in a half-turn. An example of affine transformations as FLT.

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

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

(6.3) ip2z: Take infinite point ip to the origin. x, y are dummy variables. An example of projective transformation.

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]]

(6.3) iTransform: Takes line l to the infinite line, moving the original infinite line to the x-axis. There is a fixed random vector r1 that could be replaced by another at the user’s discretion.

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]]

(6.4) complexProjectiveSingularPoints: Here we often use a very loose tolerance. It is recommended that this function be run several times with different tolerances, but even running it twice with the same tolerance will give different results. A weak tolerance, e.g. .001, will identify all possible singular points. A tight tolerance will eliminate some that are not singular but may miss some that are.

complexProjectiveSingularPoints[h_,u_,v_,eps_]:=Module[{f,w,A,a,b,sol1,sol2,solz,sing,sing1,sing2},A={RandomReal[{-1,1},3],{0,0,RandomReal[{1,2}]},RandomReal[{-1,1},3]};
f=FLT[h,A,u,v];
{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];
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=9-10 (x^2)+(x^4)+6 (y^2)+2 (x^2) (y^2)+(y^4);
complexProjectiveSingularPoints[h,x,y,1.*^-7]

(6.4) complexProjectiveIntersectionpoints:

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]

Chapter 7

(7.1) asymptotes: See infiniteRealPoints in (3.3) of this appendix for infinite points.

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}]
];

Example: rational function y=(((x^3)+2(x^2))/((x^2)-1)).

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]

Example:  Parabola

asymptotes[y-x^2,{0,1,0},x,y]

Indicates that the asymptote is the infinite line.

(7.2) conicWithFoci: Here we find the equation of a conic given foci and a “good” point q on the conic. The two points (Subscript[p, 1]), (Subscript[p, 2]) define the major axis. By “good”, q must not be on the axis, but the closest point on the axis to q must lie between (Subscript[p, 1]), (Subscript[p, 2]) and no farther from the axis than Norm[(Subscript[p, 1])-(Subscript[p, 2])]/2. This function has options, the default is “ellipse” and the other is “hyperbola”; either may be abbreviated to three or more letters.

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"]

(7.3) cTransform: Places p at infinite point {0,1,0} with the infinite line as the tangent line to the transform of f. This is deterministic; no randomization is used.

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]

(7.3) p2cTransform: The FLT P2C takes the parabola y=(x^2) to the circle (x^2)+(y^2)=1. In particular, the infinite point {0,1,0} goes to {0,1} and the vertex {0,0} goes to {0,-1}. Note that this is a global constant rather than a function.

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

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

(7.3) flts: See (6.2) in this appendix.

(7.3) Rational Normal Curve 2D case.  In my note Interpolating Rational Normal Curves (May 2021) posted on my Space Curve page I point out that the two dimensional case follows from the material in this book.  Loosely a rational normal curve is a rational parameterized curve with numerators and denominators of degree d.  The general Theorem says that d+3 points in general position determine a unique rational normal curve.  The curve is unique, not the parameterization.  In our case d=2 five points determine such a curve.  The following function gives the algebraic and parameter versions of such a curve given the 5 points.  For more details see the link on my Space Curve page.

nrc2[S_]:=Module[{l,SS3,tst,i,acurve,ct,T,g1,g2,a,p,par,f},
SS3=Subsets[S,{3}];
tst=False;
i=1;
While[i≤10&&tst==False,tst=Abs[(line[SS3[[i,1]],SS3[[i,2]],x,y]/.Thread[{x,y}->SS3[[i,3]]])]<.001;
i++];
If[i<11,Echo["Not General Position"];Abort[]];
acurve=aCurve[S,x,y];
ct=cTransform[acurve,S[[1]],x,y];
g1=FLT[acurve,ct,x,y];
g1=Expand[g1/Coefficient[-g1,y]];
p=tangentRealPoints[g1,y,x,y][[1]];
T={{1,0,-p[[1]]},{0,1,-p[[2]]},{0,0,1}};
g2=FLT[g1,T,x,y];
a=Coefficient[g2,x^2];
par={t,a t^2};
f=flts[{t,a t^2},Inverse[T.ct]];
{Chop[acurve],f}]

S={{-1,1},{2,0},{1,1},{-1,-1},{1,-1}};
{g,f}=nrc2[S]

Show[ContourPlot[g==0,{x,-3,3},{y,-3,3},ContourStyle->Orange],ParametricPlot[f,{t,-20,20},PlotStyle->Directive[Dashed,Blue]],Graphics[{Black,PointSize[Medium],Point[S]}]]

(7.4) affineInflectionPoints: Affine real inflection points of curve f=0.

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];
ips=If[Length[sol]==0,Return[{}],{x,y}/.sol];
Chop[DeleteDuplicates[ips,Norm[#1-#2]<1.*^-10&]]
];

f=((((x^2)+(y^2)))^2)-4(x^3)(y^2);
aip=affineInflectionPoints[f,x,y]

(7.4) allInflectionPoints: Finds all inflection points including infinite real inflection points.

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=((((x^2)+(y^2)))^2)-4(x^3)(y^2);
allInflectionPoints[f,x,y]

(7.5) weirstrassNormalForm: Weierstrass normal form of a cubic f given a particular inflection point p. Also returns FLT normalizing f and the j-invariant.

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}];

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

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

Chapter 8

(8.3) ghsMap: Map from projective plane to Gauss hemisphere. p is the point and r is the radius of the inner circle. The entire Gauss hemisphere has radius 3r.

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]

(8.3) inverseGHSmap: Inverse map from Gauss hemisphere to projective plane. Note that infinite points will appear as affine points with a very large Norm. The reason is that we want this function to work with symbolic values of q.

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

inverseGHSmap[{x,y},5]

(8.5) graphSD: Symmetric difference of graphs.

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]

(8.7) moebiusB: Point on a Möbius band.

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]];

(8.7) moebiusPLMB: Draws a Möbius band with an infinite line. Use with Show.

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]

(8.7) moebiusRMPP: Outline of a rectangular hemisphere.

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/5];

moebiusRMPP

(8.7) moebiusPhi: The function for a projective plane minus the negative x-axis to a rectangle hemisphere model of a projective plane. It is assumed projective points will have third coordinate 0 or 1. If any points are on negative x-axis, including {1,0,0}, then the function will return an error message and {}. As usual, to apply to a list L of points, use M=moebiusPhi/@L. You can find the position of error points in M by Position[M,{}]. Then fix these points manually. If they occur as other than the first or last points of M, you probably need to split the list before applying MB.

moebiusPhi[p_]:=If[p[[2]]≠0,
If[Length[p]==2,{ArcTan@@N[p],N[ArcTan[Norm[p]]/Pi]},{ArcTan@@N[Take[p,2]],.5}],
Which[p[[1]]>0&&Length[p]==2, {N[Pi/2],N[ArcTan[Norm[p]]/Pi]},
p[[1]]>0, {N[Pi/2],.5},p[[1]] ≤0||Length[p]>2,Echo[{},p]]];

L=moebiusPhi/@{{1,1},{2,1},{0,1,0},{1,1}}
Show[moebiusRMPP,ListLinePlot[L]]

Chapter 9

(9.2) coefficientSigns: Creates an association directly from a polynomial in one or more variables. The polynomial is f and X is a list of the variables.

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}]

(9.2) descartesDiagram: A is an association with rules {j}->e where j goes from 0 to the degree and each e is 0 or 1 depending on whether the j(^th) coefficient is positive, 0, negative, or 1. See (9.A5) for the code.

A=\[LeftAssociation]{3}->0,{2}->0,{1}->1,{0}->0\[RightAssociation];
descartesDiagram[A]

(9.2) descartesSketch: Uses polynomial values rather than signs. d is the degree, and now the association keys include {j} for j from -d to d. See 9.A7 for the code.

B=<|{-3}->1,{-2}->0,{-1}->0, {3}->0,{2}->0,{1}->1,{0}->0|>;
descartesSketch[B,3]

(9.2) dvAssoc: The function dvAssoc creates the association used in the viroDiagram function. For some purposes, you may want to start with this association and modify it. The input is an association with all non-negative keys, such as coefficientSigns for a full polynomial above and the output is an association with all keys. This is usually combined with coefficientSigns as viroAssociation. See (9.A10) for the code.

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

dvAssoc[A2]

(9.2) dDiagram: This is the general engine behind all of the two-dimensional diamond diagrams. This has options giving the colors of the dots. In addition, there is a large number of utility functions to enhance these diagrams.  Important note:  the input is an association which must contain  keys {a,b} for all   |a|+|b|≤d  for some d>0. The values assigned to these keys must be either 0 or 1.

A=\[LeftAssociation]{0,0}->0,{1,0}->1,{0,1}->0,{-1,0}->0,{0,-1}->1,{2,0}->0,{1,1}->1,{0,2}->0,{-1,1}->0,{-2,0}->0,{0,-2}->0,{-1,-1}->1,{1,-1}->0\[RightAssociation];
DDf=dDiagram[A, pcolr->Orange, ncolr->Brown]

(9.2) viroAssociation and viroDiagram: viroDiagram draws the Viro diagram of a bivariate polynomial. Important note: the polynomial f must be full of degree d, that means that each two variable monomial must appear with a non-zero coefficient.  For sparse equations use gaussDiagram below.

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]

(9.2) pickCyan: Graphics primitives to change ambiguous cells to fixed cells. Here, {a,b} is the lower-left corner of the cell.

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}]}];

f=1-x+y+x^2-x y+y^2;
Show[viroDiagram[f,x,y],pickCyan[0,-1]]

(9.2) pickMagenta:

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}]}];

f=1-x+y+x^2-x y+y^2;
Show[viroDiagram[f,x,y],pickMagenta[0,-1]]

(9.2) pickSingularity:

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}}]}]

f=1-x+y+x^2-x y+y^2;
Show[viroDiagram[f,x,y],pickSingularity[0,-1]]

(9.2) labelI, labelII, labelIII, labelIV,label13, label24:  These do not work in Mathematica 12.  They should not have been here originally anyway since labelling infinite points in a dDiagram can be simply and more accurately done with the Graphics function Text.  Note that cordinates in dDiagrams for the dots are integer points starting at 1.  See the following example, trial and error will probably be needed to get the coordinates right

f=1-x+y+x^2-x y-y^2;
Show[viroDiagram[f,x,y], Graphics[{Black,Text["A",{1.6,.7}], Text["B",{.5,-1.8}], Text["A",{-1.75,-.5}], Text["B",{-.75,1.5}]  }]]

For backwards compatibility only here are revised versions of label????, they should not be used in new work.

labelI[a_,b_,txt_]:=Graphics[{Black,Text[txt,{a,b}]}]
labelII[a_,b_,txt_]:=Graphics[{Black,Text[txt,{a,b}]}]
labelIII[a_,b_,txt_]:=Graphics[{Black,Text[txt,{a,b}]}]
labelIV[a_,b_,txt_]:=Graphics[{Black,Text[txt,{a,b}]}]
label13[a_,b_,txt_]:=Graphics[{Black,Text[txt,{a,b}],Text[txt,{-a,-b}]}]
label24[a_,b_,txt_]:=Graphics[{Black,Text[txt,{a,b}],Text[txt,{-a,-b}]}]

(9.2) toggle: Changes values of sign associations between 0, 1.

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

Af=\[LeftAssociation]{0,0}->0,{1,0}->1,{0,1}->0,{-1,0}->0,{0,-1}->1,{2,0}->0,{1,1}->1,{0,2}->0,{-1,1}->0,{-2,0}->0,{0,-2}->0,{-1,-1}->1,{1,-1}->0\[RightAssociation];
Af1=toggle[Af,{{0,0},{1,0}}]

{dDiagram[Af],dDiagram[Af1]}

(9.3) axisPoints: Axis points are points on either the x- or y-axis. Points on the x-axis are listed first.

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]

(9.3) inspectInfPoint: inspectInfPoint puts together several other functions to allow us to quickly see what is happening at an infinite point and to give affine points that directly lead to the infinite points. ip is the infinite point and r is a positive radius usually less than one. We want a circle that intersects the curve only near the infinite points. Data shown can be changed by options: “plot” just gives plot, “data” just gives data, and the default All is both. Data consists of the transformation, the equation at infinity, the boundary points on the circle around the infinite point, and the affine coordinates of these boundary points.

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, output->"data"]

(9.5) plotDiamonds: That is, equations |x|+ |y|=r, on a region plot given a list rr of d increasing positive values of r.

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]]];

(9.5) gaussAssoc, gaussDiagram: Creates an association by sampling points on diamonds of radii in the list rr of d= degree f increasing positive numbers. This is useful for sparse exact nonsingular equations symmetric about the y-axis. Note that this association considers the function value of 0 to be positive. The resulting diamond diagram is called a gaussDiagram.

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}]|>
];

gaussDiagram[f_,rr_,x,y]:=dDiagram[gaussAssoc[f,rr,x,y]];

fd=y^2+4x-x^3;
rrfd={.8,2.3,3.8};
gaussAssoc[fd,rrfd,x,y]
gaussDiagram[fd,rrfd,x,y]

(9.6) n2s: We recall the utility functions changing between sign sequences and index numbers for Newton hyperbolas.

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}]

Subroutines for Chapter 9

Here are subroutines and code for the Chapter 9 functions listed in this appendix. If the code for Chapter 9 works for all of the 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. However, if there is a problem with the code for Chapter 9, the problem could be here.

(9.A1) ddVerts: Vertices of a diamond diagram of dimension d.

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]];

(9.A2) zeroA: Zero diagram association.

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

Raw plot of diamond diagram association. Variables c1, c2 are the positive and negative colors, e.g. orange, RGBColor[1,.5,0], etc.

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]];

(9.A3) sqrSeq: Values bounding a cell.

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

sqAssoc: To associate these primitives to sqrSeq we have the following.

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 |>;

(9.A4) ps0–ps7: Graphics primitives for cells.

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}]}];

On the boundary of the diamond diagram, we have a triangle. The corresponding primitives for triangles are as follows.

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 |>;

(9.A5) dDiagram: The main engine for diamond diagrams is as follows.

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]];
If[Length[A]≠Length[ddVerts[d]], Echo["All keys {a,b}, \[VerticalSeparator]a\[VerticalSeparator]+\[VerticalSeparator]b|≤d must be included"]; Abort[]];
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]
];

(9.A6) sDiagram: Square diagram, used for Descartes’s law, and the descartesSketch function.

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]
];

(9.A7) descartesSketch:

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]]
]

(9.A8) descartesDiagram:

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]
]

(9.A9) s2A: Conversion from Sign sequence to association.

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}]]];

(9.A10) dvAssoc: Code for dvAssoc described in Chapter 9.2, A must have key {i,j} for all non-negative {i,j}.

dvAssoc[A_]:=Module[{d,K,V,v,pver},
K=Keys[A];
d=Max[Table[Total[p],{p,K}]];
If[Length[A]≠(d+2)(d+1)/2,Echo["A must have Key {i,j} for all 0≤i,j with i+j≤d"]; Abort[]];
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];

Appendix 1

Functions from Appendix 1: Topics in Numerical Linear Algebra

(A.1) gaussNewton: Method using PseudoInverse to find the least-square solution of three nonlinear equations in two variables. To use this with two equations, use h=0 as the third equation.

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;
intersectionMultiplicity[f,g,{-1/3,2/3}, 1.*^-12]
(* exact solution *)

NumberForm[{x,y}/.FindRoot[{f,g},{x,-.3},{y,.7}],10]

NumberForm[gaussNewton3[f,g,h,{-.3,.7},x,y],10]

(A.3) fAssoc: Puts a bivariate polynomial f into association format using CoefficientRules but adding keys for terms with coefficient zero. It returns the association and degree.

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=10+15 x+9 (x^2)-9 y;
FA=fAssoc[f,x,y]

(A.3) shiftFA: Multiplication of the association form of the polynomial by a monomial given as a pair of exponents.

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=10+15 x+9 (x^2)-9 y;
FR=shiftFA[FA[[1]],{1,2},FA[[2]]+3]
FromCoefficientRules[Normal[FR],{x,y}]
Expand[f*x y^2]

(A.3) sylD1: Sylvester matrix of one bivariate polynomial f with n as the order.

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]}]]

f=10+15 x+9 (x^2)-9 y;
sylD1[f,3,x,y]//MatrixForm

(A.3) sylvesterMatrix: Sylvester matrix of f, g bivariate polynomials with n as the order.

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

f=10+15 x+9 (x^2)-9 y;
g=2+3x+4y;
S3=sylvesterMatrix[f,g,2,x,y];
S3//MatrixForm

(A.3) invSyl: Inverse Sylvester matrix converts each row of Sylvester matrix S of order d to a polynomial in x and y. It is a subroutine for dDivide and sqFree.

(* using S3 directly above *)
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]

(A.3) nDivide: Checks if g is a factor of h and returns quotient if true to tolerance tol.

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];

f=10+15 x+9 (x^2)-9 y;
g=2+3x+4y;
h=Expand[f*g]

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

(A.3) nGCD: Numerical greatest common divisor of ff, hh using the Kaltofen–Dayton algorithm (step 2 of the Zeng–Dayton algorithm). tol is the tolerance.

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=20+60 x+63 (x^2)+27 (x^3)+22 y+33 x y+36 (x^2) y-36 (y^2);
nGCD[ff,h,x,y, 1.*^-12]

Expand[%/%[[1]]*2]

(A.3) sqFree: Square-free decomposition of numerical polynomial f to tolerance tol.

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]

(A.4) difOp: Differential operator where f is the bivariate function, pt is the point, i is the order of x, and j is the order of y. Variables x, y are assumed.

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

f=10+15 x+9 (x^2)-9 y;
difOp[f,{-1/3,2/3},1,1]
difOp[x*y*f,{-1/3,2/3},3,1]

(A.4) dAss: A recursive subroutine.

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=10+15 x+9 (x^2)-9 y;
dAss[f,{-1/3,2/3},0]
dAss[f,{-1/3,2/3},2]

(A.4) macaulayD1: A subroutine for a matrix of a single polynomial.

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

f=10+15 x+9 (x^2)-9 y;
macaulayD1[f,{-1/3,2/3},2]//MatrixForm

(A.4) macaulayMatrix: Macaulay matrix of pair of functions f, g at point p with order n.

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

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

(A.4) intersectionMultiplicity: Intersection multiplicity of intersection point p of curves f, g. A default tolerance of (10^(-9)) is suggested. This version does not, by default, print out intermediate values of null space, the depth, or the breadth, but the option HF->On will give this information.

Options[intersectionMultiplicity]={HF->Off};
intersectionMultiplicity[f_,g_,p_,tol_,OptionsPattern[]]:=Module[{M,c,r,m,n,del, maxm,i,L},
maxm=tDeg[f,x,y]+tDeg[g,x,y];
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 =13+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]

(A. 5) p2aMatrix: A subroutine for converting polynomial parametric functions to algebraic curves.

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]

(A.5) p2aEq: A subroutine that converts a sequence of coefficients to an equation. The typical use is to apply it to a row of the left null space of p2aMatrix.

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}]];

(* uses matrix PA directly above *)
NS=NullSpace[Transpose[PA]]

p2aEq[NS[[1]],x,y]

(A.6) singularFactor: Factors curve f when it contains too many singularities. Its use requires functions from Chapter 5, and it may be dependent on the tolerance tol.

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]

(A.6) singularFactorInf: Treats the case where some singularities are infinite.

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]