Mathematica 12.3 ,13             December 2021

Global functions for Surfaces

Note the naming conventions:   All global functions defined in the Surface Book begin with a lowercase letter, compound names will capitalize first letters  of subsequent words, (camel casing).  This avoids confusion with built in Mathematica functions. Also curve functions with polynomial and/or point arguments will  end in  2D, 3D,   or MD depending on whether they work in 2,3 or all dimensions.  Surface functions may end in S, NS or RS.  This makes clear what  the  arguments are and distinguishes these functions from my Plane Curve functions so both sets can be initialized together without conflict, however most plane curve functions that you may need are contained here with 2D designation.  Note that functions with suffix MD take  variables as a list, but members of the list should be atomic variables, e.g. not X[[2]] but possibly x[2].   CAUTION:  This is still a work in progress and the version on my website may change without warning eg. note change in hBasis, hbasisQ and FLTMD effective 11/2019.  Note also bug fix in 62  nDivideMD  in January 2021.  A new procedure for intersecting lines given parametrically in space and some related functions were added February 2021.  ClosestRealPointMD added April 2021.

Before using  functions Evaluate Initialization Cells.  To test functions Evaluate Notebook. There should be no Mathematica errors but there will be some warnings from the functions defined here.  Caution-some functions may not run correctly on earlier versions of Mathematica than 12.3

The author makes no representations, express or implied, with respect to this documentation or software it describes, including, without limitation, any implied warranties of merchantability, interoperability or fitness for a particular purpose, all of which are expressly disclaimed.Use of Mathematica and other related software is subject to the terms and conditions as described at www.wolfram.com/legal.
In addition to the forgoing, users should recognize that all complex software systems and their documentation contain errors and omissions.Barry H.Dayton, and Wolfram Research a) shall not be responsible under any circumstances for providing information or corrections to errors and omissions discovered at any time in this  software; b) shall not be liable for damages of any kind arising out of the use of (or inability to use) this software; c) do not recommend the use of the software for applications in which errors or omissions could threaten life, or cause injury or significant loss.

Original material in this notebook is covered by Creative Commons license by - nc - sa/3.0.
https://creativecommons.org/licenses/by-nc-sa/3.0/

1; A good (Overscript[default,  ]) tolerance, a constant.

dTol:=1.*^-12;

dTol

2; Total degree of a polynomial in X.   Rounding and chopping polynomials.

tDegMD[f_,X_]:=Max[Total/@Keys[CoefficientRules[f,X]]];

tDegMD[x^3-x y^2-y z+3+0 z^5+1.*^-15 x^4,{x,y,z}]

Note monomials with coefficient 0 are ignored however monomials with very small coefficients are counted.  If f is a numerical polynomial, eg. at least one coefficient is a decimal, you can use the following to ignore monomials with very small coefficients.

tDegMD[Chop[x^3-x y^2-y z+3+0 z^5+1.*^-15 x^4,dTol],{x,y,z}]

3

Occasionally we need to round polynomials .  Typically Mathematica will do this automatically however will leave coefficients as machine numbers.  If the coefficients appear to be integers and you want exact integer coefficients  use  tolerance 1

roundPolyMD[P_,X_,tol_]:=FromCoefficientRules[Normal[Round[<|CoefficientRules[P,X]|>,tol]],X];

P=9.000000000000002` -9.999999999999957` (x^2)+1.0000000000000049` (x^4)-9.999999999999712` (y^2)+1.9999999999999667` (x^2) (y^2)+0.999999999999971` (y^4)+5.999999999999873` (z^2)+2.000000000000007` (x^2) (z^2)+1.9999999999999287` (y^2) (z^2)+0.9999999999999508` (z^4);

roundPolyMD[P,{x,y,z},1.*^-5]

roundPolyMD[P,{x,y,z},1]

While the built in Round will not work directly with polynomials note that  Chop will . Caution, chop may change the degree but  roundPolyMD usually will not.

