Asymptotes, Conics, and Cubics

In this chapter, we apply our fractional linear transformations to some classical algebraic geometry. A major goal is to parameterize these low-degree curves. In the last section, we give Arthur Caley’s proof that a curve of genus 0 is rational by giving the parameterization using simple cases.

7.1 Asymptotes and Rational Function Curves

We have touched on the subject of asymptotes briefly in the past, but in this section, we will carefully define them and give a procedure to calculate them.

Asymptotes are an affine concept but are associated with infinite points. If the infinite point is regular, there is a single asymptote, which is the tangent line to the curve at the infinite point. However, if the infinite point is singular, then there may be several asymptotes.

The standard example of the latter case is a rational function curve with several vertical asymptotes.

f=y(x^2-4)x-1;

ContourPlot[{f==0,x==0,x==2,x==-2,y==0},{x,-3,3},{y,-6,6},ContourStyle->{Black,Directive[Orange,Dashed],Directive[Orange,Dashed],Directive[Orange,Dashed],Directive[Green,Dashed]}]

By elementary algebra considerations, we have three vertical asymptotes associated with the infinite point {0,1,0} and one horizontal asymptote associated with the infinite point {1,0,0}.

A=ip2z[{1,0,0},x,y];
B=ip2z[{0,1,0},x,y];
f1=FLT[f,A,x,y]
f2=FLT[f,B,x,y]

{ContourPlot[{f1==0,y==0},{x,-1,1},{y,-1,1},ContourStyle->{Blue,Dashed}],ContourPlot[{f2==0,y==0},{x,-1,1},{y,-1,1},ContourStyle->{Blue,Dashed}]}

In these plots, the dashed line represents the affine infinite line and the point {0,0} represents the original infinite point. In the left plot, we see that there is a regular point with what appears to be the vertical tangent line x=0.

FLT[x,Inverse[A],x,y]

FLT[x,Inverse[ip2z[{1,0,0},x,y]],x,y]

It is otherwise known by the affine line y=0. The right plot shows an ordinary singularity at {0,0}, which represents the original infinite point {0,1,0}.

f2min = -2x^3+8x y^2;
Factor[f2min]

FLT[x,Inverse[B],x,y]

FLT[x-2y,Inverse[B],x,y]

FLT[x+2y, Inverse[B],x,y]

These are otherwise given by x=0, x=2, x=-2, respectively.