Chop[.3-.0001x+1.3*1*^-12 x^2]
roundPolyMD[0.3` x-0.0001` (x^2)+1.3000000000000001`*^-12 (x^3),{x},10.^-5]

3;  linearSetMD , including lines.  P  is a list of  points in (ℝ^n).  The smallest linear set of (ℝ^n) containing P is given.  The first point must be affine, i.e. length n, other points may be affine or projective of length n+1.  To get the equations defining the linear vector space spanned by these points as vectors add the point of length n  {0,…,0} this vector space may have dimension one higher than the linear set without the zero point.. Note the actual dimension is n minus the number of linear equations given.  This routine returns the empty list when the linear set has dimension n which means it is all of (ℝ^n).

linearSetMD[P_,X_]:=Module[{n,XX,t,H,M,U,S,V,T,eq0,eq1,r,rk},
n=Length[X];
XX=Append[X,t];
r=Length[P];
If[Length[P[[1]]] ≠ n,Echo[n,"First point must be affine, length "];Abort[]];
M=Table[Switch[Length[P[[i]]],n,Append[P[[i]],1],n+1, P[[i]],_,Echo[i,"illegal point "];Abort[]],{i,r}];
{U,S,V}=SingularValueDecomposition[N[M]];
rk=Length[Select[Diagonal[S],#>1.*^-10&]];
If[rk>n,Return[{}]];
T=Take[Transpose[V],rk-n-1,n];
eq0=T.X;
eq1=eq0/.Thread[X->P[[1]]];
eq0-eq1]

lineMD[p_,q_,X_]:=linearSetMD[{p,q},X];

linearSetMD[{{1,0,0},{0,1,1},{2,1,3}},{x,y,z}]
linearSetMD[{{1,0,0},{0,1,1},{2,1,3},{.5,.5,.5},{1,1,2}}, {x,y,z}]

These two linear sets have dimension 3-1=2

linearSetMD[{{0,0,0},{1,0,0},{0,1,0},{0,0,1}},{x,y,z}]

This linear set obviously has dimension 3.

lineMD[{1,1,1},{1,2,3,0},{x,y,z}]

Find the implicit equations of a rational parametric line in (ℝ^3) .  A simple way to find the implicit equations is to calculate 2 points and use lineMD.   But sometimes to get more accuracy or to find  rational points on a surface if the 2 points are rational we would like an equation system with rational coefficients.   But lineMD returns floating point numbers as do the methods in section 1.3.  A simple routine is here.

ratLine3D[p_,q_]:=Module[{form, formp,formq,sol},
form={ x-a  y+b,x-c z+d};
formp=form/.Thread[{x,y,z}->p];
formq=form/.Thread[{x,y,z}->q];
sol=Solve[formp==0&&formq==0][[1]];
form/.sol]

Note that it is assumed that the variables are x, y, z and that x is a parameter, meaning the two points p, q have distinct first component.  If not rename the variables, run then name them back again.  It is somewhat surprising that the equation solved appears to be underdetermined, but Solve apparently needs the extra variable. Anyway we only need one solution so if the Solve returns several we  only use the first.  Here is an example:

p={-(14/15),-(17/15),0};
q={(1/13),-(11/13),(11/13)};

l=ratLine3D[p,q]

Test: Note that r1 p + r2 q will be in the line through p,q for any r1+r2=1

r=3/7 p+ 4/7q

l/.Thread[{x,y,z}->r]

4; Quick test for line or plane in (R^3),P a list of 4 affine points, if possibly collinear a line is returned, if possibly planar a plane is returned, otherwise {} is returned .

planar3D[P_]:=Module[{v,w,n,c,l,res},
If[Length[P]<4 ||Dimensions[P][[2]]≠3, Print["need 4 points in (R^3)"];Abort[]];
n=Cross[P[[2]]-P[[1]],P[[3]]-P[[1]]];
If[Norm[n]<.0001,
l=lineMD[P[[2]],P[[3]],{x,y,z}],
c=n.P[[1]];
l=n[[1]]x+n[[2]]y+n[[3]]z-c];
res=Echo[l/.Thread[{x,y,z}->P[[4]]],"Residue = "];
If[Abs[res]>.0001, l={}];
l];

planar3D[{{1,0,0},{0,1,0},{0,0,0},{1, 1,1}}]

5; Tangent vector in 3-space at point p.

Options[tangentVector3D]={tol->1.*^-8};
tangentVector3D[{f_,g_},p_,{x_,y_,z_},OptionsPattern[]]:=Module[{n1,n2,bi},
If[Norm[{f,g}/.Thread[{x,y,z}->p]]>OptionValue[tol],Echo[p,"not a point "];Return[Fail]];
n1={D[f,x],D[f,y],D[f,z]}/.Thread[{x,y,z}->p];
n2={D[g,x],D[g,y],D[g,z]}/.Thread[{x,y,z}->p];
bi=N[Cross[n1,n2]];
If[Norm[bi]<.0001,Echo[p,"No tangent vector at "]; bi,Normalize[bi]]]

F={y-x^2,z-x^3};
tangentVector3D[F,{1,1,1},{x,y,z}]

6; Critical points of naive curve or surface in 3 - space using random linear optimization function, see below, will return randomly different points each run.

criticalPoints3D[{f_,g_},{x_,y_,z_}]:=Module[{J,ob},
ob=RandomReal[{.7,1.3},3].{x^2,y^2,z^2};
J=D[{f,g,ob},{{x,y,z}}];
{x,y,z}/.NSolve[{f,g,N[Det[J]]},{x,y,z},Reals]]

F={x^2+.25y^2+.81z^2-1,x+y-z};
cpF=criticalPoints3D[F,{x,y,z}]

Critical Points of arbitrary algebraic set, not necessarily a curve, using optimization with a random linear objective function similar to an algorithm by Wang and Bican.  F is a set  of polynomials, X is the list of variables.  Can have more or fewer equations than variables, but if the original has k equations in n variables the resulting square system has size n+k.   May not find points in all components, in particular will usually not find isolated points of curves.  Because of the random objective function different runs may give different numbers of critical points, run several times. If you are persistently not getting enough solutions add the option solutions->Complexes and look for real points.  This may not find all components, eg. isolated points.

Options[criticalPointsMD]={solutions->Reals};
criticalPointsMD[F_,X_,OptionsPattern[]]:=Module[{uv,n,k,wbg,b},
n=Length[X];
k=Length[F];
uv=Table[u[i],{i,k}];
b=RandomReal[{-1,1},{n,1}];
Echo[X.b-RandomReal[{-1,1}],"Objective Function "];
wbg=Flatten[Expand[uv.Grad[F,X]-b]];
X/.NSolve[Join[F,wbg],Join[X,uv],OptionValue[solutions]]]

F={x^2+.25y^2+.81z^2-1,x+y-z};
criticalPointsMD[F,{x,y,z}]

C4={x+y+z+w, x y+y z+z w+w x, x y z+ y z w+ z w x+ w x y, x y z w-1};

cps=criticalPointsMD[C4,{w,x,y,z},solutions->Complexes]

Chop[cps,10*^-6]

Check that {-1,1,-1,1} is a real point on C4.  Again, rerun for more points.

Milnor critical points

Based on Milnor' s 1963 paper “Betti numbers of real varieties”.  Again note randomness.  You may get an empty set on some runs, run this several times to get the largest set.  F is a single equation in 3 variables X.

milnorCritPtsNS[F_,X_]:=Module[{A,G,G1,G2,sol},
A=Orthogonalize[RandomReal[{-1,1},{4,4}]];
G=FLT3D[{F},A,X][[1]];
G1=-D[G,X[[1]]]/D[G,X[[3]]];
G2=-D[G,X[[2]]]/D[G,X[[3]]];
sol=NSolveValues[{G,G1,G2},X,Reals,WorkingPrecision->2];
fltMD[#,Inverse[A]]&/@sol
]

cps=milnorCritPtsNS[x^3+y^3+z^3-1,{x,y,z}]

Show[ContourPlot3D[x^3+y^3+z^3==1,{x,-2,2},{y,-2,2},{z,-2,2},Mesh->None],Graphics3D[{Red,Table[Ball[cps[[i]],.1],{i,Length[cps]}]}]]

This function can find singular points, eg

milnorCritPtsNS[ x^2+y^2-z^2,{x,y,z}]

7;  Closest point
 from a point in 3-space to a naive curve.  Here {f,g}  is the naive curve, p is the point.

closestPoint3D[{f_,g_},p_,{x_,y_,z_}]:=Module[{J,sol},
J=D[{f,g,(x-p[[1]])^2+(y-p[[2]])^2+(z-p[[3]])^2},{{x,y,z}}];
sol={x,y,z}/.NSolve[{f,g,N[Det[J]]},{x,y,z},Reals];
MinimalBy[sol,Norm[#-p]&][[1]]
]

F={x^2+.25y^2+.81z^2-1,x+y-z};
closestPoint3D[F,{0,0,0},{x,y,z}]

From Real point to a real curve.  Here F is an  arbitrary real curve, p is a real point.  This uses the powerful built in minimizer FindMinimum.  Use option distance->True to get the distance from point to curve.

Options[closestRealPointMD]={distance->False};
closestRealPointMD[F_,p_,X_,OptionsPattern[]]:=Module[{con,k,eq,eqs,obj,s,min,q},
k=Length[F];
s=Length[X];
con=ToString[F[[1]]==0,InputForm];
Do[eqs=ToString[F[[i]]==0,InputForm]; con=con<>"&&"<>eqs,{i,2,k}];
con=ToExpression[con];
obj=Total[Table[(X[[i]]-p[[i]])^2,{i,s}]];
{min,q}=FindMinimum[{obj,con},X];
If[OptionValue[distance],Echo[min,"distance"]];
X/.q]

F={(y^2)-x z,x y-z,(x^2)-y};
closestRealPointMD[F,{1,2,3},{x,y,z}]

You might want to check to see if answer is on curve, if there is no real solution this function  may return a point not on the curve.  The following is typical an actual real solution:

q=closestRealPointMD[F,{1,2,3},{x,y,z}, distance->True]
F/.Thread[{x,y,z}->q]

Here there is no real solution.

q=closestRealPointMD[{x^2+y^2+5},{1,2},{x,y}]
(x^2+y^2+5)/.Thread[{x,y}->q]

For complicated cases there may be many warning messages!

8;  Find path (piecewise linear)  between two points on an arc of a naive curve in direction of tangent vector.  s is step size.  The maximum number of points can be changed by option maxit.  To change the direction of the path one can use option dir→-1

Options[pathFinder3D]={maxit->30,tol->1.*^-8,dir->1};
pathFinder3D[{f_,g_},p_,q_,s_,{x_,y_,z_},OptionsPattern[]]:=Module[{k,p0,p1,tv1,tv,L},
p0=p;
L=Reap[Sow[p];
k=0;
While[Norm[q-p0]>2s && k<OptionValue[maxit],
tv1=OptionValue[dir]*tangentVector3D[{f,g},p0,{x,y,z},tol->OptionValue[tol]];
If[tv1.(q-p0)>0,tv=tv1,tv=-tv1];
p0=closestPoint3D[{f,g},p0+s*tv,{x,y,z}];
Sow[p0];
k++];
If[k≥OptionValue[maxit],Print["Warning, iteration limit reached"]];
Sow[q]];
L[[2,1]]];

F={x^2+.25y^2+.81z^2-1,x+y-z};
pathFinder3D[F, cpF[[1]],cpF[[3]],.4,{x,y,z}]

9; forms of monomial in n variables  (3D case no longer used).

formMD[f_,k_,X_]:=FromCoefficientRules[Select[CoefficientRules[f,X],Total[#[[1]]]==k&],X];
maxFormMD[f_,X_]:=formMD[f,tDegMD[f,X],X];

f = x^2 + y^3 - z^3 - x y z^2 + x^2 y^2;
formMD[f, 3, {x, y, z}]
maxFormMD[f, {x, y, z}]

12;  Infinite real points for naive curve in (ℝ^3).  Does not generalize, see 66 for general case.

infiniteRealPoints3D[{f_,g_},{x_,y_,z_}]:=Module[{sol},
sol={x,y,z}/.NSolve[{maxFormMD[f,{x,y,z}],maxFormMD[g,{x,y,z}],x^2+y^2+z^2-1},{x,y,z},Reals];Append[#,0]&/@Tally[sol,Norm[#1+#2]<.0001&][[All,1]]]

F={x^2-y^2-1, x+y+z-1};
infiniteRealPoints3D[F,{x,y,z}]

13; Finds singular real points of naive 3-space curve.   Does not generalize.

singularRealPoints3D[{f_,g_},{x_,y_,z_}]:=Select[Tally[criticalPoints3D[{f,g},{x,y,z}], Norm[#1-#2]<.0001&],#[[2]]>1&][[All,1]]

F={y^2-x z-z,z-x^2+1}
singularRealPoints3D[F,{x,y,z}]

14; Homogeneous exponent lists of degree d in n variables.

hExpsMD[n_,d_]:=Module[{hps},
hps[0]={Table[0,{n}]};
hps[m_]:=hps[m]=DeleteDuplicates[Flatten[Table[ReplacePart[p,i->(p[[i]]+1)],{p,hps[m-1]},{i,n}],1]];
hps[d]];
hExps[n_,d_]:=hExpsMD[n,d];

hExpsMD[4,2]

Actual Monomials

mhExpsMD[d_,X_]:=Take[mExpsMD[d,X],-Length[hExpsMD[Length[X],d]]]

mhExpsMD[2,{x,z,y}]

15; Affine exponent list of degree d in n variables.

expsMD[n_,d_]:=Drop[hExpsMD[n+1,d],None,1];
exps[n_,d_]:=expsMD[n,d];

expsMD[3,2]

16  Gives CoefficientRules as an association, returning also the degree.  Missing terms are given coefficient zero.  Used in constructing Macaulay and Sylvester matrices

fAssocMD[f_,X_]:=Module[{A,d,n,FA},
n=Length[X];
A=Association[CoefficientRules[f,X]];
d=Max[Table[Total[p],{p,Keys[A]}]];
FA=Association[Table[If[MissingQ[A[p]],p->0,p->A[p]],{p,exps[n,d]}]];
{FA,d}]

FA=fAssocMD[1+2x y-3z ,{x,y,z}]

17; Gives homogeneous coefficient rules as association adding in zero coefficients.  Homogeneous degree given.

hAssocMD[h_,X_]:=Module[{A,d,n,FA},
n=Length[X];
A=Association[CoefficientRules[h,X]];
d=Max[Table[Total[p],{p,Keys[A]}]];
FA=Association[Table[If[MissingQ[A[p]],p->0,p->A[p]],{p,hExps[n,d]}]];
{FA,d}]

hAssocMD[w+2 x y-3 z w,{x,y,z,w}]

18; The shift operation for constructing Sylvester and Macaulay matrices.  Note q is a list of length n.

shiftFAMD[FA_,q_,d_]:=Module[{S,K,n}, 
K=Keys[FA];
n=Length[K[[1]]];
S=Association[Table[p+q->FA[p],{p,K}]];
Association[Table[If[MissingQ[S[p]],p->0,p->S[p]],{p,exps[n,d]}]]];

shiftFAMD[FA[[1]],{1,0,1},4]

19; Sylvester matrix  of order m for one polynomial.

sylMD[f_,m_,X_]:=Module[{FA,d},
n=Length[X];
{FA,d}=fAssocMD[f,X];
If[d>m,Print["Degree error in syl"];Abort[]];
Table[Values[shiftFAMD[FA,q,m]],{q,exps[n,m-d]}]];

sylMD[1+x y-y^2,3,{x,y,z}]//MatrixForm

20; Sylvester matrix for list of polynomials of order m.

sylvesterMD[F_,m_,X_]:=Flatten[Table[sylMD[F[[i]],m,X],{i,Length[F]}],1];

sylvesterMD[{x-1,y+2},2,{x,y}]//MatrixForm

21;  Macaulay matrix of order m for one polynomial at point p.  Important note: if any row is all zero that row is eliminated.

macaMD[f_,m_,p_,X_]:=Module[{M,fp,FA,d,n},
fp=Expand[f/.Thread[X->X+p]];
n=Length[X];
{FA,d}=fAssocMD[fp,X];
M=Table[Values[shiftFAMD[FA,q,m]],{q,expsMD[n,m-1]}];
Select[M,AnyTrue[#,Abs[##]>0&]&]]

macaMD[1+x y-y^2,2,{0,0,0},{x,y,z}]//MatrixForm

22; Macaulay Matrix of order m at point p.  Important note: Because zero rows are eliminated one can not predict the size of a Macaulay matrix in advance.

macaulayMD[F_,m_,p_,X_]:=Flatten[Table[macaMD[F[[i]],m,p,X],{i,Length[F]}],1];

macaulayMD[{x-1,y+2},2,{1,-2},{x,y}]//MatrixForm

23; Numerical reduced row echelon form of matrix M with tolerance epsilon.  Gives list, first entry is list of pivot columns

nrref[M_, eps_]:= Module[{p,P,j, R,mn,n,r,s, U,S,V}, 
{U,S,V}= SingularValueDecomposition[N[M],Tolerance->eps];
r=Length[Select[Diagonal[S],#>0&]]; (* rank *)
R = Take[Transpose[V],r];   (* row space of M *)
 mn=Dimensions[R];
n=1;
While[Norm[Take[R,All,{n}]]<eps,n++];
p = {n};
For[j=n,j>0,j++,
If[mn[[1]]<=Length[p],Break[],
p=Append[p,j];
P=Check[R[[All,p]],Abort[]];  
s=Length[Select[SingularValueList[N[P],Tolerance->eps],#>0&]];
If[s<Length[p], 
p=Drop[p,-1];, Null];
];
];
P=R[[All,p]];
{p,Chop[Inverse[P].R]} ]

M={{1,3,-4,3},{-2,1,3,-4},{-.99, 4.01,-1.01,-1}};
rrM=nrref[M,.001];
rrM[[1]]
rrM[[2]]//MatrixForm

24  Approximate reverse row echelon form, first entry is pivot list.

arref[M_, eps_]:= Module[{p,P,j,r,s, R,mn,n,U,S,V}, 
{U,S,V}= SingularValueDecomposition[N[M],Tolerance->eps];
r=Length[Select[Diagonal[S],#>0&]]; (* rank *)
R = Take[Transpose[V],r];   (* row space of M *)
 mn=Dimensions[R];
n=mn[[2]];
While[Norm[Take[R,All,{n}]]<eps,n--];
p = {n};
For[j=n-1,j>0,j--,
If[mn[[1]]<=Length[p],Break[],
p=Prepend[p,j];
P=Check[R[[All,p]],Abort[]];   
s=Length[Select[SingularValueList[N[P],Tolerance->eps],#>0&]];
If[s<Length[p], 
p=Drop[p,1];, Null];
];
];
P=R[[All,p]];
{p,Chop[Check[Inverse[P],Abort[]].R]} ]

M={{1,3,-4,3},{-2,1,3,-4},{-.99, 4.01,-1.01,-1}};
arrM=arref[M,.0001];
arrM[[1]]
arrM[[2]]//MatrixForm

25; Dual of Sylvester or Macaulay matrix.

dualMatrix[A_,tol_]:=Module[{ns,r,c, U,S,V},
If[Length[Dimensions[A]]!=2,Echo["DualMatrix: not a matrix"];Abort[]];
c=Dimensions[A][[2]];
If[c==0,Return[IdentityMatrix[First[Dimensions[A]]]]];
{U,S,V}=SingularValueDecomposition[N[A],Tolerance->tol];
r=Length[Select[Diagonal[S],#>0&]]; (* numerical rank *)
Take[V,All, r-c]]

DS=dualMatrix[sylvesterMD[{x^2,y^2},2,{x,y}], dTol]

26; Local dual matrix or "left dual"

localDualMatrix[A_,tol_]:=With[{dim=Dimensions[A]},
If[dim!=2,Echo["localDual: not a matrix"];Abort[]];
If[dim[[1]]==0,Return[IdentityMatrix[dim[[2]]]]];
Transpose[dualMatrix[Transpose[A],tol]]];

localDualMatrix[DS, dTol].mExpsMD[2,{x,y}]

A utility function is matrixrank .

matrixrank[M_,tol_]:=Module[{s,k,l},
s=SingularValueList[N[M],Tolerance->0];
If[s[[1]]<tol,Return[0]];
l=Length[s];
s=s/s[[1]];
k=1;
While[k<=l,If[s[[k]]<tol,Return[k-1],k++]];
k-1];

Note that one may with to run this with a very loose tolerance to compare with tight tolerance, agreement means the matrix is probably well conditioned.

Here is a nice application to elementary linear algebra,  the intersection of 2 vector subspaces as done in section 2.5.3 .  The actual vector spaces are the row spaces of V,W.  Both V,W have same number of columns.

vectorSpaceIntersection[V_,W_,tol_]:=Module[{dV,dW,dd},
dV=dualMatrix[V,tol];
dW=dualMatrix[W,tol];
dd=Join[dV,dW,2];
localDualMatrix[dd,tol]]

Example :

V={{2,0,0,-2,1},{-2,1,3,-2,0},{1,2,3,4,5}};
W={{2,-3,2,-3,-3},{0,-7,1,2,-9},{4,1,-8,5,4}};
VW=vectorSpaceIntersection[V,W,dTol]

Check :

{matrixrank[V,dTol],
matrixrank[W,dTol],
matrixrank[Join[V,W],dTol],
matrixrank[VW,dTol],
matrixrank[Join[V,VW],dTol],
matrixrank[Join[W,VW],dTol]}

The first and last two 3' s show that VW is contained in both V and W, hence in the intersection.  But the 4 shows V ≠W so the intersection can have at most dimension 2 which VW has.  Hence VW is the intersection.

This function can be extended to any number of subspace,s for example

vectorSpaceIntersection3[V_,W_,U_,tol_]:=Module[{dV,dW,dU,dd},
dV=dualMatrix[V,tol];
dW=dualMatrix[W,tol];
dU=dualMatrix[U,tol];
dd=Join[dV,dW,dU,2];
localDualMatrix[dd,tol]]

27;  All monomials in variables X of degree d or less.

mExpsMD[d_,X_]:=Module[{n},
n=Length[X];Table[FromCoefficientRules[{p->1},X],{p,expsMD[n,d]}]];

mExpsMD[3,{x,y}]

Homogeneous version

mhExpsMD[d_,X_]:=Table[FromCoefficientRules[{r->1},X],{r,hExpsMD[Length[X],d]}];

mhExpsMD[3,{x,y,z}]

28;  Calculates Hilbert function from nrref data.

Options[hilbertFunctionMD]={diff->False};
hilbertFunctionMD[p_,m_,n_,OptionsPattern[]]:=Module[{h},h=Table[Binomial[d+n-1,n-1]-Length[Select[p,Binomial[d+n-1,n]< # ≤Binomial[d+n,n]&]],{d,0,m}];
If[OptionValue[diff], Differences[Prepend[h,0]],h]]

s=nrref[macaulayMD[{x^2,y^2},3,{0,0},{x,y}],dTol][[1]]
hilbertFunctionMD[s, 3,2]

29; Multiplicity of isolated point of, square,non-linear system in arbitrary number of variables, including 2.  For non-square systems it will run but answer must be interpreted correctly.

Old version, gives hilbert function, suggestion: use expected depth as m.  An error message could indicate a larger tolerance should be used.

multiplicity0MD[F_,m_,p_,X_,tol_]:=Module[{M,n,l,A,h},
n=Length[X];
M = macaulayMD[F,m,p,X];
{l,A}=nrref[M,tol];
h=hilbertFunctionMD[l,m,n];
Echo[h,"hilbert Function"];
Echo[Length[Select[h,#>0&]]-1,"Depth"];
If[h[[m+1]]>0, Echo[h[[n+1]],"Warning: possibly not isolated, otherwise use higher m"]];
Total[h]]

multiplicity0MD[{x^2,y^2},4,{0,0},{x,y},dTol]

Note : The second term of this local Hilbert function when p is non-isolated gives the local dimension of the solution at p.  For example using the twisted cubic

multiplicity0MD[tBasis3,4,{2,4,8},{x1,x2,x3},dTol]

So this is locally a point on a curve.  On the other hand for the solution set of {x z, y z} in 3 space

multiplicity0MD[{x z,y z}, 3,{0,0,1},{x,y,z},dTol]

multiplicity0MD[{x z,y z}, 3,{1,1,0},{x,y,z},dTol]

These two are different since this algebraic set is the union of a line and a plane.

Faster version in most cases

multiplicityMD[F_,p_,X_,tol_]:=Module[{ttd,svdl,cols,rnk,k,M,h,dh,hk},
ttd=Total[tDegMD[#,X]&/@F];
k=tDegMD[F[[1]],X];
dh=1;
h=0;
While[k≤ ttd && dh>0,
M=macaulayMD[F,k,p,X];
cols=Last[Dimensions[M]];
rnk=Length[SingularValueList[N[M],Tolerance->tol]];
hk=cols-rnk;
dh=hk-h;
h=hk;
k++];
h]

Experimental alternate

multiplicity2MD[F_,p_,X_,tol_]:=Module[{ ttd,M, svdl,cols, rnk,h},
ttd=Total[tDegMD[#,X]&/@F];
M=macaulayMD[F,ttd,p,X];
cols=Last[Dimensions[M]];
rnk=Length[SingularValueList[N[M],Tolerance->tol]];
cols-rnk]

30;  (Changed 5/2020) This function takes a non-linear system and order m then attempts to find small H-Basis in the sense of Macaulay. Will give a correct H-Basis only only if m is large enough.  Unfortunately, a-priori one can’t tell what m is large enough.  This version works well with exact or well conditioned systems, otherwise one must play with the the tolerance.  An unintended feature of this version is that it will often return integer or obvious rational coefficients, even for numerical systems , if possible.  We have eliminated some of the options, instead this version flashes a progress report if it will take a long time  so the user can decide to Abort or not.

Options[hBasisMD]={quiet->False};
hBasisMD[F_,m_,X_,tol_,OptionsPattern[]]:=Module[{n,p, S,A,a,Sa,H,H1,r,s,s1,k,temp},
n=Length[X];
H={};
H1={};
S=sylvesterMD[F,m,X];
{p,A}=arref[S,tol];
If[!OptionValue[quiet],Echo[hilbertFunctionMD[p,m,n],"Initial Hilbert Function "]];
r=Length[p];
H1={A[[1]].mExpsMD[m,X]};
H=H1;
S=sylvesterMD[H,m,X];
s=Length[SingularValueList[S,Tolerance->tol]];
k=2;
While[k≤r,
H1=Append[H,A[[k]].mExpsMD[m,X]];
Sa=sylvesterMD[H1,m,X];
s1=Length[SingularValueList[Sa,Tolerance->tol]];
If[s1>s,s=s1;H=H1];
If[r>30&&Mod[k,10]==0,temp=PrintTemporary["hBasis: At equation ",k, " of ",r];Pause[3];NotebookDelete[temp]];
k++];
S=sylvesterMD[H,m+1,X];
{p,A}=arref[S,tol];
If[!OptionValue[quiet],Echo[hilbertFunctionMD[p,m,n],"Final Hilbert Function "]];
H];

{f1,g1}={x+y,x-z};
f=Expand[f1*z+g1]
g=Expand[f +f1]
hBasisMD[{f,g},2,{x,y,z},dTol]
hBasisMD[{f,g}, 4,{x,y,z},dTol]

The second run with m=4 gives the correct H-basis for this quadratic system but not m<5.

Javanbakht and Sauer (2019) give the following example zero dimensional example known as Weispfenning94 which has H-basis  consisting of 11 polynomials of degree 6.  Since it is zero dimensional we use option  zeroD→True.

FJS={x^2-2 x y+y^2+y^4+x y^2 z+z^2,-2 x y-x^3 y^2+y^4+x y^2 z+x y z^3,-3-2 x^2 y+x y^4+y z^4};

hBasisMD[FJS,6,{x,y,z},dTol]

This does not agree with Javanbakht and Saur but here is a correct run, suppressed. It takes a a very long time to run but you could change cell property to Evaluatable.

FJSB=hBasisMD[FJS,9,{x,y,z},dTol]

Javanbakht and Saur are apparently wrong, see below.

One useful application of hBasisMD is to convert slightly inconsistent systems to consistent systems.  We use an example from 3.2.5 of an inconsistent linear system of 4 equations in 3 variables.

lsys={0.277262174208273` +0.5144436627966619` x+0.10598605713201434` y+0.8045124985078014` z,0.7202433507070195` -0.626115747655119` x+0.27531975154910765` y+0.1158776108984591` z,-0.258819002786987`-0.3361601677996709` x-0.0989359825030034` y+0.9001226231727609` z,0.4685805085964169` -0.849601020102557` x+0.17911927833946` y-0.16287018674812506` z}

NSolve[lsys]

So this system in inconsistent.  But

hsys=hBasisMD[lsys,1,{x,y,z},1.*^-8]

{x,y,z}/.NSolve[hsys]

31 hBasisMDQ : (note hBasisQ has been replaced by hBasisMDQ, but still exists) This checks possible H-Basis H for F by comparing with DegreeLex Groebner basis which is known to be an H-Basis but will probably be larger with large integer coefficients if the system F has integer coefficients.  This works by checking each member of the Groebner basis  of tDegree k to see if it can be expressed as a polynomial combination of H using terms of degree k or less.  If it fails then it will give the degree and position of the first Groebner basis member to fail the test. hBasisQ is numerical, hence a tolerance is needed, the input system can be integer or numeric, in the latter case the Groebner Basis may not compute correctly.  In the form  hBasisQ[F,F,X,tol] it checks F to see if F it is already a H-Basis, in the form hBasisQ[F,H,X,tol]  it checks that H is an H-Basis for F.  For the convenience  of the user the Groebner basis used sorted by degree is saved under the name lastHBGroebner.  One can use hBasisQ then to construct a more efficient H-Basis from the Groebner Basis, see example below.  An option, useF->True tests H against a known H-Basis F.  This is slightly confusing since, in either case,  the first argument gives the known H-basis and it is always  the second  argument is the one to be tested.

Warnings: This function may return false negatives, take a long time or crash, output True should be fairly reliable. This function may not work for some numerical systems.

Options[hBasisMDQ]={useF->False};
hBasisMDQ[F_,H_,X_,tol_,OptionsPattern[]]:=Module[{G,m ,j,degG,degH,selG,SH,SG,r1,r2},
G=If[OptionValue[useF],G=F,G=GroebnerBasis[F,X,MonomialOrder->DegreeLexicographic]];
m=Max[tDegMD[#,X]&/@G];
degG=Sort[DeleteDuplicates[tDegMD[#,X]&/@G]];
G=SortBy[G,tDegMD[#,X]&];
lastHBGroebner=G;
If[MemberQ[degG,0],Echo["F not proper ideal"];Return[False]];
Echo[{Length[G],degG},"{size of Groebner Basis, degrees}"];
degH=Sort[DeleteDuplicates[tDegMD[#,X]&/@H]];
If[degG[[1]]<degH[[1]],Echo[degG[[1]],"No poly in H of degree "];Return[False]]; 
Catch[Do[SH=sylvesterMD[Select[H,tDegMD[#,X]≤k&],k,X];
r1=Length[SingularValueList[N[SH],Tolerance->tol]];
r2=r1;
selG=Select[G,tDegMD[#,X]==k&];
j=Length[selG];
i=0;
While[r1==r2&&i<j,
i++;
SG=sylMD[selG[[i]],k,X];
r2=Length[SingularValueList[N[Join[SH,SG]],Tolerance->tol]]];
If[r1==r2,Continue[],Echo[{k,i},"Problem at poly i degree k "];Throw[Return[False]]],
{k,degG}]];
True]
Options[hBasisQ]={useF->False};
hBasisQ[F_,H_,X_,tol_,OptionsPattern[]]:=hBasisMDQ[F,H,X,tol]

F1={x-z+x z+y z,2 x+y-z+x z+y z};
H1={x+y,x-z};
hBasisMDQ[F1,F1,{x,y,z},dTol]
hBasisMDQ[F1,H1,{x,y,z},dTol]

Here is a difficult quadratic  space curve example in 5 variables called due to Liu (See Javanbakht-Sauer)

Liu={y z-z w-x-u,-x y+y w-z+u,x w-z w-y+u, -x y+x z-w+u};

hBasisMDQ[Liu,Liu,{x,y,z,w,u},dTol]

We append the missing polynomial to the original system and check :

hLiu=Append[Liu,lastHBGroebner[[5]]]

Timing[hBasisMDQ[Liu,hLiu,{x,y,z,w,u},dTol]]

The option useF→True is also useful for the membership question for example one may want to know if

f=-(u^2)-(u^3)-3 (u^2) x-3 u (x^2)-(x^3)-2 u w y+2 u x y-(w^2) (y^2)+2 w x (y^2)-(x^2) (y^2)+2 u z-3 (u^2) w z-6 u w x z-3 w (x^2) z+3 (u^2) y z+2 w y z-2 x y z+6 u x y z+3 (x^2) y z-(z^2)-3 u (w^2) (z^2)-3 (w^2) x (z^2)+6 u w y (z^2)+6 w x y (z^2)-3 u (y^2) (z^2)-3 x (y^2) (z^2)-(w^3) (z^3)+3 (w^2) y (z^3)-3 w (y^2) (z^3)+(y^3) (z^3);

is in the span of hLiu

hBasisMDQ[{f},hLiu,{x,y,z,w,u},1.*^-10,useF->True]

Testing our result for the Javanbakht and Saur  example

hBasisMDQ[FJS,FJS,{x,y,z},1.*^-9]

32 Tangent vector in non-naive case.  Fixed tolerance can be changed inside function. Note: For square systems of breadth one (see Dayton-Li-Zeng) isolated points mimic curves and have unique tangent vectors.  We first show a simple version using just the Jacobian.  Note that if there is no unique tangent vector a warning is displayed and the zero vector is returned for continuity of output.

tangentVectorJMD[F_,p_,X_]:=Module[{J,ns},
If[Norm[F/.Thread[X->p]]>1.*^-7,Echo["Large residue, use tangentVectorMD"];Return[Fail]];
J=D[F,{X}]/.Thread[X->p];
ns=NullSpace[J];
If[Length[ns] ==1,Return[ns[[1]]],Echo[p,"no unique tangent vector at"]];
Table[0,{Length[X]}]]

tangentVectorJMD[{-y+x z,-x-(x^2)+y z,-1-x+(z^2)},{0,0,1},{x,y,z}]

But this version may not work in various cases.   Iit does not distinguish regular points on curves from isolated points.  In particular for square systems an isolated multiple zero of breadth 1 [Dayton-Li-Zeng] may locally mimic a regular curve point, also in 2D it may not recognize a singular point.

 Here is a more complicated version using the Macaulay Matrix, an important benefit is that  It returns also the Hilbert function which should be all 1’s for a regular, non-isolated point on a curve but stabilizes at the singular point multiplicity for singular points.  If there is a zero in the Hilbert function the point is likely isolated, although there may be a tolerance issue, default tolerance 1.×(10^(-7)) may be changed by  option tol→whatever. If all ones then it probably is not,  but this is not proof because a zero could turn up with higher order.  The option order->k will calculate the first k+1 terms of the Hilbert function.  k must be at least 2, setting the order at 2 will make this faster but gives little information about whether the point is isolated.  To leave off Hilbert Function use option hilbertFunction→False.

Options[tangentVectorMD]={tol->1.*^-7,ord->4,hilbertFunction->True};
tangentVectorMD[F_,p_,X_,OptionsPattern[]]:=Module[{M2,n,pv,orth,J,hf},
If[OptionValue[ord]<2, Echo["ord must be at least 2"];Abort[]];
n=Length[X];
{pv,M2}=nrref[macaulayMD[F,OptionValue[ord],p,X],OptionValue[tol]];
If[AnyTrue[Flatten[Take[M2,All,1]],Abs[#]>OptionValue[tol]&],Echo[p,"Not a solution, p = "];Return[]];
hf=hilbertFunctionMD[pv,OptionValue[ord],n];
If[OptionValue[hilbertFunction],Echo[hf,"Hilbert Function"]];
If[hf[[OptionValue[ord]+1]]==0,Echo[ "point may be isolated","Warning"]];
If[hf[[2]]==1,Return[Take[Orthogonalize[Append[M2[[1;;n-1,2;;n+1]],RandomReal[{-1,1},n]]],-1][[1]]],
Echo[p,"No unique tangent vector at "]];
Null];

This  example is a zero dimensional system - all isolated points and not a curve.

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

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

tangentVectorMD[{x z-y, y z-x^2-x,z^2-x^2-1},{0,0,3},{x,y,z}]

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

The next two examples are from the Shen - Yaun curve

SY={-3-7 x-5 (x^2)-(x^3)+7 y+9 x y+3 (x^2) y-5 (y^2)-3 x (y^2)+(y^3),-(x^2)-(x^3)+2 x z+3 (x^2) z-3 x (z^2)+(z^3),-3 y+4 (y^2)-(y^3)-2 y z+3 (y^2) z+6 (z^2)-3 y (z^2)+(z^3)};
rdpt1={7.72614233653476`,14.291303050025972`,5.0724460564306995`};