form[f_,k_,x_,y_]:=FromCoefficientRules[Select[CoefficientRules[f,{x,y}],Total[#[[1]]]==k&],{x,y}];
MaxForm[f_,x_,y_]:=Module[{A,Bd},
A=<|CoefficientRules[f,{x,y}]|> ;
d=Max[Table[Total[k],{k,Keys[A]}]];
Form[f,d,x,y]];
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]];

Now we can find asymptotes. This simple procedure assumes the point p is an infinite point of the curve f. The user may wish to add checks on these assumptions. Recall that the graph of a form will be a (possibly empty) union of lines through the origin. We intersect this union with a random line that does not go through the origin. The equation of each line is a factor of the form.

In Chapter 3, we defined a function finding infinite real points using a random real line. If working with a live notebook, the following code should be initialized now.

rnl7=-1+0.6916206531506335` x-0.6706463898598312` y

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

asymptotes[f_,p_,x_,y_]:=Module[{g,gm,A,sol,k},
A=ip2z[p,x,y];
g=FLT[f,A,x,y];
gm=pointMinForm[g,{0,0},x,y];
sol=DeleteDuplicates[{x,y}/.NSolve[{gm,rnl7},{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}]
];

We first look at some rational function curves.

f=y x (x+1)(x-2) -(2x^4-1);
inf=infiniteRealPoints[f,x,y]

asy1=asymptotes[f,inf[[1]],x,y]
asy2=asymptotes[f,inf[[2]],x,y]

a1=Expand[asy1[[1]]/Coefficient[asy1[[1]],y]]
a2=Expand[asy2[[1]]/Coefficient[asy2[[1]],x]]
a3=Expand[asy2[[2]]/Coefficient[asy2[[2]],x]]
a4=Expand[asy2[[3]]/Coefficient[asy2[[3]],x]]

ContourPlot[{f==0,a1==0,a2*a3*a4==0},{x,-5,8},{y,-10,20},ContourStyle->{Black,Dashed,Dashed}]

g= y((x^2)+1)-(x^4);
inf=infiniteRealPoints[g,x,y]

asy=asymptotes[g,inf[[1]],x,y]

In Chapter 1, we mentioned the line c=0 was the infinite line. So we have the one infinite point {0,1,0} with infinite asymptote.

Before leaving rational function curves, we should note that from the rational function point of view, the function is continuous and differentiable for all points (except when the denominator is zero, in which case there are no affine crossing points), but making the rational function into an algebraic curve by multiplying by the denominator can, in some cases, introduce a vertical component that meets other components at an affine point. For example, a curve from the family of curves known as the Trident of Newton can have an affine singularity of this type.

f= 3x-3x^2+2(x^3)-x y
{Plot[(3x-3x^2+2x^3)/x,{x,-4,4},PlotRange->{0,8},Axes->False,Frame->True, AspectRatio->1],ContourPlot[f==0,{x,-4,4}, {y,0,8}]}

There is a singularity at {0,3} as well as an infinite singular point at {0,1,0}.

h=0.9` -2 (x^2)+(x^4)-2 (y^2)-2 (x^2) (y^2)+(y^4);

inf=infiniteRealPoints[h,x,y]

asy1=asymptotes[h,inf[[1]],x,y]
asy2=asymptotes[h,inf[[2]],x,y]

{ContourPlot[{h==0,asy1[[1]]==0,asy1[[2]]==0}, {x,-3,3},{y,-3,3},ContourStyle->{Thick,Dashed,Dashed}],ContourPlot[{h==0,asy2[[1]]==0,asy2[[2]]==0},{x,-3,3},{y,-3,3},ContourStyle->{Thick,Dashed,Dashed}]}

Finally, we remind the reader to not confuse rational function curves with rational curves. A rational curve is a curve that has a rational parameterization. Of course, a rational function curve y g[x]-f [x] (where f, g are polynomials) has parameterization ϕ[t]={t,((f[t])/(g[t]))}, so it is a rational curve, but the converse is not true; we will see in the next section that all conics have a rational parameterization. In the last section of this chapter, we will discuss a theorem of Caley, which shows that irreducible curves with a maximum number of singularities (all double points) are rational curves.

7.2 Euclidean Geometry of Conics

In elementary algebra, real affine nonsingular conics appear in the form

((x^2)/(a^2))+((y^2)/(b^2))=1 for an ellipse,

((x^2)/(a^2))-((y^2)/(b^2))=1 for a hyperbola or

y=a(x^2) for a parabola.

In this section, we see that every nonsingular conic (exact or numeric) is equivalent to a conic in one of these forms. In an advanced course for mathematics teachers in the early ’90s, I attempted an exact example by hand, which ended up taking an entire class period. Or see the recent book by Thomas Garrity (2013), where in Section 1.2, no actual example was even given—neither in the text nor as an exercise. In precalculus or calculus, it is often shown how conics in the previous forms can be described in terms of foci and directrix. I have not seen a description of how to do this for a general conic, but since we’re using FLT, we can get the actual transformation as well as the final form, and by taking the inverse, we can find the foci and directrix.

We start with a conic f=(Subscript[a, 1])(x^2)+(Subscript[a, 2])x y+(Subscript[a, 3])(y^2)+(Subscript[a, 4])x+(Subscript[a, 5])y+(Subscript[a, 6])  where the coefficients could be exact rational numbers or numerical machine numbers . It is not a bad idea to plot the conic and typically we can see what sort of conic we have . An empty plot may indicate a curve with only complex points or a non - square free curve . Review Section 1.2 . Alternately if one wanted to, say, write a procedure to do this automatically one could first look for a real singular point, if a warning message is returned that there is an infinite solution set use sqFree to factor . Otherwise a single singular point could be found indicating two intersecting lines . If there is no singular point then look for infinite points, if there are none you have an ellipse, one generally indicates a parabola or possibly 2 parallel lines . Two infinite points give a hyperbola . These 3 latter cases are the interesting ones and will be discussed separately below .

Parabola

Let the infinite point be {u,v,0}. Then vx-uy=0 is the equation of the line through the infinite point and origin. Using A1=klRotation[vx-uy,x,x,y], the line b x-ay is rotated to the y-axis and the resulting conic is in the form g=a(x^2)+bx+c. Use

tangentRealPoints[f_,l_,x_,y_]:=With[{J=Det[D[{f,l},{{x,y}}]]},{x,y}/.NSolve[{f,J},{x,y},Reals]]

with tangentRealPoints[g,y,x,y] to find the vertex p and A2=PQTranslation[p,{0,0}] to move it to the origin. To get the equation in the form by=a(x^2), divide by b; note that this will not change the transformation.

The single FLT transformation A=A1.A2 is the transformation from the original conic to the conic y=a (x^2). Letting q=(1/(4a)), this parabola has focus {0,q} and directrix y=-q. Applying the FLT given by Inverse[A] gives the focus and directrix for the original conic.

f= 3.3244920269981044` -5.643417116045773` x+0.0836822330181507` x^2+4.251598647968187` y-0.7266575251177807` x y+1.5774888520713963` y^2

ipts=infiniteRealPoints[f,x,y]

ipt=ipts[[1]];
k=ipt.{-y,x,1}

A1=klRotation[k,x,x,y]

g=FLT[f,A1,x,y]

p=tangentRealPoints[g,y,x,y][[1]]

A2=pqTranslation[p,{0,0}];
h1=FLT[g,A2,x,y]

A=A2.A1

h=Expand[FLT[f,A,x,y]/-Coefficient[h1,y]]

q=1/4/Coefficient[h,x^2]

F=flt[{0,q},Inverse[A]]

d=FLT[y+q,Inverse[A],x,y]

plot1:=ContourPlot[{f==0,d==0},{x,-4,6},{y,-6,4}];
plot2:=ListPlot[{F}, PlotStyle->{Red}, PlotMarkers->{Automatic,Medium}];
Show[plot1,plot2]

Ellipse

We first put the ellipse in the standard form above; then use the classical formulas to find the eccentricity, a focus and the directrix; and then transform it back.

To find the standard form, we first find the center by intersecting two diagonals, that is lines connecting points with parallel tangents. After transforming the center to {0,0}, the default critical points of Chapter 3 give the closest and farthest points from the origin and hence the major and minor axes, which we can then rotate to the x and y axes, respectively, to get the standard form.

f=2x^2+2 x y+3 y^2-4x-2y-3;
ContourPlot[f==0,{x,-2,4},{y, -4,2}]

tp1=tangentRealPoints[f,x,x,y]

tp2=tangentRealPoints[f,y,x,y]

We use the line function from Chapter 1 to find the major and minor axes.

l1= line[#1,#2,x,y]&@@tp1
l2= line[#1,#2,x,y]&@@tp2

p={x,y}/.NSolve[{l1,l2}][[1]]

B1=pqTranslation[p ,{0,0}]

f1=FLT[f,B1,x,y]

pqTranslation[{1.`,0.`},{0,0}]

FLT[-3-4 x+2 (x^2)-2 y+2 x y+3 (y^2),pqTranslation[{1.`,0.`},{0,0}],x,y]

We need to rotate this ellipse to put the major axis on the x-axis.

cps=SortBy[criticalPoints[f1,x,y],Norm[#]&]

l3= Chop[line[#1,#2,x,y]&@@ Take[cps,-2]]

B2= klRotation[l3,y,x,y]
f2=FLT[f1,B2,x,y]
B=B2.B1;
FLT[f,B,x,y]

g=Expand[f2/f2[[1]]]

a=1/Sqrt[-Coefficient[g,x^2]]
b=1/Sqrt[-Coefficient[g,y^2]]

We note that a>b. If that is not the case, then we rotated incorrectly.

F1={Sqrt[a^2-b^2],0}

e= Sqrt[a^2-b^2]/a

d1=x-a/e

Show[ContourPlot[{f2==0, d1==0},{x,-3,3},{y,-3,3}],ListPlot[{F1}, PlotMarkers->Automatic]]

F=flt[F1,Inverse[B]]
d=FLT[d1,Inverse[B],x,y]

plot3:=ContourPlot[{f==0,d==0},{x,-2,4},{y,-3.5,2.5}];
plot4:=ListPlot[{F}, PlotStyle->{Red}, PlotMarkers->{Automatic,Medium}];
Show[plot3,plot4]

Hyperbola

We can proceed in much the same way as the ellipse, but the center will be the intersection of the two asymptotes. There will be only two real default critical points, but rotating the line through these to the x-axis (y=0) will give the desired standard form. The focus of ((x^2)/a)-((y^2)/b)=1 is {c,0}, the eccentricity is e=(c/a) and the directrix is x=(a/e)= ((a^2)/c) where c=Sqrt[(a^2)+(b^2)].

h=1-2x-3y-x^2+x y +2y^2;

Hint: this parabola transforms to

((x^2)/((1.9021130325903055`)^2))-((y^2)/((1.1755705045849454`)^2))=1,

where the focus is {0.558859,0} and the directrix is x=0.26509.

A converse to this problem is to start with the foci and find the conic. For the ellipse and hyperbola, the symmetry shows that there are actually two foci. The line containing them is the major axis and the conic is symmetric about this line. It turns out that the foci and one point on the conic determines the conic. Our goal is to find the equation from this information and the knowledge of whether this conic should be an ellipse or hyperbola. We should note that in the classical situation, the extra point q is where the semilatus rectum meets the ellipse or hyperbola, that is, a point q on the conic where the closest point on the major axis to q is a focus.

The additional point cannot be completely random. To be sure to get a conic with real coefficients, we restrict ourselves to good points. Given our desired foci (Subscript[p, 1]), (Subscript[p, 2]), a point q will be called good if it is not on the major axis but the closest point on the major axis l to q, call it (Subscript[q, c]), is between points (Subscript[p, 1]),(Subscript[p, 2]) and no farther from l than half the distance from (Subscript[p, 1]),(Subscript[p, 2]). This latter quantity is known as the linear eccentricity c.

Our technique is to transform our data (Subscript[p, 1]),(Subscript[p, 2]),q to standard position {0,c},{0,-c},(Subscript[q, 1]) using an FLT with matrix M. We then solve a nonlinear system to find a,b where, in the ellipse case, (c^(2 ))=(a^2)-(b^2) and (Subscript[q, 1]) satisfies the equation ((x^2)/(a^2))+((y^2)/(b^2))=1. In the hyperbola case, we use the equations (c^2)= (a^(2 ))+(b^2) and -((x^2)/(a^2))+((y^2)/(b^2))=1. We then apply the FLT with matrix (M^(-1)) to the latter equation to get our desired equation.

The details are similar to the earlier calculations in this section, but now we have a function to do this automatically. In this function, the default case is the ellipse case, which we will consider in the next chapter, an option conic→"hyperbola" will give the hyperbola case. The function will check to make sure for the given (Subscript[p, 1]),(Subscript[p, 2]) that q is good.

p1={1,3};
p2={3,5};
q={2,5};
f=conicWithFoci[p1,p2,q,x,y]

ContourPlot[{f==0,line[p1,p2,x,y]==0}, {x,0,4},{y,2,6}, Epilog->{PointSize[Large],Red,Point[{p1,p2,q}]}]

The focus is a Euclidean concept only. When using Euclidean transformations for plotting conics, it is important that the x-range is the same as the y-range. Otherwise, plotting introduces an affine, non-Euclidean transformation that distorts the ellipse or hyperbola. For example, doing the same plot with a different scale makes the major axis not look like the major axis!

{ContourPlot[{f==0,line[p1,p2,x,y]==0}, {x,.7,3.5},{y,1,7}, Epilog->{PointSize[Large],Red,Point[{p1,p2,q}]}],
ContourPlot[{f==0,line[p1,p2,x,y]==0}, {x,.7,3.5},{y,1,7}, Epilog->{PointSize[Large],Red,Point[{p1,p2,q}]}, AspectRatio->2.1]}

7.3 Projective Conics and Parameterization

All nonsingular projective conics are equivalent under FLT. We can use this to parameterize any nonsingular conic, either by a rational or trigonometric parameterization.

A key to this, both here and in the coming sections, is the cTransform[f,p,x,y] function.

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

To find an FLT taking a nonsingular conic f=0 to the parabola y=(x^2), first find a real point p on the conic and use A1=CTransform[f,p,x,y]. The result of this transform will be a parabola with vertical axis f1=a (x^2)+b x+c-d y.

Next, find the vertex q using q=tangentrealpoints[f1,y,x,y][[1]], which will transform by A2=pqTranslation[p,{0,0}] to f 2=a (x^2)+b y. We now use A3=homothety[r,1],r=±Sqrt[Abs[a/b]] to get a (x^2)-ay, which is equivalent to (x^2)-y.

h= 1-2x-3y-x^2+x y +2y^2;
cps= criticalPoints[h,x,y]

p=cps[[1]]

A1=cTransform[h,p,x,y]
h1=FLT[h,A1,x,y]

tp= tangentRealPoints[h1,y,x,y]

A2= pqTranslation[tp[[1]],{0,0}]

f2= FLT[h,A2.A1,x,y]

r=Sqrt[Coefficient[f2,x^2]/Coefficient[f2,y]]

A3=homothety[r,-1]

f3=FLT[h,A3.A2.A1,x,y]

A=A3.A2.A1

Rational Parameterization

Once we have an FLT given by matrix A, which takes our curve to the parabola y-(x^2), we can parameterize the latter by F t={t,(t^2)}, and apply the FLT given by Inverse[A] to get our rational parameterization.

A technical detail to note is that the function flt expects a pair of real numbers rather than variables.

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

F = flts[{t,t^2}, Inverse[A]]

ParametricPlot[F, {t,-10,10},Frame->True, Axes->False]

Limit[F,t->Infinity]

It is interesting to note that we have a converse.

If a parameterized curve has a parameterization of the form

F={(((Subscript[a, 1,1])t+(Subscript[a, 1,2])(t^2)+(Subscript[a, 1,3]))/((Subscript[a, 3,1])t+(Subscript[a, 3,2])(t^2)(Subscript[a, 3,3]))), (((Subscript[a, 2,1])t+(Subscript[a, 2,2])(t^2)+(Subscript[a, 2,3]))/((Subscript[a, 3,1])t+(Subscript[a, 3,2])(t^2)(Subscript[a, 3,3])))},

then letting t=x,(t^2)=y , we see this parameterization is in the form of an FLT with matrix

A=(({{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[a, 13])}, {(Subscript[a, 21]), (Subscript[a, 22]), (Subscript[a, 23])}, {(Subscript[a, 31]), (Subscript[a, 32]), (Subscript[a, 33])}})).

This curve can be obtained in the form FLT[(x^2)-y,A,x,y] and is generally a conic.

F={3t-t^2+1,-2t+t^2-2}/(1+t+t^2)

A=N[{{3,-1,1},{-2,1,-2},{1,1,1}}];

g=FLT[x^2-y,A,x,y]

P1:=ParametricPlot[F,{t,-10,10},PlotStyle->Black];
P2:=ContourPlot[g==0,{x,-4,3},{y,-3,3}];

{P1,P2,Show[P2,P1]}

Trigonometric Parameterization

An FLT that takes the parabola y=(x^2) to the circle (x^2)+(y^2)=1 is given by:

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

Therefore, start with any nonsingular conic and follow the transformation to y-(x^2)  by P2C. We get an FLT transformation from the conic to the standard unit circle. The unit circle can be parameterized by F={Cos[t],Sin[t]}. Applying the inverse transformation parameterizes our original conic by cosine and sine.

h =1-2 x-(x^2)-3 y+x y+2 (y^2);

A={{0.17152143109999518`,-1.4700447703390642`,1.4019961464680566`},{1.8134718990760177`,-1.4123120357497518`,1.422889116329615`},{0.4696313544127659`,-0.3657437507250866`,0.8035408513311424`}};

FLT[h,P2C.A,x,y]

F=flts[{Cos[t],Sin[t]},Inverse[P2C.A]]

ParametricPlot[F,{t,0,2 Pi}]

7.4 Inflection Points

A point p is an inflection point of an algebraic curve if the intersection multiplicity at p of the curve and its tangent line is 3 or higher. Inflection points are important in elementary mathematics for understanding graphs of single variable polynomial functions, which can be considered a special case of an algebraic curve. Lines and conics do not have inflection points. We will see in the next section that inflection points play an important role in the theory of cubic curves. In general, however, there may be many inflection points, but they are hard to compute and do not play a significant role in real curve theory.

Since we are mostly interested in inflection points of conics, the most efficient method of finding inflection points is the classical one. Intersect the curve with its Hessian curve. The Hessian curve is best viewed as a projective curve, even if we are only interested in real inflection points, so the first step is to homogenize the equation of the curve f=0 using the following (from Section 5.2):

fh=homog[f,x,y,z]

Next, we calculate the Hessian matrix:

(({{(∂_(x,x))fh, (∂_(x,y))fh, (∂_(x,z))fh}, {(∂_(y,x))fh, (∂_(y,y))fh, (∂_(y,z))fh}, {(∂_(z,x))fh, (∂_(z,y))fh, (∂_(z,z))fh}}))

Classical textbooks, e.g. (Fischer, 2001), often give shortcuts for hand calculation since fh is a homogeneous polynomial, but with Mathematica, we can simply use D[fh,{{x,y,z},2 to get the Hessian matrix. We then do a rare use of the determinant in numerical math to get the Hessian curve.

f=x^3-5x^2+7x-y^2

fh=homog[f,x,y,z]

H=Det[D[fh,{{x,y,z},2}]]

inflH={x,y}/.NSolve[{f,H, z-1},{x,y,z},Reals]

ha=H/.{z->1}

Show[ContourPlot[{f==0,ha==0},{x,-1,4},{y,-4,4}],ListPlot[inflH,PlotStyle->Red, PlotMarkers->Automatic]]

Note that in this case, the Hessian is of degree 3. More generally, the degree of the Hessian curve is 3(d-2) where d is the degree of the original curve, so there can be as many as 3d(d-2) inflection points for a curve of degree d, but as in this example, many are complex and some may be infinite.

affineInflectionPoints[f_,x_,y_]:=Module[{H,fh,d,z},
fh=Homog[f,x,y,z];
H=Det[D[fh,{{x,y,z},2}]];
{x,y}/.NSolve[{fh,H,z-1},{x,y,z},Reals]
];

To find all projective inflection points, we could use the trick from Chapter 6 for finding projective singular points.

rinfl= 8.390338886809685` x+1.572094833417954` y-1.8964087979890598`;
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=homog[rinfl,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 = 1-4x^3+2y+y^2;
allInflectionPoints[f,x,y]

7.5 Cubics

As mentioned in Chapter 1, most discussions of cubic curves are curves of the form

g= u[x]- (y^2) where u=u[x] =(x^3)+ r (x^2) + s x + t

and r, s, t are integers. In this book, and with the Wolfram Language, we can use machine numbers as well.

({{Description, Type, Zeros of u[x], Example, Picture}, {one oval, nonsingular, one real, 2 complex, (y^2)=(x^3)+x, }, {two ovals, nonsingular, a<b<c  real, (y^2)=(x^3)-x, }, {oval and isolated point, singular, a=b<c r eal, (y^2)=(x^3)-(x^2), }, {node, singular, a<b=c  real, (y^2)=(x^3)+(x^2), }, {cusp, singular, a=b=c  real, (y^2)=(x^3), }})

Reducible cubics are either the union of three lines or are the union of a line and a nonsingular conic. If the coefficients are integers, use Factor to get this reduction; otherwise, find points and use aCurve for the equation. The main result of this section is to show that every irreducible exact or numerical cubic can be reduced to the form u[x]-(y^2) shown previously or to the even simpler Weierstrass normal form (y^2)=(x^3)+s x+t.

This is a classical result, but it is usually described theoretically only or using examples that have been purposely chosen so that there are obvious inflection points with rational coefficients, e.g. (Garrity, 2013).

f=newtonHyperbola[836,x,y]

ips=allInflectionPoints[f,x,y]

We can pick any of these points. If we had an infinite one, we would give preference to that one. Here we just pick the first.

B1=cTransform[f, ips[[1]],x,y]
f1=FLT[f,B1,x,y]

This is what we expect: the x(y^2), (x^2)y and (y^3) terms are missing. If any were present at this step, they would probably have very small coefficients, which we can use Chop to shorten.

k=Expand[Coefficient[f1,y]/Coefficient[f1,y^2]/2+y]

B2={{1,0,0},{Coefficient[k,x],1,k[[1]]},{0,0,1}}
f2=FLT[f1,B2,x,y]

This is already up to a constant in the form (y^2)=u[x] , so we could stop here if our only interest was in parameterizing the curve with an FLT. The next two steps are fairly easy, so we continue next with a homothety.

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

Getting rid of the (x^2) term can be done with a linear change-of-variables trick that goes back to Gerolamo Cardano in 1539.

B4={{1,0,Coefficient[f3,x^2]/3},{0,1,0},{0,0,1}}
f4=FLT[f3,B4,x,y]

{ContourPlot[f==0,{x,-25,20},{y,-20,20}],ContourPlot[f4==0,{x,-4,5},{y,-10,10}]}

B=B4.B3.B2.B1

F1=flts[{x,Sqrt[10.136815015201705` -8.926609369242879` x+1.0000000000000002` (x^3)]},Inverse[B]]

F2=flts[{x,-Sqrt[10.136815015201705` -8.926609369242879` x+1.0000000000000002` (x^3)]},Inverse[B]];

ParametricPlot[{F1,F2},{x,-3.44,40}, PlotRange->20, Frame->True, Axes->False]

We note that in the interval -3.44≤x≤40, there is an interval where the cubic under the radical is negative and a point where the denominator is zero. However, ParametricPlot takes care of this for us.

As noted previously, we could have also used f 2 and B=B2.B1.

Since there are no other choices, we can put these steps together in an almost black-box algorithm. This is difficult with numerical algorithms since unexpected numerical errors may crop up.

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

In this implementation, we calculate the inflection points separately and choose one. We also return not only the normal form and FLT B but a number called the j-invariant (which we explore below). The user could change the last line and only get the normal form, if preferred.

In general, for a nonsingular cubic in Weierstrass form (y^(2 ))=(x^3)+a x+b, the j-invariant is

j-1728((4(a^3))/(4(a^3)-27(b^2))).

weierstrassNormalForm[f,ips[[1]],x,y]

weierstrassNormalForm[f,ips[[2]],x,y]

weierstrassNormalForm[f,ips[[3]],x,y]

We get different normal forms and FLT, but the j-invariant is always 874.894. That is one reason we call it an invariant; it does not depend on the inflection point chosen.

A={{0.07485986157215319`,0.9693047762451474`,0.24950067263711073`},{0.6094077133922`,-0.4539225635800057`,-0.5247878943188331`},{-0.02373700112965471`,0.3027242271117272`,-0.048647756985034984`}};
g=FLT[f,A,x,y]

ipg=allInflectionPoints[g,x,y]

weierstrassNormalForm[g,ipg[[2]],x,y]

We still get the same j-invariant. Therefore, cubics that differ by a fractional linear transformation have the same j-invariant. The converse is partly true. Suppose we have two cubics in Weierstrass form with the same j-invariant:

(Subscript[C, 1]): (y^2)=(x^3)+ (Subscript[a, 1])x + (Subscript[b, 1])
(Subscript[C, 2]): (y^2)=(x^3)+ (Subscript[a, 2])x +(Subscript[b, 2])

Assuming (Subscript[b, 1]) is not zero, I claim that the FLT homothety[(Subscript[d, 1]),(Subscript[d, 2])] where

(Subscript[d, 1])=(([(Subscript[b, 2])/(Subscript[b, 1])])^(1/3))
(Subscript[d, 2 ])=-((Subscript[d, 1])^(3/2))

takes cubic (Subscript[C, 1]) to cubic (Subscript[C, 2]). I leave the easy case (Subscript[b, 1])=0, (Subscript[a, 1])≠0 to the reader.

Now if (Subscript[B, 1 ,])(Subscript[B, 2]) are FLT taking cubics f to (Subscript[C, 1]) and g to (Subscript[C, 2]), then the FLT Inverse[(Subscript[B, 2])]. homothety[(Subscript[d, 1]),(Subscript[d, 2])]. (Subscript[B, 1]) takes f to g.

One technical detail is that if (Subscript[b, 1]),(Subscript[b, 2]) are of opposite signs, then (Subscript[d, 1]) is negative, so (Subscript[d, 2]) is imaginary. Thus, the claim that having the same j-invariant implies the cubics are FLT equivalent works only over the complex numbers. In particular, the formulas for a, b real, b≠0,

(Subscript[C, 1]): (y^2)=(x^3)+ ax + b
(Subscript[C, 2]): (y^2)=(x^3)+ ax -b,

have the same j-invariant but are not equivalent via a real FLT because they are equivalent by homothety[-1,-i ].

It is also true that given any real number j (other than 0), there is a nonsingular cubic with that j-invariant. One can pick j and b to be arbitrary real numbers and use NSolve to solve for a in the single variable polynomial equation j(4(a^3)-27(b^2))-6912(a^3)=0. Because it is of degree 3 in a, one zero must be real.

7.6 Caley’s Theorem

For the curves in this chapter, a recurring problem is to find a parameterization. Although the techniques are somewhat different, we end this chapter with Caley’s theorem from 1865, which says that an irreducible curve of any degree with a maximal number of singularities, each of which is a double point, is a rational curve, that is, the curve has a rational parameterization. This is a special case of a more general theorem in modern algebraic geometry that says a genus 0 curve is a rational curve. The interesting thing about Caley’s proof is that it is constructive. Using the Wolfram Language, we can actually find this parameterization, following the argument in (Popescu-Pampu, 2016).

First, we say what we mean by a double point. It is a point p on a curve f=0 such that the intersection multiplicity with a random line through this point is 2. Since a regular point on a curve locally behaves like its tangent line, this also means that the intersection multiplicity of any curve g that intersects f at p at a regular point of g has intersection multiplicity 2 with f at p.

{Labeled[ContourPlot[y^2-x^3-x^2==0,{x,-1,1},{y,-1,1}],"node"], Labeled[ContourPlot[y^2-x^3==0,{x,-1,1},{y,-1,1}],"cusp"], Labeled[ContourPlot[(y-x^2)(y+x^4)==0,{x,-1,1},{y,-1,1}], "common tangent"],Labeled[ContourPlot[x y (x-y)==0,{x,-1,1},{y,-1,1}],"triple point"]}

l=RandomReal[{-1,1},2].{x,y}

{IntersectionMultiplicity[y^2-x^3-x^2,l,{0,0}],
IntersectionMultiplicity[y^2-x^3,l,{0,0}],
IntersectionMultiplicity[(y-x^2)(y+x^4),l,{0,0}],
IntersectionMultiplicity[x y (x-y),l,{0,0}]}

In particular, nodes and cusps are double points.

Caley’s theorem is that if an irreducible plane curve of degree d has (((d-1)(d-2))/2) double points, then it can be parameterized by rational parameters.

We know that an irreducible plane curve cannot have more than (((d-1)(d-2))/2) singularities. A curve with this many singularities has genus 0.

Let f=0 be our irreducible curve of degree d with (((d-1)(d-2))/2) double points. We start our proof/algorithm by forming a set 𝒮 consisting of the singularities and 2d-3 additional random points on f. These should not be too bunched up.

Expand[(d-1)(d-2)/2+2d-3]

This can be written as (((d+1)d)/2)-2, which is exactly one less than the number of points needed (see the end of Chapter 2) for interpolating a unique bivariate polynomial of degree d-1. We now choose (randomly but with some care) two points: (Subscript[q, 1]),(Subscript[q, 2]) not on f, and, in practice, not too near f. Then the sets 𝒮⋃{(Subscript[q, 1])} and 𝒮⋃{(Subscript[q, 2])} will produce two curves (Subscript[g, 1]),(Subscript[g, 2]) of degree d-1 by interpolation. The randomness should be enough to require (Subscript[g, 1])≠(Subscript[g, 2]). Nonetheless, both curves should contain all the points in 𝒮, but because they contain (Subscript[q, 1]) or (Subscript[q, 2]) and are of smaller degree, they are also distinct from f.

For each real number, let t, (Subscript[g, t])=(Subscript[q, 1])+t (Subscript[q, 2]). This is, again, a curve of degree d-1 containing the points of 𝒮.

Our second serendipity is that by applying Bézout, each (Subscript[g, t]) intersects f in d(d-1) points counting intersection multiplicities.

Expand[2((d-1)(d-2)/2)+2d-3]

Note that this is exactly one less than the Bézout total. The intersection of f, (Subscript[g, t]) contains exactly one point not in 𝒮, so we can name this point Φ[t], but since t is a real variable, a small change in t gives a small perturbation in (Subscript[g, t])=(Subscript[q, 1])+t (Subscript[q, 2]), so the solution set will only have a small change, but 𝒮 does not change, so instead, Φ[t] has a small change. Thus, Φ[t] will be a continuous function giving a continuous parameterization of at least part of the curve f=0.

f=Expand[3 ((x^2)+(y^2))^2+8 x (-(x^2)+3 (y^2))+6 ((x^2)+ (y^2))-1]

cp=criticalPoints[f,x,y]

S=cp[[{1,2,5}]]

rl=RandomReal[{-1,1},2].{x-S[[1,1]],y-S[[1,2]]};
intersectionMultiplicity[f, rl,S[[1]],1.*^-10]

S=Join[S,cp[[{3,4,6}]],{x,y}/.NSolve[{x-.5,f},{x,y},Reals]]

S1={{1.1354546845775721`,-0.3327801966326215`},{-0.27241452246256326`,1.1514296886703472`}};

g1=aCurve[Append[S,S1[[1]]],x,y]
g2=aCurve[Append[S,S1[[2]]],x,y]

{ContourPlot[{f==0,g1==0},{x,-1.2,1.2},{y,-1.2,1.2}], ContourPlot[{f==0,g2==0},{x,-1.2,1.2},{y,-1.2,1.2}]}

g:=Expand[g1+# g2]&;

g[.5]

S
sol=DeleteDuplicates[{x,y}/.NSolve[{f,g[.5]}],Norm[#1-#2]<1.*^-9&]

Show[ContourPlot[{f==0,g[.5]==0},{x,-1.2,1.2},{y,-1.2,1.2}],ListPlot[{sol[[-1]]},PlotMarkers->Automatic]]

ϕ[t_,S_]:=Module[{sol},
sol={x,y}/.NSolve[{g[t],f}];
Select[sol,Total[Boole[Table[Norm[p-#]<1.*^-10,{p,S}]]]==0&][[1]]];

L=Table[ϕ[t,S],{t,0,2,.1}];

Show[ContourPlot[f==0,{x,-2,2},{y,-2,2},ContourStyle->Cyan],ListLinePlot[L, PlotStyle->Black]]

This example shows we have a parameterization of the hypocycloid, but we need to show that ϕ is a rational function of t. To do so, we must introduce the classical resultant.

Prior to the twenty-first century, solving nonlinear polynomial systems was, in general, an ad-hoc process, but there was an established method for 2×2 systems using the resultant. The reader who is unfamiliar may look this up in many references in the theory of equations, algebra or algebraic geometry. The resultant is implemented in Mathematica and, in short, if f and g are bivariate polynomials in x and y, then Resultant[f,g,y] is a polynomial in the one variable x whose zeros are all the x-coordinates of solutions of the system. By substituting these back in f and g and solving for y, one can get all the solutions to the system. Typically, the resultant is the product of degree( f  ) and degree( g), which almost proves Bézout’s theorem except for the problem of infinite solutions.

The key to Caley’s method is the well-known fact in elementary algebra: for a monic polynomial p in the one variable x of degree d, the coefficient of (x^(d-1)) is the sum of the zeros. More generally, the sum of the zeros is -Coefficient[p, (x^((d - 1)))]/Coefficient[p, (x^d)] after dividing by the leading coefficient. However, as we have seen above, all but one zero of the intersection of (Subscript[g, t]) and f  are known, so all but one zero of the resultants Resultant[f,(Subscript[g, t]) x] and Resultant[f,(Subscript[g, t]),y] are known. Finding the last one just requires an addition and a subtraction. As we see in the example, however, the coefficients are functions of t, so we get a rational expression in t. This is our rational parameterization!

TSx=-Total[Flatten[Take[S,All,1]]]
TSy=-Total[Flatten[Take[S,All,-1]]]

Rx=Resultant[f,SetPrecision[g[t],20], y];
Ry=Resultant[f,SetPrecision[g[t],20],x];

ϕx=TSx- Coefficient[Rx,x^11]/Coefficient[Rx,x^12];
ϕy=TSy-Coefficient[Ry,y^11]/Coefficient[Ry,y^12];
dϕx=N[Together[ϕx],6]
dϕy=N[Together[ϕy],6]

lx={dϕx[[1]],dϕy[[1]]}

Show[ContourPlot[f==0,{x,-2,2},{y,-2,2},ContourStyle->Cyan],ParametricPlot[{ϕx,ϕy},{t,-10,10}, PlotStyle->Black], ListPlot[{lx},PlotMarkers->Automatic]]