tangentVectorMD[SY,rdpt1,{x,y,z}]

tangentVectorJMD[SY,{-1,1,0},{x,y,z}]

tangentVectorMD[SY,{-1,1,0},{x,y,z}]

The following is an isolated point which may be difficult to distinguish from a regular point on a curve.

tangentVectorMD[{x+y+z-x^4,x+y+z-y^4,x+y-z-z^4},{0,0,0},{x,y,z}]

tangentVectorMD[{x+y+z-x^4,x+y+z-y^4,x+y-z-z^4},{0,0,0},{x,y,z},ord->7]

35; Subroutines for FLTMD

fVecMD[f_,m_,X_]:=Module[{n,FA,d},
n=Length[X];
{FA,d}=fAssocMD[f,X];
Values[shiftFAMD[FA,0,m]]]

fMatMD[F_,m_,X_]:=Table[fVecMD[F[[i]],m,X],{i,Length[F]}];

gMapMD[T_,m_,X_,Y_]:=
fMatMD[Expand[mExpsMD[m,Y]/.Thread[Y->T]],m,X]

36; Apply a FLT to a numerical point

fltMD[p_,A_]:=TransformationFunction[A][p]

A={{1 ,2,-2,0},{2,-1,4,0},{3,.25,.7,1}};
fltMD[{1,2,3},A]

37 Homogenize affine polynomial with last variable t

homogMD[f_,X_,t_]:=Module[{h,fass,hrl,deg},
fass=Association[CoefficientRules[f, X]];
deg=Max[Total[#]&/@Keys[fass]];
hrl=Table[Append[k,deg-Total[k]]->fass[k],{k,Keys[fass]}];
FromCoefficientRules[hrl,Append[X,t]]];

homogMD[2+3x y+4z^3, {x,y,z},t]

38. Fractional Linear Transformations:   Fractional Linear Transformation applied to curve.

FLT3D  First we have a simple version for  non-singular transforms, it was written for Chapter 1 for naive curves but will actually work for any n, including 2,  so the 3D designation is mostly to distinguish this function from FLTMD.  Works on equations one by one, so returns the same number of equations as input.

FLT3D[F_,A_,X_]:=Module[{B,d,g,h,t,n},
n=Length[X];
If[Dimensions[A]≠{n+1,n+1},Echo[{n+1,n+1},"need A to be of size"];Abort[]];
If[MatrixRank[A]≠ n+1,Echo["A must be invertible"];Abort[]];
B=Inverse[A].Append[X,t];
Reap[Do[
d=tDegMD[f,X];
g=Expand[t^d(f/.Thread[X->X/t])];
h=Expand[g/.Thread[Append[X,t]->B]];
Sow[Chop[h/.{t->1},dTol]],{f,F}]][[2,1]]]

This will take a naive curve to a naive curve, more generally convert any system equation by equation.

The general version from (ℝ^n)⟶(ℝ^s) uses a (s+1)×(n+1) matrix A in FLT form,  m order of Sylvester matrix to use, tol is tolerance.  m may need to be found by trial and error, result may need human interpretation or cleaning.  This function may take a while to run, option timing ->True  will notify user as to progress.  (Note: Change 11/2019: Hilbert Function will not be reported unless option hilbertReport→True is used because reporting will slow down function ). If calculating Hbasis is too slow that can be turned off, hBasis→False, to get many raw equations but be careful to check these for consistency with fltMD.  One can expect better results if the original F is a hBasis.

Options[FLTMD]={timing->False, hilbertReport->False,hBasis->True,quiet->False};
FLTMD[F_,A_,m_, X_, Y_,tol_,OptionsPattern[]]:=Module[{H,XH,YH,T,S,DS,G,TDS,ST,B0,B1,B,n,s,time},
time=TimeUsed[];
n=Length[X];
s=Length[Y];
If[Dimensions[A]≠{s+1,n+1},Print[Style["Dimension Error A",Orange]];Abort[]];
XH=Append[X,#x];
YH=Append[Y,#y];
H=Table[homogMD[f,X,#x],{f,F}];
T=A.XH;
G=gMapMD[T,m,XH,YH];
S=sylvesterMD[H,m,XH];
If[OptionValue[timing], Echo[TimeUsed[]-time,"Start Dual"]]; 
DS=dualMatrix[N[S],tol];
TDS=G.DS;
If[OptionValue[timing],Echo[{Dimensions[TDS],MatrixRank[TDS]},"Dim TDS,rank TDS"]]; 
ST=localDualMatrix[TDS,tol];
If[Length[ST]==0,Print[Style["Fail, try larger m",Orange]];Abort[]];
B0=ST.mExpsMD[m,YH];
If[!OptionValue[hBasis],Return[B0/.{#y->1}]];
If[OptionValue[timing], Echo[TimeUsed[]-time,"Start HBasis"]]; 
B1=hBasisMD[B0,m,YH,tol,quiet->OptionValue[quiet]];
B=B1/.{#y->1};
If[OptionValue[timing], Echo[TimeUsed[]-time,"Total Time"]]; 
B];

G=FLTMD[{x^2+2x y+4z^2-9,x-y+z-1}, fCompProj[3,3],3,{x,y,z},{x,y},1000*dTol]

Check:

cp= criticalPoints3D[{x^2+2x y+4z^2-9,x-y+z-1},{x,y,z}]

G[[1]]/.Thread[{x,y}->fltMD[cp[[2]],fCompProj[3,3]]]

checks!

40;  Rotation about origin taking vector u to vector v, see Wolfram help for RotationMatrix.

(* see Mathematica RotationMatrix[{u,v}] for description, rotates about origin only *)
uvRotationMD[u_,v_]:=Module[{R,R1,n},
n=Length[u];
If[Length[v]≠n, Print[Style["Dimension Error u,v",Orange]];Abort[]];
R1=RotationMatrix[{N[u],N[v]}];
R=Table[0,{n+1},{n+1}];
R[[1;;n,1;;n]]=R1;
R[[n+1,n+1]]=1;
R]

41; FLT that takes infinite Point in 3D to origin of 3-space

ip2z3D[ip_]:=Module[{p,A},p=Take[ip,3];
A={{1,0,0,0},{0,1,0,0},{0,0,0,1},{1,1,1,0}};
If[Norm[Take[p,2]]<1.*^-6,A,A.uvRotationMD[p,{0,0,1}]]]

B=ip2z3D[{1,3,1,0}]

For work with plane curves

ip2z2D[ip_]:=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=line2D[p,{0,0},#x,#y];
B.klRotation2D[ral,#x,#x,#y]]

ip2z2D[{0,1,0}]

Needs (see plane curve book documentation)

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

42; Like fltMD but allows infinite points for input or output.  Important Note: fltMD will take a set of points for input but fltiMD will not.  Use fltiMD[#,A]&/@S to transform a set  S of points.

fltiMD[p_,A_]:=Module[{m,n,u},{m,n}=Dimensions[A]-{1,1};
If[Length[p]==n,u=A.Append[p,1],u=A.p];
If[Abs[u[[m+1]]]<1.*^-6,u[[m+1]]=0;u,Take[Expand[u/u[[m+1]]],m]]]

fltiMD[{1,3,1,0},B]
fltiMD[{0,0,0}, Inverse[B]]

43  Several utility functions are useful.

m2TM[M_] := With[{dim = Dimensions[M]}, Join[Append[M, Table[0, {dim[[2]]}]], Append[Table[{0}, {dim[[1]]}], {1}], 2]]
tM2M[T_] := With[{dim = Dimensions[T]}, Take[T, dim[[1]] - 1, dim[[2]] - 1]]

P = RandomInteger[{-9, 9}, {2, 4}]

T=m2TM[P]

tM2M[T]

factorFLTA   transformation matrix T of dimensions {n + 1, k + 1} will be called an FLT projection if n < k and Rank T = n + 1. Every such transformation can be factored into an affine linear FLT projection and an invertible FLT.The importance is more theoretical than practical as it allows us to concentrate on affine projections. But it is easier to give the code than explain.Note that since the argument to this function is a matrix, not points or equations the suffix MD is not needed.

factorFLT[A_]:=Module[{n,k,m,tab1,tab2,A1,A2,A3,B,B1,B2,B3,P,M},
{n,k}=Dimensions[A]-{1,1};
m=k+1;
tab1=Table[{i,m}->#1[i],{i,n}];
B1=ReplacePart[IdentityMatrix[m],tab1];
A1=A.B1;
B1=ReplacePart[IdentityMatrix[m],(tab1/.Solve[Take[A1,n,-1]==0])[[1]]];
A1=A.B1;
B2=ReplacePart[IdentityMatrix[m],{{m,m}->A1[[n+1,m]]^-1}];
A2=A1.B2;
B3=ReplacePart[IdentityMatrix[m],Table[{m,i}->-A2[[n+1,i]],{i,k}]];
A3=A2.B3;
B=N[B1.B2.B3];
{Chop[A3],Chop[Inverse[B]]}];

Consider the transformation matrix

A={{8,5,5,9,8},{3,4,5,5,8},{2,2,4,9,2}};
MatrixRank[A]
A//MatrixForm

3
(({{8, 5, 5, 9, 8}, {3, 4, 5, 5, 8}, {2, 2, 4, 9, 2}}))

{P,B}=factorFLT[A];
P//MatrixForm
B//MatrixForm

A-P.B

Note that we can almost write down the factorization without calculating but the last column of B does require work.

44  default pseudorandom projections (R^3) to (R^2)

fprd3D={{-0.30519764945947847`,0.9522890290055899`,0.`,0.`},{-0.14191095867181538`,-0.045480825358668514`,0.9888340479238873`,0.`},{0.`,0.`,0.`,1.`}};
prd3D=Take[fprd3D,2,3];

fprd3D
prd3D

45;  FLT for projection  (ℝ^n)→ (ℝ^(n-1))  removing (i^th) component

fCompProj[i_,n_]:=Module[{F},
If[i>n,Abort[]];
F=IdentityMatrix[n+1];
Delete[F,{i}]];

46; Black box function projecting curve F=0, not necessarily naive, from (ℝ^3) to (ℝ^2) where projection pr is given in FLT form, m Sylvester order and rng is the range of a symmetric box about origin

showProjection3D[F_,pr_,m_,X_,{u_,v_},rng_]:=Module[{PRT,AXS,marks},
PRT=FLTMD[F,pr,m,X,{u,v},dTol,quiet->True];
Echo[PRT,"projection Function"];
AXS:=ListLinePlot[{{{0,0},fltMD[{1,0,0},pr]},{{0,0},fltMD[{0,1,0},pr]},{{0,0},fltMD[{0,0,1},pr]}}, PlotStyle->Orange,PlotRange->All];
marks:=ListPlot[{{fltMD[{1.2,0,0},pr]},{fltMD[{0,1.2,0},pr]},{fltMD[{0,0,1.2},pr]}},PlotMarkers->{"x= 1","y=1","z=1"},PlotStyle->Black];
Show[ContourPlot[PRT==0,{u,-rng,rng},{v,-rng,rng}],AXS,marks,Frame->False]]

showProjection3D[{x^2+2x y+4z^2-9,x-y+z-1},fprd3D,3,{x,y,z},{x,y}, 3]

47;  Critical points and forms in 2D

criticalPoints2D[f_,x_,y_]:={x,y}/.NSolve[{f, y D[f,x]-x D[f,y]},{x,y},Reals]

tDeg2D[f_,x_,y_]:=Total[MaximalBy[CoefficientRules[f,{x,y}][[All,1]],Total[#]&][[1]]];
form2D[f_,k_,x_,y_]:=FromCoefficientRules[Select[CoefficientRules[f,{x,y}],Total[#[[1]]]==k&],{x,y}];
maxForm2D[f_,x_,y_]:=form2D[f,tDeg2D[f,x,y],x,y];

f=3+x y-y^3;
criticalPoints2D[f,x,y]
form2D[f,2,x,y]
maxForm2D[f,x,y]

48;  pFiberMD. Mainly subroutine for fFiberMD, gives fiber over a projection parametrically in t. Change 6/20: the affine projection with matrix P is no longer assumed to be orthogonal,  only n×k of rank k-1.

pFiberMD[P_,p_,t_]:=Module[{n,k,P1,ns,q},
{n,k}=Dimensions[P]; 
If[n ≠k-1||MatrixRank[P]≠k-1 ,Echo["not valid Projection"];Abort[]];
P1=Append[P,RandomReal[{-3,3},k]];
ns=NullSpace[P][[1]];
q=Inverse[P1].Append[p,RandomReal[{-3,3}]];
q+t*ns]

P=RandomReal[{-1,1},{2,3}]

pFiberMD[P,{2,3},t]

The following utility programs are for orthogonal projections only.

The projective viewpoint of a projection can be deduced by the direction vector of a fiber, which is independent of the point p. Thus just using the first part of the above we can isolate that point.  In fact all fibers go through this projective point.

projectiveViewPoint[P_]:=Module[{PC,q,m,n},
{m,n}=Dimensions[P];
If[!OrthogonalMatrixQ[P],Print[Style["P not orthogonal",Orange]];Abort[]];
If[m≠ n-1,Print[Style["P not projection matrix",Orange]];Abort[]];
PC=Orthogonalize[N[Append[P,RandomReal[{-1,1},n]]]];
Append[PC[[n]],0]];

projectiveViewPoint[prd3D]

Conversely we may construct an orthogonal projection with given viewpoint by the following.  Note that answer is not unique.

vp2Proj[vp_]:=Module[{P,A,A0,n,v},
If[vp[[-1]]!=0,Echo[v,"not infinite point"];Abort[]];
v=Drop[vp,-1];
n=Length[Drop[v]];
A=Prepend[RandomReal[{-1,1},{n-1,n}],v];
A0=Orthogonalize[A];
Take[A0,1-n]]

P=vp2Proj[{.3,.2,.1,0}];
vp=projectiveViewPoint[P]
Expand[.3/vp[[1]]*vp]

49; P is a projection in non-FLT form from (ℝ^(n ))-> (ℝ^(n-1)), that is a n-1 ×n  matrix, F is curve in (ℝ^n) and p is a point in (ℝ^(n-1)).  This important function attempts to find  points of (ℝ^n) in the fiber over p.  If a finite set it returns all.  Often one needs to use a very loose tolerance.  If more than one point is found it generally implies that p is an artifactual singular point of the projection curve.  If no points are found this function returns {}.  Informational material is Echoed, if no point in fiber is found the number given indicates where in the code this decision is made.  Often this function is used to lift a point path L in (ℝ^(n-1))  to (ℝ^n)   using form Flatten[fFiberMD[F,P,#, X,tol]&/@L,1]

Options[fFiberMD]={complex->False}
fFiberMD[F_,P_,p_,X_,tol_,OptionsPattern[]]:=Module[{Pf,FF,FFs,sol,sol0,sol1,k,n,l,q,u,j,t0},
n=Dimensions[P][[2]];
k=Length[F];
Pf=pFiberMD[P,p,t734];
FF=Chop[Expand[F/.Thread[X->Pf]],tol];
t0=RandomReal[{-1,1}];
FF=SortBy[FF,(#/.{t734->t0})==0&];
If[AllTrue[FF,#==0&],Print["inf many sols at",p];Return[Fail]];
FF=Chop[FF,tol];
If[OptionValue[complex],sol=NSolve[FF[[1]]],sol=NSolve[FF[[1]],t734,Reals]];
If[Length[sol]==0,Echo[p,"(1) no point in fiber at"];Return[{}]];
sol0=t734/.sol;
j=2;
While[j≤k && Length[sol0]>0&&(FF[[j]]/.{t734->t0})≠0,
If[OptionValue[complex],sol=NSolve[FF[[j]]],sol=NSolve[FF[[j]],t734,Reals]];
If[Length[sol]==0,Echo[p,"(2) no point in fiber at"];sol0={};Break[]];
sol1=t734/.sol;
sol0=
Flatten[Reap[Do[If[Norm[q-u]<tol,Sow[q]],{q,sol0},{u,sol1}]][[2]]];
j++];
sol0=DeleteDuplicates[sol0,Norm[#1-#2]<tol&];
If[Length[sol0]==0,Echo[p,"(3) no point in fiber at "]];
If[Length[sol0]>1, Echo[p,"multiple fiber points"]];
Pf/.{t734->#}&/@sol0
]

F={x^2+2x y+4z^2-9,x-y+z-1};
P={{1,0,0},{0,1,0}};
p={-1.1400796343876143`,-2.741829341611946`};
fFiberMD[F,P,p,{x,y,z}, dTol]

L={{-1.1400796343876143`,-2.741829341611946`},{-1.447893202594263`,-2.3770934318145547`},{-1.43809658353943`,-1.7482423269869396`},{-1.2072130512981776`,-1.102629513190312`},{-0.8608698685364934`,-0.5010734718497363`},{-0.25830535044405833`,0.2467487508294853`}};
K=Flatten[fFiberMD[F,P,#, {x,y,z}, dTol]&/@L,1]

50;  Converts a projection, decreasing dimension by 1, to FLT form.

proj2flt[P_]:= Module[{n,T,v},
n=Dimensions[P][[2]];
T=Append[P,Table[0,{n}]];
v=Partition[Append[Table[0,{n-1}],1],1];
Join[T,v,2]]

pr={{1,0,0},{0,1,0}};
proj2flt[pr]

51  A few functions to work with plane curves

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

closestPoint2D[y-x^2, {1,4},x,y]

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

gtVec2D[y-x^2, {2,4},x,y]

52 Path Finder 2D, traces plane curve using  closest point algorithm.  Maximum iterations default is 20 but can be changed with option .

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

pathFinder2D[y-x^2,{0,0},{-2,4}, .1,x,y, maxit->10]

54; line in plane, original formulation.

line2D[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,RandomReal[{-1,1},3]};
({a,b,c}/.(Solve[M.{a,b,c}=={0,0,1}][[1]])).{x,y,1}
];

l=line2D[{1,1},{1,2,0},x,y];
Expand[l/Coefficient[l,y]]

55. PathfinderT2D  - the more traditional and faster path finding algorithm.  s is step size, by default maximum iterations is 30 but can be changed by option maxit.

Options[pathFinderT2D]:={maxit->30};
pathFinderT2D[f_,p_, q_,s_,x_,y_,OptionsPattern[]]:= Module[{L,p0,p1,tv,u,l,k},
p0=p;
L=Reap[Sow[p];
k=0;
While[Norm[q-p0]>2s && k<OptionValue[maxit],
tv=gtVec2D[f,p0,x,y];
u=p0+s*tv;
l=line2D[u,{tv[[2]],-tv[[1]]},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≥OptionValue[ maxit], Print["Warning: iteration limit reached"]];
L];

pathFinderT2D[y-x^2, {0,0},{-2,4},.5,x,y,maxit->10]

56; (DegreeMD discarded, the degree of a space curve should be the degree of a generic plane projection.)

56. Multiplicity of a plane singular point p of curve f=0.  Note that p can be complex, but not infinite.

singPointMult2D[f_, p_,x_,y_,tol_]:=Module[{l},
l=line2D[p,Append[RandomReal[{-1,1},2],0],x,y];
multiplicityMD[{f,l},p,{x,y},tol]
]

f=Expand[-(x-3.1)^2-(x-3.1)^3+(((y+1.4))^2)]
singPointMult2D[f,{3.1,-1.4},x,y,dTol]

57; One step multivariable Gauss Newton, accepts non - square systems, gives least squares solution.

gaussNewton1MD[F_,X0_,X_]:=Module[{Y,Z,J,n,X1,Z1,B},
n=Length[X];
J= D[F,{X}]/.Thread[X->X0];
B=N[F/.Thread[X->X0]];
Y=LeastSquares[J,B];
Z=X0-Y;
(*Echo[{Norm[Y],Norm[F/.Thread[X->Z]]},"{Change, Residue}= "];*)
Z]

gaussNewton1MD[{y^3-x^4,2x+3y}, {.02,.03},{x,y}]
gaussNewton1MD[{y^3-x^4,2x+3y,-2x^3-y^2}, {.02,.03},{x,y}]

k steps of Gauss Newton giving information on convergence.

gaussNewtonMD[F_,X0_,X_,k_]:=Module[{Y,Z,J,n,X1,Z1,B},
Z=N[X0];
Do[
n=Length[X];
J= D[F,{X}]/.Thread[X->Z];
B=N[F/.Thread[X->Z]];
Y=LeastSquares[J,B];
Z=Z-Y;
Echo[Z,"Z"];
Echo[{ Norm[Y],SingularValueList[N[J],n],Norm[F/.Thread[X->Z]]},"{Change, SVL, Residue}= "],{k}];
Z]

gaussNewtonMD[{y^3-x^4,2x+3y,x-y}, {.06,-.03},{x,y},3]

58 Plane affine curve given by full interpolation.

aCurve2D[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=exps[2,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}]
];

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

59  Change of Center subroutine for dual Interpolation

c2zMD[q_,n_]:=Module[{m,Tn,ss,bi,bj,r1,C,s,pow},
pow[a_,m_]:=If[m<=0,1,a^m];
s=Length[q];
Tn=expsMD[s,n];
ss=Length[Tn];
ss=Length[Tn];
C={};
Do[bj=Tn[[j]];
C=Append[C,Table[Product[Binomial[Tn[[i]][[k]],bj[[k]]]*pow[q[[k]],(Tn[[i]][[k]]-bj[[k]])],{k,s}],
{i,ss}]],
{j,ss}];
Transpose[C]]

60; Given a curve and some points this will give a possibly new curve, perhaps an algebraic component, perhaps a different system of defining equations of the original curve.  For our convenience we give this as a black box.  This algorithm works better doing each step separately with separate tolerances.  Thus it is often necessary to try this with different degree m, tolerances and numbers of points to get a useful answer.  Checking is a must.

Options[dualInterpolationMD]:={hBasis->True}
dualInterpolationMD[F_,P_,m_,X_,tol_,OptionsPattern[]]:=Module[{M,DM,DSi,DS,S,G,i,np},
np=Length[P];
DS={{}};
For[i=1,i≤np,i++,
M=macaulayMD[F,m,P[[i]],X];
DM=dualMatrix[M,tol];
DSi=c2zMD[P[[i]],m].DM;
DS=Join[DS,DSi,2]];
S=localDualMatrix[DS,tol];
If[Dimensions[S][[1]]==0,Print["no curve, try larger m"]; Abort[]];
G=S.mExpsMD[m,X];
If[OptionValue[hBasis],Return[Chop[hBasisMD[G,m,X,tol],tol]],Return[G]];
]

(* see example 31 - get Hilbert Basis *)
f=x-z+x z+y z;
g=2 x+y-z+x z+y z;

dualInterpolationMD[{f,g}, {{0,0,0},{-1,1,-1}},5, {x,y,z},1.*^-8]

(* Find a component *)
f=Expand[(x^2+y^2-1)(x-y)];
dualInterpolationMD[{f},{{0,1},{-1,-0}},3,{x,y},1.*^-10]

(*Strip Ramification from Blowup *)
f=y^2-x^3-x^2;
F={f,y-x z};
Φ=Append[#,#[[2]]/#[[1]]]&;
P=Reap[Do[rx=RandomReal[{0,3}];
ry=RandomChoice[{-1,1}]*Sqrt[rx^3+rx^2];
Sow[Φ[{rx,ry}]],{5}]][[2,1]];
{tm,G}=Timing[dualInterpolationMD[F,P,4,{x,y,z},1.*^-7]]

When one gets to larger degrees and variables the  biggest problem with this algorithm is the slowness of the last step, the H-basis algorithm.  The earlier steps still run fairly quickly.  The option hBasis→False will skip the last step.  A very large output will result but individual equations may be approximately true.  With input above

{tm,B}=Timing[dualInterpolationMD[F,P,4,{x,y,z}, 1.*^-7,hBasis->False]];
tm
Length[B]

61; Infinite points of plane curve

infiniteRealPoints2D[f_,x_,y_]:=Module[{fm,rnd,sol},
fm=maxForm2D[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]}]
];

Normalize/@infiniteRealPoints2D[x^3-2x y^2+2x^2 y-y^2, x,y]

62; Numerical attempt to divide multivariable polynomial by another.  Works in all dimensions, but mostly useful in 2D.  Must specify a tolerance.  Used as subroutine for nGCDMD below.  In earlier versions there was a bug, it should be fixed here (1/2021).

nDivideMD[h_,g_,X_,tol_]:=Module[{n,l,m,d1,d2,P,S,f,ex,t},
n=Length[X];
d1=tDegMD[g,X];
d2=tDegMD[h,X];
If [d1>d2,Print["Does Not Divide"];Return[Fail]];
P=PseudoInverse[N[sylMD[g,d2,X]],Tolerance->tol];
S=Chop[sylMD[h,d2,X].P];
ex=expsMD[n,d2-d1];
l=Length[ex];
f=FromCoefficientRules[Table[ex[[i]]->S[[1,i]], {i,l}],X];
t=Expand[f*g-h];
If[NumberQ[t],
If[Abs[t]<d2*tol,Return[f],Print["Does not divide at this tolerance"];Return[Fail]]];
If[Norm[Flatten[sylMD[Expand[f*g-h],d2,X]]]>d2*tol,Print["Does not divide at this tolerance"];Return[Fail]];
f];

g=x^2-3y z;
h=(x^3)+(x^2) y-(x^2) z-3 x y z-3 (y^2) z+3 y (z^2);
nDivideMD[h,g,{x,y,z},dTol]

Numerical GCD by Kaltofen - Dayton algorithm, mostly useful in 2D but works for polynomials of all dimensions, used as subroutine to sqFreeMD below

nGCDMD[ff_,hh_,X_,tol_]:=Module[{nX,S,Sk,U,V,u,v,su,d1,d2,dv,f,h,n,m,k,W,lnull,ex,l},
nX=Length[X];
lnull[M_]:=Dimensions[M][[1]]-Length[SingularValueList[N[M],Tolerance->tol]];
d1=tDegMD[ff,X];
d2=tDegMD[hh,X];
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=sylvesterMD[{f,h},k,X]];
{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+nX,nX]],tol];
ex=expsMD[nX,dv];
l=Length[ex];
u=FromCoefficientRules[Table[ex[[i]]->su[[i]], {i,l}],X];
nDivideMD[h,u,X,tol]];

f=-1.912931182772389` x+1.912931182772389` y+3.605551275463989` (x^2) y-3.605551275463989` x (y^2)-5.738793548317167` z+10.816653826391967` x y z;
h=-1.912931182772389` (x^2)+3.605551275463989` (x^3) y+5.738793548317167` y z-10.816653826391967` x (y^2) z;
g=nGCDMD[f,h,{x,y,z},1.*^-10]
nDivideMD[f,g,{x,y,z}, 1.*^-10]
nDivideMD[h,g,{x,y,z}, 1.*^-10]

63 Square Free reduction of a single multivariable polynomial.  Mostly useful in 2D but works in general. In higher dimensions the concept of square free is replaced by the more complicated notion of  a radical ideal.  That all equations of a curve are square free does not imply that the system is radical.  Dual interpolation may be a way to find a radical system of equations of a curve.  Geometrically in higher dimensions one can plot the radical curve by projecting sufficiently randomly to the plane, using this square free reduction and fiber lifting.

sqFreeMD[f_,X_,tol_]:=Module[{rr,h,g,u,n},
n=Length[X];
rr=RandomReal[{-1,1},n];
h=Grad[f,X].rr;
g=nGCDMD[f,h,X,tol];
If[tDegMD[g,X]==0,Echo["Square Free"]; Return[f],Echo[g," excess factor"]; u=nDivideMD[f,g,X,tol]]];

{f,g}={x z-y,y^2-x^3-x^2 };
sqFreeMD[f,{x,y,z},dTol]
sqFreeMD[g,{x,y,z},dTol]

h=FLTMD[{f,g},fprd3D, 6, {x,y,z},{x,y},100*dTol][[1]]

sh=sqFreeMD[h,{x,y},100*dTol]

64 Functions from plane curves specifically to calculate projective plane singular points.

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

Note FLT2D can be replaced by FLT3D[{f}, A, {x,y}], (see ¶38) needs invertible 3×3  A.

FLT2D[f_,A_,x_,y_]:=Module[{fh,gh, AI,B,d,z},
d=tDeg2D[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]
];

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

complexProjectiveSingularPoints2D[h_,u_,v_,eps_]:=Module[{f,w,A,a,b,sol1,sol2,solz,sing,sing1,sing2},A={RandomReal[{-1,1},3],{0,0,1},RandomReal[{-1,1},3]};
f=FLT2D[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=flti2D[#,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]];

65  random projection to plane for FLT from n-space.

randomOrthogonalFLTplaneProjectionMD[n_]:=Module[{P,FP},
P=RandomReal[{-1,1},{n,n}];
P=Orthogonalize[P];
FP=Append[Take[P,n-1],ConstantArray[0,n]];
Join[FP,Append[ConstantArray[0,{n-1,1}],{1}],2]
]

randomOrthogonalFLTplaneProjectionMD[4]

66. Infinite Real points MD.  Needs subroutine sysMaxForms. May give error message or complex points, use output of sysMaxForms and proceed manually.  A tolerance is last variable, loosening  may delete complex infinite points.

sysMaxForms[F_,X_,tol_]:=Module[{n,d,S,A},
n=Length[X];
d=Max[Table[tDegMD[f,X],{f,F}]];
S=sylvesterMD[F,d,X];
A=arref[S,tol][[2]];
A[[All,Binomial[d+n-1,n]+1;;Binomial[d+n,n]]].Drop[mExpsMD[d,X],Binomial[d+n-1,n]]
]

infiniteRealPointsMD[F_,X_,tol_]:=Module[{maxforms,ssq,sol,ans},
maxforms=sysMaxForms[F,X,tol];
ssq=Total[Table[vv^2,{vv,X}]];
sol=Chop[X/.NSolve[Append[maxforms,ssq-1],X],10*tol];
ans=Append[#,0]&/@Tally[sol,Norm[#1+#2]<.0001&][[All,1]];
Select[ans,Norm[Im[#]]<10*tol&]
]

L=lineMD[{2,3,5,3},{1,0,0,1,0},{x,y,z,w}];
infiniteRealPointsMD[L,{x,y,z,w},dTol]

Some times (eg Bezout' s Theorem) we want complex points too.

infinitePointsMD[F_,X_,tol_]:=Module[{maxforms,ssq,sol,ans},
maxforms=sysMaxForms[F,X,tol];
ssq=Total[Table[vv^2,{vv,X}]];
sol=Chop[X/.NSolve[Append[maxforms,ssq-1],X],10*tol];
If[sol==X,Return[{}]];
Append[#,0]&/@Tally[sol,Norm[#1+#2]<.0001&][[All,1]]]

F2={5-11 (x^2)-3 y-17 x y-17 (y^2)+4 z+2 x z+17 y z-2 (z^2),1+5 x+41 (x^2)-2 y+59 x y+53 (y^2)+4 z-8 x z-59 y z+8 (z^2),1+3 x+9 (x^2)+3 y-5 x y-31 (y^2)+5 z-4 x z+5 y z+4 (z^2)};

infinitePointsMD[F2,{x,y,z},dTol]

67. p2aRawMD -  a multivariable approach to finding equations for polynomial parametric curves; generalization of A.5 in plane curve book.   Takes a parametric curve, F, in n-variables and returns a set of n-variable  equations for hypersurfaces that contain this parametric curve.  This set may be redundant or insufficient as a system describing the curve.  Additional work is necessary to weed out redundant equations (eg. hbasisMD) .  Often d can be the maximum degree of a component of F, but may need to be larger.  Note that this is an improvement of the routines of A.5 for 2 variables.   But this assumes the polynomials are exact or accurately described numerically, otherwise you may want to re-write this with a numerical Nullspace finder.

p2aRawMD[F_,d_,t_,X_]:=Module[{n,TB,cr,ar,SA,NSA,FA},
n=Length[X];
TB=Expand[Table[m/.Thread[X->F],{m,mExpsMD[d,X]}]];
cr=<|CoefficientRules[#,{t}]|>&/@TB;
ar=Append[Flatten[Table[Table[{i,First[k]+1}->cr[[i]][k],{k,Keys[cr[[i]]]}],{i,Length[cr]}],1],{_,_}->0];
SA=SparseArray[ar];
NSA=NullSpace[Transpose[SA]];
If[Length[NSA]<n-1, Echo["Fail, Try higher d"];Abort[]];
FA=NSA.mExpsMD[d,X];
Echo[Table[Expand[FA[[j]]/.Thread[X->F]],{j,Length[FA]}],"Residues "];
FA]

The following example is from   [L.Shen, C.Yaun, Implicitization using Univariate Resultants, J Sys Sci Complex (2010) 23, pp.804 - 814.]

p1=Expand[t^2(t-2)]
p2=Expand[(1-t)^2(t+1)]
p3=Expand[t(t-1)(t-2)]

FR=p2aRawMD[{p1,p2,p3},3,t,{x,y,z}]

The solution set of this system contains the parameterized curve because the residues all vanish on the curve.  Here is a nice implicitization

H=hBasisMD[FR,4,{x,y,z},1.*^-12]

Chop[Expand[H/.Thread[{x,y,z}->{p1,p2,p3}]],1*^-12]

So the implicit curve is contained in the solution of H.  Shen-Yaun claim that the point {-1,1,0}  satisfying their implicitization  is not on the original parameterized curve.   If we take the first two equations of H

Chop[Take[H,2]/.Thread[{x,y,z}->{-1,1,0}], 1*^-12]

This point is numerically on the naive curve defined by the first two equations but since {-1,1,0} is not in H[[3]] it is not on the parameterized curve verifying the Shen-Yaun claim.

The following constants which describe an H-Basis for the parametric curves {t,…,(t^d)} may be useful

tBasis2={(x1^2)-x2};
tBasis3={(x2^2)-x1 x3,x1 x2-x3,(x1^2)-x2};
tBasis4={(x3^2)-x2 x4,x2 x3-x1 x4,(x2^2)-x4,x1 x3-x4,x1 x2-x3,(x1^2)-x2};
tBasis5={(x4^2)-x3 x5,x3 x4-x2 x5,(x3^2)-x1 x5,x2 x4-x1 x5,x2 x3-x5,(x2^2)-x4,x1 x4-x5,x1 x3-x4,x1 x2-x3,(x1^2)-x2};

Note we verify the formulas above for n=2

hBasisMDQ[p2aRawMD[{t,t^2},2,t,{x1,x2}],tBasis2, {x1,x2},dTol]
hBasisMDQ[p2aRawMD[{t,t^2,t^3},5,t,{x1,x2,x3}],tBasis3,{x1,x2,x3},dTol]
hBasisMDQ[p2aRawMD[{t,t^2,t^3,t^4},5,t,{x1,x2,x3,x4}],tBasis4,{x1,x2,x3,x4},dTol]
hBasisMDQ[p2aRawMD[{t,t^2,t^3,t^4,t^5},5,t,{x1,x2,x3,x4,x5}],tBasis5,{x1,x2,x3,x4,x5},dTol]

68. ndReduction
Numerical dimension reduction ndReduction takes an n×d matrix M of rank r   and writes it as a product of an n×r and r×d matrix consisting of the independent rows of M.  Geometrically this writes a linear transformation as a projection onto r components of the range followed by an embedding.

Options[ndReduction] := {tol -> dTol};
ndReduction[M_, OptionsPattern[]] := Module[{S, basis, b, B, NB, NNB, R},
  If[! MatrixQ[M], Echo["Argument must be Matrix"]; Abort[]];
  {n, d} = Dimensions[M];
  b = 1;
  basis = {};
  R = Reap[Do[
     If[Norm[M[[j]]] < OptionValue[tol], Continue[]];
     If[b == 1, Sow[{j, b} -> 1]; AppendTo[basis, j]; b++; Continue[]];
     B = Append[Table[M[[k]], {k, basis}], M[[j]]];
     NB = NullSpace[Transpose[B], Tolerance -> OptionValue[tol]];
     If[Length[NB] == 0, Sow[{j, b} -> 1]; AppendTo[basis, j]; b++; Continue[]];
     NNB = NB[[1]]/NB[[1, -1]];
     Do[Sow[{j, i} -> -NNB[[i]]], {i, Length[NNB] - 1}],
     {j, n}];
    b = Length[basis]];
   {Normal[SparseArray[R[[2, 1]], {n, b}]],M[[basis]]}]

Example:

M={{1,2,3,4},{-2,1,3,5},{-1,3,6,9},{5,0,-2,-3},{4,3,4,6}};
ndM=ndReduction[M]
ndM[[1]].ndM[[2]]

Application to parametric curves:

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

PSM={{-2,1},{-3,-1},{1,1},{0,1}};
{ndP1,ndP2}=ndReduction[PSM]

A1=Join[ndP1,{{1},{2},{3},{1}},2]
A2=Append[Join[ndP2,{{0},{0}},2],{0,0,1}]

Writes P as the composition of two FLT applied to the parabola {t,(t^2)}, the first an invertible linear transformation of (ℝ^2) followed by a projective transformation from  (ℝ^2) to (ℝ^3) showing the resulting ellipse is just a projective transformation from a polynomial curve and must be planar since the second FLT has domain in (ℝ^2.)

P2=fltMD[{t,t^2},A2]
ParametricPlot[P2,{t,-4,4},ImageSize->Tiny,Ticks->None]

Simplify[fltMD[P2,A1]]

ParametricPlot3D[P,{t,-10,10},ImageSize->Tiny,Ticks->None]

69. Bézout’s Theorem: Vectors in nullspace of Sylvester Matrix of order large enough m associated with  point of a system,  aVecMD for affine point, iVecMD for infinite point

aVecMD[p_,m_,X_]:=mExpsMD[m,X]/.Thread[X->p]
iVecMD[p_,m_,X_]:=Module[{lS,lh},
lS=Length[expsMD[Length[X],m]];
lh=Length[hExpsMD[Length[X],m]];
Join[Table[0,{lS-lh}],mhExpsMD[m,X]/.Thread[Append[X,#t]->p]]]

Example : Conjectural  Approximate multivariable Bezout theorem:

F2={5-11 (x^2)-3 y-17 x y-17 (y^2)+4 z+2 x z+17 y z-2 (z^2),1+5 x+41 (x^2)-2 y+59 x y+53 (y^2)+4 z-8 x z-59 y z+8 (z^2),1+3 x+9 (x^2)+3 y-5 x y-31 (y^2)+5 z-4 x z+5 y z+4 (z^2)};

S26=sylvesterMD[F2,6,{x,y,z}];
aSol={x,y,z}/.NSolve[F2]
iSol=infinitePointsMD[F2,{x,y,z},dTol]

nullsp1=Table[aVecMD[p,6,{x,y,z}],{p,aSol}];
nullsp2=Table[iVecMD[p,6,{x,y,z}],{p,iSol}];
nullsp=Join[nullsp1,nullsp2];
MatrixRank[nullsp]
SingularValueList[S26.Transpose[nullsp]]

Thus Transpose[nullsp] is an approximate nullspace for the Sylvester matrix of order 6!

This gives the following  theorem:

Suppose F is an zero dimensional  system of r non-linear real or complex polynomial equations in n≤r  variables X={(Subscript[x, 1]),…,(Subscript[x, n])}  with no multiple solutions  Suppose the equations have degrees  (Subscript[d, 1]),…,(Subscript[d, n]) and m=(Subscript[d, 1])+(Subscript[d, 2])+⋯+(Subscript[d, n])   Let (Subscript[b, a]) be the number of affine solutions and (Subscript[b, i ])be the number of infinite solutions.  Then b=(Subscript[b, a])+(Subscript[b, i]).   Further let k≥m and  for   each affine solution (Subscript[y, j]) let (Subscript[v, i])= aVecMD[(Subscript[y, i]),k,X] and for each infinite solution  (Subscript[z, j]) let (Subscript[w, j]) = iVecMD[(Subscript[z, j]),k,X].  Then (Subscript[v, 1]),…,(Subscript[v, (Subscript[b, a])]),(Subscript[w, 1]),…,(Subscript[w, (Subscript[b, i])]) as column vectors, form an approximate basis for the nullspace of the Sylvester matrix of F of order k.

More generally the dimension of the null space of the Sylvester matrix of high enough order  for a zero dimensional  system counts the number of complex projective solutions counted by multiplicity.

70. nsQSIC3D

This function takes a non-singular QSIC Q with point q and constructs a cubic plane curve h and two rational mappings (not FLT) Ω : (ℝ^3)⟶(ℝ^2) and ℧:(ℝ^2)⟶(ℝ^3) so that Ω maps Q to h and ℧ maps h to Q.  We have Ω[℧[p]]=p for p in Q except for a few points that go to infinity.  These are not inverse functions in the other direction.  This uses a classical trick.

nsQSIC3D[Q_,p_,{x_,y_,z_}]:=Module[{p0,A,F,h,L,M,R,S,Ω,℧},
p0=Normalize[Append[p,1]];
A=Reverse[Orthogonalize[Prepend[RandomReal[{-1,1},{3,4}],p0]]];
F=FLT3D[Q,A,{x,y,z}];
L=formMD[F[[1]],1,{x,y,z}];
M=formMD[F[[2]],1,{x,y,z}];
R=formMD[F[[1]],2,{x,y,z}];
S=formMD[F[[2]],2,{x,y,z}];
h=Expand[L*S-R*M]/.{z->1};
Ω=Take[fltMD[#,A],2]/(fltMD[#,A][[3]])&;
℧=Take[Inverse[A].Join[#,{1,(-R/L)/.Thread[{x,y,z}->Append[#,1]]}],3]/Last[Inverse[A].Join[#,{1,(-R/L)/.Thread[{x,y,z}->Append[#,1]]}]]&;
{h,Ω,℧}]

The preferred syntax is

{hn,Ωn,℧n}=nsQSICMD[Q,q,{x,y,z}];

where n is an integer or other symbol.  Each run of this function gives a different answer so you should distinguish the output.  We use the semicolon because we would rather not see the output.

Caution : the point q is not in the stable domain of Ω in practice one should stay away from this point.

f=x^2+y^2+z^2-16;
g=57-12x+4x^2+y^2-64z+16z^2;
Q={f,g};
q={3.496687222735451`,0.`,1.942467108179341`};
{h1,Ω1,℧1}=nsQSIC3D[Q,q,{x,y,z}];
h1

Note Ω1 is an FLT, ℧1 is more complicated  but each coordinate of ℧[{t,s}] will be a rational function where the numerator and denominator will be a degree 2 polynomial function in t,s. with all denominators are the same.

F=Ω1[{x,y,z}]

p=RandomChoice[criticalPoints2D[h1,x,y]]

℧1[p]
Q/.Thread[{x,y,z}->℧1[p]]

p-Ω1[℧1[p]]

Note that Ω[q] is not in h, the denominator of Ω is zero at this point.

h1/.Thread[{x,y}->Ω1[q]]

Δ=Denominator[F[[1]]]

Δ/.Thread[{x,y,z}->q]

70.1 Weierstrass Normal form for plane cubic (modified from plane curve book).  f is any plane cubic in x,y and ip is an inflection point (see allInflectionPoints2D below).  This returns a cubic in form w=(x^3)+a x+b-(y^2) and an FLT2D transformation matrix taking f to w.

weierstrassNormalForm2D[f_,ip_,x_,y_]:=Module[{B1,B2,B3,B4,B,f1,f2,f3,f4,j,k,c4,c6,cy2},
B1=cTransform2D[f,ip,x,y];
f1=Chop[FLT2D[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=FLT2D[f1,B2,x,y];
B3=homothety2D[CubeRoot[-(Coefficient[f2,x^3]/Coefficient[f2,y^2])],1];
f3 = FLT2D[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=FLT2D[f3,B4,x,y];
B=B4.B3.B2.B1;
{f4,B}];

f=x^3+2x^2 y-y^2x+y^3-2;
ainflptsf=allInflectionPoints2D[f,x,y]
wnf=weierstrassNormalForm2D[f,ainflptsf[[1]],x,y]
ContourPlot[{f==0,wnf[[1]]==0},{x,-4,4},{y,-4,4},ImageSize->Small]

Subroutines for weierstrassNormalForm

cTransform2D[f_,p_,x_,y_]:=Module[{fh,ph, nh, t,cs,A,d},
d=tDegMD[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}];

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

allInflectionPoints2D[f_,x_,y_]:=Module[{H,fh,d,z,l,ips,sol},
fh=Expand[z^tDegMD[f,{x,y}](f/.Thread[{x->x/z,y->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
];

Routine for formatting radical expressions in parameterization of QSIC.

specialExpand[w_,u_,s_,sgn_]:=Module[{w1},
w1=Expand[w/.{s^2->u}];
Collect[w1,s]/.{s->sgn*Sqrt[u]}]

rexp=2+3s+4t-5s^2-6 s t+7t^2
specialExpand[rexp, 3+2t-t^3, s,1]
specialExpand[rexp, 3+2t-t^3, s,-1]

71  Random real regular points in 2D.  Set a bounding rectangle {{a,b},{c,d}} giving x,y ranges.

randomRealRegularPoints2D[f_,{{a_,b_},{c_,d_}},x_,y_,n_]:=Module[{rx,ry,px,py,l,k},
k=1;
Reap[While[k<=n,
rx=RandomReal[{a,b}];
ry=RandomReal[{c,d}];
{px,py}=closestPoint2D[f,{rx,ry},x,y];
If[Abs[Det[D[{f,rx x-ry y},{{x,y}}]]/.Thread[{x,y}->{px,py}]]>.01,Sow[{px,py}];k++]]][[2,1]]
]

f=-(x^2)-(x^3)+(y^2);
rf=randomRealRegularPoints2D[f,{{-2,5},{-5,5}},x,y,5]
f/.Thread[{x,y}->#]&/@rf

72  m2Tm, tM2M switch between matrices and TransformationMatrices  (Now at 43)

Twisted Cubic  - constant

twCubic={x z-y^2,y-x^2,z-x y};

73. Convert Parametric line in n - space to equations.  l is line,  t is parameter, X variables for n-space.

pl2eqMD[l_,t_,X_]:=Module[{n,p1,p2,eq,e1,e2},
n=Length[X];
p1=l/.{t->0};
p2=Append[Table[Coefficient[l[[j]],t],{j,n}],0];
lineMD[p1,p2,X]]

l={3+4t,-2-t,4t};

eq=pl2eqMD[l,t,{x,y,z}]

p3=l/.{t->3}

eq/.Thread[{x,y,z}->p3]

74. Automatic Implicitization for cubics (or less).  P is a list of n+1  polynomials of degree 3 or less in t where X is a list of n-variables. This will represent a rational curve where the n+1 coordinate of P is the common denominator.  This is 1 if it is a polynomial parameterization.

autoImplicitCubicMD[P_,t_,X_,tol_]:=Module[{CF},
CF=Reverse[CoefficientList[#,t,4]]&/@P;
FLTMD[tBasis3,CF,3,{x3,x2,x1},X,tol,quiet->True]]

P={2+t^2,t+t^3,1-t^3,1+t};

This represents rational function

Q=Take[P,3]/P[[4]]

autoImplicitCubicMD[P,t,{x,y,z},dTol]

75. intersection of parametric lines L1, L2 in parameter t and n affine parameters X.  May return { } for no complex projective intersection, a complex affine point of length n or a complex projective point of length n+1 ending in 0.  It also may Fail.  The point of this function is that attempting to find the intersection of two lines given by numerical affine equations using NSolve will often fail since the system needed is overdetermined.  This function uses numerical linear algebra, not NSolve.  (modified 12/21).

nullspace[M_,tol_]:=Take[SingularValueDecomposition[N[M]][[3]],All,-(Dimensions[M][[2]]-matrixrank[M,tol])]

pLineIntersectionMD[L1_,L2_,t_,X_,tol_]:=Module[{n,cr1,cr2,p1,p2,v1,v2,eq1,eq2,S,r,ans},
n=Length[X];
If[Length[L1]≠n,Echo["Line 1 error"];Abort[]];
If[Length[L2]≠n,Echo["Line 2 error"];Abort[]];
p1=Chop[L1/.{t->0}];
v1=Append[Chop[(L1-p1)/.{t->1}],0];
eq1=lineMD[p1,v1,X];
p2=Chop[L2/.{t->0}];
v2=Append[Chop[(L2-p2)/.{t->1}],0];
eq2=lineMD[p2,v2,X];
S=sylvesterMD[Join[eq1,eq2],1,X];
r=matrixrank[S,tol];
If[r<n,Return[{0}]];
If[r>n,Return[{}]];
ans=Flatten[nullspace[S,tol]];
If[Abs[ans[[1]]]<tol,RotateLeft[Chop[ans,tol],1],Take[ans/ans[[1]],-n]]
]

If one only wants to know whether there is an intersection and what type one can use

Length[pLineIntersectionMD[L1, L2, t, X, tol]] which will return 0 for skew lines, 1 for equal lines, 3 for affine point and 4 for infinite point.

parl1={1+2t,3+4t,5+6t,7+8t};
parl2={4.857409341135908` -0.44145714818604453` t,11.00848757377149` +0.02627731817798784` t,16.725434723479104` -0.850047295730326` t,23.0195911974063` +0.0606538619390089` t};

pLineIntersectionMD[parl1,parl2,t,{x,y,z,w},dTol]

The following two lines intersect in the infinite plane in a looser tolerance

l26={1.` -2.` t,0.15469892331131768` +1.9999999999246267` t,-2.3093978465294582`};
l27={1.` -2.` t,-0.42264892323998693`+1.999999999999317` t,-1.1547021535217088`};
pLineIntersectionMD[l26,l27,t,{x,y,z},dTol]
pLineIntersectionMD[l26,l27,t,{x,y,z},1.*^-10]

To show lines are skew one should use a very loose tolerance, eg .003

76. Miscellaneous Transformations

Affine Translation

pqTranslationMD[P_,Q_]:=With[{n=Length[P]},ReplacePart[IdentityMatrix[n+1],Table[{i,n+1}->Q[[i]]-P[[i]],{i,n}]]]

S={{3,2},{4,5},{0,0},{1,2,0},{-3,0,0}}

fltiMD[#,pqTranslationMD[Take[S[[1]],2],{0,0}]]&/@S

Note 2 variable affine translation leaves infinite points fixed.

tangentRealPoints2D[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,{}] ];

homothetyMD[v]:=DiagonalMatrix[Append[v,1]];

Construct Transformation taking set A of d independent vectors in (ℝ^d) to set B in (ℝ^d).  Returns d×d matrix.  IMPORTANT:  the vectors in question are the COLUMNS of A,B.  Function does dimension and rank tests.

getLinTransform[A_,B_]:=Module[{d},
If[Dimensions[A]!=Dimensions[B] || Dimensions[A][[1]]!=Dimensions[A][[2]],Echo["Dimension error, A, B must be square of same size"];Abort[]];
If[matrixrank[A,.0005]!=d,Echo["A singular"];Abort[]];
B.Inverse[A]]

Construct transformation function matrix taking d+1 points in general position in (ℝ^d) to d+1 given other points.  Note points are given as row vectors of A,B.  The function uses the transposes and treats the transformation function as an affine transformation

getTransformMatrix[At_,Bt_]:=Module[{A,B,d,F,G,a,b,AA,BB,L},
A=Transpose[At];
B=Transpose[Bt];
d=Dimensions[A][[1]];
If[!gpTestMD[Transpose[A],d,.0003],Echo["Not general Position"];Abort[]];
a=Take[A,All,-1];
G=Append[Join[IdentityMatrix[d],-a,2],Append[Table[0,{d}],1]];
AA=Transpose[Table[Transpose[A][[i]]-Flatten[a],{i,d}]];
BB=Transpose[Table[Transpose[B][[i]]-Transpose[B][[4]],{i,3}]];
L=BB.Inverse[AA];
F=Append[Join[L,Take[B,All,-1],2],Append[Table[0,{d}],1]];
F.G]

A=RandomInteger[{-5,5},{4,3}]
B=RandomInteger[{-5,5},{4,3}]

M=getTransformMatrix[A,B];
M//MatrixForm

fltMD[{-3,-1,-5},M]
fltMD[{4,-1,0},M]

and so on.

77. General position test.  A set of   points in (ℝ^d)  possibly including infinite points, is in general position if no k+1 points lie in a k-1 dimensional linear set. (Recall that generically k+1 points determine a k dimensional set.)  tol is the tolerance, usually quite loose.   If in general position test returns True.

gpTestMD[S_,d_,tol_]:=Module[{PS,tst,i,k,n},
PS=Table[Switch[Length[s],d,Append[s,1],d+1,s,_,{}],{s,S}];
Do[If[Length[s]==0,Echo["Bad Point"];Abort[]],{s,PS}];
tst=True;
n=Min[d,Length[S]-1];
SS=Subsets[PS,{n+1}];
k=Length[SS];
i=1;
While[i≤k&&tst,
tst=(matrixrank[SS[[i]],tol]==n+1);
i++];
If[tst==False, Echo[SS[[i-1]],"gpTest Fails at "];Return[False],Return[True]];
]

S={{-2,-9,0},{5,-9,4},{6,9,4},{2,2,-7},{9,-9,6`},{1,2,3}};

gpTestMD[S,3,.003]

S={{-2,-9,0},{5,-9,4},{6,9,4},{2,2,-7},{9,-9,6.285714285714281`},{1,2,3}};
gpTestMD[S,3,.003]

78.  Harris' rational normal curve theorem.  In my Mathematica Journal Article Degree vs Dimension…  I identify three numbers for a rational curve, n,d and r, n being the dimension of the  ambient space (R^n), d being the maximal degree of a numerator or the denominator of the parameterization and r the rank of the space of polynomials generated by these numerators and denominators.  If n=d=r  the curve is a rational normal curve, as opposed to the rational normal curve {t,(t^2),…,(t^n)}.  Joe Harris’ Theorem 1.18 says that given d+3 points in general position there is a unique rational normal curve through these 3 points.  the following function gives a parameterization for this curve.  S is a set of d+3  affine or infinite points in general position, k is any non-zero number which will be the parameter value for the last point.   If the points are all affine then the function returns a pair giving first a list of d+3 parameter values and then a rational normal parameterized function passing through these points at the given parameter values.  Note to evaluate at t = ∞ you need to take limit.

S={{9,2,1},{2,-7,8},{-7,-6,4},{-6,5,-2},{-7,-1,3},{5,-9,-8}}

{τ,F}=rncInterpolate[S,2,3]

{F/.{t->τ[[1]]},F/.{t->τ[[2]]},F/.{t->τ[[3]]},F/.{t->0},Limit[F,t->∞],F/.{t->2}}

Like gpTestMD this function will accept points in affine or infinite form.  To get general position at most d can be actually infinite.  Two short routines to help check are included to be used  when option ProjectiveForm->True is used.

Options[rncInterpolate]={ProjectiveForm->False,tolerance->.0003};
rncInterpolate[S_,k_,d_,OptionsPattern[]]:= Module[{A,F0,p,q,p1,q1,τ,PF0,PF1,F},
If[Length[S]≠ d+3,Echo["S must be of length d+3"];Abort[]];
If[gpTestMD[S,d,OptionValue[tolerance]]≠True,Echo["Not General Position"];Abort[],
A=Transpose[Table[Switch[Length[S[[i]]],d,Append[N[S[[i]]],1],d+1,N[S[[i]]],_,Echo["Bad Data"];Abort[]],{i,d+1}]]];
p1=If[Length[S[[d+2]]]==d,N[S[[d+2]]],Take[S[[d+2]]/S[[d+2,d+1]],d]];
p=TransformationFunction[Inverse[A]][p1];
q1=If[Length[S[[d+3]]]==d,N[S[[d+3]]],Take[N[S[[d+3]]]/S[[d+3,d+1]],d]];
q=TransformationFunction[Inverse[A]][q1];
τ=Table[-k(p[[i]]-q[[i]])/q[[i]],{i,d}];
F0=Table[p[[i]]t/(τ[[i]]+t),{i,d}]; (*Special case *)
denom=Expand[Product[(t-τ[[i]]),{i,d}]];
num[i_]:=nDivideMD[Numerator[F0[[i]]]*denom,(t-τ[[i]]),{t},1.*^-10];
PF0=Append[Table[num[i],{i,d}],denom];
PF1=Expand[A.PF0];
F=If[OptionValue[ProjectiveForm],PF1,Table[PF1[[i]]/PF1[[d+1]], {i,d}]];
{Join[τ,{0,∞,k}],F}
]

projectiveLimitMD[F_,rule_]:=Module[{a},
a=Limit[Normalize[F],rule];
If[Abs[a[[-1]]]<.00001,Chop[a],Drop[a/a[[-1]],-1]]]

rncRecoverAffineForm[F_]:=With[{n=Length[F]},Table[F[[i]]/F[[n]],{i,n-1}]]

rncGetProjectiveForm[f_]:=With[{n=Length[f]},Append[Table[Numerator[f[[i]]],{i,n}],Denominator[f[[1]]]]]

Note : in rncGetProjectiveForm all denominators must be exactly the same, including sign .

Examples, first an affine example

S={{9,2,1},{2,-7,8},{-7,-6,4},{-6,5,-2},{-7,-1,3},{5,-9,-8}};
{τ,F}=rncInterpolate[S,1,3]

Checking

{F/.{t->τ[[1]]},F/.{t->τ[[2]]},F/.{t->τ[[3]]},F/.{t->0},Limit[F,t->∞],F/.{t->1}}

Now an example with infinite points and affine points given in projective form

S={{1,5,3,0},{2,1,3},{5,2,-3},{2,5,2},{1,8,6,7},{8,7,9}};
{τ,H}=rncInterpolate[S,1,3,ProjectiveForm->True]

chk=Table[projectiveLimitMD[H,t->τ[[i]]],{i,6}]

In the two cases of disagreement it is because the projective form of the input is not unique.  For S[[1]] this is projective so the check is returning

Normalize[N[S[[1]]]]

In the case of S[[5]] this is an affine point so the check is returning the standard affine form

Drop[N[{1,8,6,7}/7],-1]

Once the result is checked one may want to revert to the affine form of F, one could redo the calculation with the default option or use the following where F is the projective output above.

Ha=rncRecoverAffineForm[H]

79. pfMatMD  Matrix of parametric rational function with parameter t, degree d and common denominator.  Will return condition number.  A different choice of k will change the condition number.  Implicitization with a poorly conditioned matrix (condition number very large) will not work or give bad results

pfMatMD[f_,t_,deg_]:=Module[{n,svd,cd,mat,coeflist},
coeflist:=With[{cl=CoefficientList[#,t]},If[Length[cl]==deg+1,cl,Join[cl,Table[0,{deg+1-Length[cl]}]]]]&;
n=Length[f];
mat=Append[Table[Reverse[coeflist[Numerator[f[[i]]]]],{i,n}],Reverse[coeflist[Denominator[f[[1]]]]]];
svd=SingularValueList[N[mat]];
If[Length[svd]<n+1, Echo["nonsingular matrix"];Abort];
Echo[svd[[1]]/svd[[-1]],"condition number"];
mat]

f={((-1.4957685796193565`+7.866280247681459` t-3.8500000000000014` (t^2))/(2.301182430183626` -5.136194206239249` t+1.` (t^2))),((-5.062601346403977`+8.07490885060259` t-3.5500000000000003` (t^2))/(2.301182430183626` -5.136194206239249` t+1.` (t^2)))}

M=pfMatMD[f,t,2]

fltMD[{t^2,t},M]

This is very well conditioned.  To find implicit equation

eq=FLTMD[tBasis2,M,2,{x2,x1},{x,y},dTol][[1]]

eq/.Thread[{x,y}->#]&/@{{-4.6`,3.85`},{-0.95`,0.2`},{-0.65`,-2.2`},{-3.85`,-3.55`},{4.1`,4.05`}}

This shows that the example f above passes through all 5 points.

79.  Find implicit equation for polynomial  parameterized surface

pol2affNS[F_,d_,md_,T_,X_]:=Module[{n,TB,ar,cr,SA,AK,mon,ncr,nak,NSA,medNSA,FA,SAA},
n=Length[X];
If[Length[F] ≠ n, Echo["Dimension mismatch F,X"];Abort[]];
TB=Expand[Table[mon/.Thread[X->F],{mon,mExpsMD[md,X]}]];
cr=<|CoefficientRules[#,T]|>&/@TB;
ncr=Length[cr];
AK=exps[2,md*d];
nak=Length[AK];
SAA=Reap[For[i=1,i<=ncr,i++,For[j=1,j<=nak,j++,If[KeyExistsQ[cr[[i]],AK[[j]]],Sow[{i,j}->cr[[i]][AK[[j]]]]]]]][[2,1]];
SA=Transpose[SparseArray[SAA]];
NSA=NullSpace[SA];
If[Length[NSA]==0,Return["Fail, try higher md"],Echo[Length[NSA],"Number of equations"]];
medNSA=Median[Abs[Flatten[NSA]]]+1;
N[NSA/medNSA].mExpsMD[md,X]
]

For a rational function F of degree d in  n space use following trick based on the following  idea, there is a connection between stripped  polynomial parametric functions in 4 space and rational parametric functions in 3 space.  Here is an example

Clear[s,t,F];
F={((4+(t^3))/(1+(s^3))),((5+2 (t^2))/(1+(s^3))),((-2+3 t)/(1+(s^3)))};

Stripped means all constants are removed and the denominator becomes the 4th component

strF={t^3,2t^2,3t,s^3};

We can recover F from strF by an FLT

A=Join[IdentityMatrix[4],{{4},{5},{-2},{1}},2];
A//MatrixForm

Then

fltMD[strF,A]

So to find an implicit equation for F

raw=pol2affNS[strF,3,3,{s,t},{x,y,z,w}];

stEqs=hBasisMD[raw,3,{x,y,z,w},1.*^-9]

And finally

eq=FLTMD[stEqs,A,3,{x,y,z,w},{x,y,z},1.*^-9][[1]]

Check:

Chop[Simplify[eq/.Thread[{x,y,z}->F]],dTol]

Closest point to a surface

closestSurfacePointRS[F_,p_,X_]:=FindMinimum[{Norm[X-p],F==0},X];

80  iTransform3D  (Temporarily disabled for repairs)

Sends a given plane to the infinite plane.  Temporarily only works for certain planes.

iTransTemp[a_]:={{1,0,0,0},{0,1,0,0},{0,0,0,1},{1,1,1,a}}

Sends the plane x+y+z+a=0 to infinite plane.

Example :

f=x^2/9+y^2/4+z^2-1;

g=FLT3D[{f},iTransTemp[2.12],{x,y,z}][[1]]

{ContourPlot3D[{f==0,x+y+z+0.32},{x,-4,4},{y,-4,4},{z,-4,4},ContourStyle->Opacity[.6],Mesh->None, ImageSize->Small],ContourPlot3D[g==0, {x,-4,4},{y,-4,4},{z,-4,4},Mesh->None]}

We include the 2 dimensional versions of this idea.

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

ip2z2D[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[[2]]]<1.*^-6,Return[B]];
ral=lineMD[{0,0},p,{x,y}][[1]];
B.klRotation2D[ral,x,x,y]]

gaussPlot2D[f_,r_,pts_,{x_,y_}]:=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]];

Options[inspectInfPoint2D]={output->All};
inspectInfPoint2D[f_,ip_,r_,{x_,y_},OptionsPattern[]]:=Module[{A,g,bp,ap,sol,GP},
A=ip2z2D[ip,x,y];
g=FLT3D[{f},A,{x,y}][[1]];
sol=NSolve[{g,x^2+y^2-r^2},{x,y},Reals];
If[Length[sol]>0,bp={x,y}/.sol;
ap=TransformationFunction[Inverse[A]][#]&/@bp,
bp={}; ap={}];
GP:=Show[gaussPlot2D[g,r,bp,{x,y}], 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}],"ap",Print[GP];Return[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]

81. Hemisphere transform in d variables .  Polar angles are preserved but norm maxes out at  3r.  Helps visualize curves or surfaces in projective 3-space.

ghsMapMD[p_,r_,d_]:=Switch[Length[p],d,N[6*r/Pi*ArcTan[Norm[p]/6/r*Pi]]*Normalize[p],d+1,3*r*N[Normalize[Take[p,d]]],_,Fail];

Show[ContourPlot3D[x^2+y^2+z^2==9,{x,-3,3},{y,-3,3},{z,-3,3},Mesh->None,ContourStyle->Opacity[0.5]],ParametricPlot3D[ghsMapMD[{s,t,s^2+t^2-2},1,3],{s,-30,30},{t,-30,30},PlotRange->All,Mesh->None,PlotStyle->Cyan,MaxRecursion->5],Axes->False,Boxed->False]

82 Normal vectors, normal lines,  tangent planes

For Naive implicit surfaces

tangentPlaneNS[f_,p_,X_]:=(Grad[f,X]/.Thread[{x,y,z}->p]).(X-p)
normalLineNS[f_,p_,X_]:=With[{G=Grad[f,X]/.Thread[{x,y,z}->p]},If[Norm[G]>dTol,lineMD[p,Append[G,0],X],Echo["not regular point"];{}]]

Example:

ts2=-48.`+80.` x-25.` (x^2)+48.` z-80.` x z+25.` (x^2) z-16.` (z^2)+16.` (z^3);
p={2.2928203230275512`,2,-0.5`};

ts2/.Thread[{x,y,z}->p]
Expand[tangentPlaneNS[ts2,p,{x,y,z}]]
normalLineNS[ts2,p,{x,y,z}]

f=-6 (x^2)+4 (x^3)+6 y-16 x y-6 (x^2) y-10 (y^2)+2 (y^3)+12 z+16 x z-2 (z^2);
q={-3,0,-9}
Expand[tangentPlaneNS[f,q,{x,y,z}]]
normalLineNS[f,q,{x,y,z}]

For Rational surfaces

normalVectorRS[F_,st0_,st_]:=Module[{pv,vs,vt},
vs=D[F,st[[1]]]/.Thread[st->st0];
vt=D[F,st[[2]]]/.Thread[st->st0];
Cross[vs,vt]]

tangentPlaneRS[F_,st0_,st_,X_]:=Module[{nv,p},
p=F/.Thread[st->st0];
nv=normalVectorRS[F,st0,st];
N[Expand[nv.(X-p)]]]

F1={s,t,s^2-t^2};
normalVectorRS[F1,{3,2},{s,t}]
tangentPlaneRS[F1,{3,2},{s,t},{x,y,z}]

83 : The following function is given in section 1.3. 2  of the Surface Story to calculate possible naive surfaces, that is polynomial, that contain a given parametric surface.

par2affRS[F_,d_,md_,T_,X_]:=Module[{n,TB,ar,cr,SA,AK,mon,ncr,nak,NSA,medNSA,FA,SAA},
n=Length[X];
If[Length[F] ≠ n, Echo["Dimension mismatch F,X"];Abort[]];
TB=Expand[Table[mon/.Thread[X->F],{mon,mExpsMD[md,X]}]];
cr=<|CoefficientRules[#,T]|>&/@TB;
ncr=Length[cr];
AK=exps[2,md*d];
nak=Length[AK];
SAA=Reap[For[i=1,i<=ncr,i++,For[j=1,j<=nak,j++,If[KeyExistsQ[cr[[i]],AK[[j]]],Sow[{i,j}->cr[[i]][AK[[j]]]]]]]][[2,1]];
SA=Transpose[SparseArray[SAA]];
NSA=NullSpace[SA];
If[Length[NSA]==0,Return["Fail, try higher md"],Echo[Length[NSA],"Number of equations"]];
medNSA=Median[Abs[Flatten[NSA]]]+1;
N[NSA/medNSA].mExpsMD[md,X]
]

For example consider

F={s+t,s^2-2t,s^3-3t s+t};

To find a naive surface containing this of degree at most 3 we do

H=par2affRS[F,3,3,{s,t},{x,y,z}]

Note F is contained in H numerically.

Simplify[H[[1]]/.Thread[{x,y,z}->F]]

In this case the cubic containing F is essentially unique .  If we looked for quartics containing F

H=par2affRS[F,3,4,{s,t},{x,y,z}]

we would get 4 possible surfaces   Incidentally, these are singular surfaces.