### Section 5: Mathematica Code

Refer to the pages above for more information about these functions. This page is currently under construction. When completed it will contain the complete code for implementing the Newton Polygon algorithm in Mathematica.

$$\newcommand{\bint}{\displaystyle{\int\hspace{-10.4pt}\Large\mathit{8}}} \newcommand{\res}{\displaystyle{\text{Res}}} \newcommand{\wvalx}{\underbrace{z^{\lambda_4}(c_4+w_5)}_{w_4}} \newcommand{wvalxx}{\underbrace{z^{\lambda_3}(c_3+\wvalx)}_{w_3}} \newcommand{wvalxxx}{\underbrace{z^{\lambda_2}\{c_2+\wvalxx\}}_{w_2}} \newcommand{wvalxxxx}{z^{\lambda_1}\big(c_1+\wvalxxx\big)}$$

1. ### {normalizedFunction,lambda,mu}=getNormalizedFunction[function]

This routine normalizes a function. For example:

theFunction=(1+z)+(z+z^2)w+z^3 w^2;
{normalizedFunction,lambda,mu}=getNormalizedFunction[theFunction];


returns
{(z+z^2)+(1+z)w+w^2,2,1}

#### Function code

convertToNormalForm[p_]:=Module[{mymu,temp,flag,mylambda,thepoly,degree
,currentmu,mychecks,mynegs,theclist,theOrders},
theclist=CoefficientList[p,w];
degree=Length[theclist]-1;
mylambda=0;
theOrders=Table[getTheOrder[theclist[[n+1]]],{n,0,degree}];
If[theOrders[[degree+1]]!=0,{mymu[n_,\[Lambda]_]:=n \[Lambda]-theOrders[[degree+1]];
flag=False;
While[flag==False,mylambda++;
currentmu=mymu[degree,mylambda];
If[currentmu>0,{mychecks=Table[currentmu+theOrders[[i]]-(i-1) mylambda,{i,2,degree}];
mynegs=Select[mychecks,#<0&];
If[Length[mynegs]>0,flag=False;,flag=True;];},flag=False;];];
thepoly=Expand[z^currentmu (p/.w->w/(z^mylambda))];},{thepoly=p;}];
{thepoly,mylambda,currentmu}];

getTheOrder[p_]:=Module[{clist,theorder,j},clist=CoefficientList[p,z];
theorder=\[Infinity];
If[Length[clist]!=0,{j=0;
While[clist[[j+1]]==0&&(j+1)<=Length[clist],{j++;}];
theorder=j;},theorder=\[Infinity];];
theorder];


2. ### {newtonPoints, hullPoints, gHull, gPoints, gnewtonLeg, gslopeLines, maxx, maxy} = getPolygonSetup[normalizedFunction, type]; (and associated functions) 

Returns the Newton polygon for a normalized function. If type=1 then segments with zero slope are included with the returned data otherwise only legs with negative slopes are returned. Example:

theFunction=(1+z)+(z+z^2)w+z^3 w^2;
{normalizedFunction,lambda,mu}=getNormalizedFunction[theFunction];
{newtonPoints, hullPoints, gHull,gPoints, gnewtonLeg, gslopeLines, maxx, maxy} =
getPolygonSetup[normalizedFunction, 1];

returns:
newtonPoints={point1,point2,...,pointn} (* points on the lower Newton leg *)
hullPoints={hull1,hull2,...,hullm}      (* all points making up the convex hull *)
gHull=graphics polygon of hull
gPoints=graphics points of data set
gnewtonLeg=graphics of Newton leg
gslopeLines=graphics of slopes of segments
maxx=maximum x for plot
maxy=maximum y for plot


#### Function code

Needs["ComputationalGeometry"];

hexToRGB=RGBColor@@(IntegerDigits[ToExpression@StringReplace[#,"#"->"16^^"],256,3]/255.)&

polyForm[poly_,var_]:=Module[{coeffs=CoefficientRules[poly,var]//Sort},Interpretation[Row[Table[Row[{"(",coeff[[2]],")",var^coeff[[1,1]]/.1->""}],{coeff,coeffs}],"+"],poly]]
myyLine[x_,x1_,y1_,x2_,y2_]:=(y2-y1)/(x2-x1) (x-x1)+y1;

lowerConvexHull[ls_]:=
Module[{
v=ConvexHull[ls],ks,min,max,left,right,pp},ks=ls[[v]];
min=Min[ks[[All,1]]];
pp=Position[ks,{min,_}];
left=Last[pp][[1]];
Append[Take[ks,left-Length[ks]-1],First[ks]]];

getPolygonSetup[fn_,thek_]:=Module[{mypoints,lcp,lcpsave,slope,mylines,newtonleg,sg,convexhull,pgp,maxx,maxy,wdegree},
wdegree=Exponent[fn,w];
mypoints=CoefficientRules[fn ,w]/.Rule[u_,v_]:> {First[u],Exponent[v,z,Min[#&]]};
lcp=lowerConvexHull[mypoints];
lcpsave=lcp;
For[n=1,n<=Length[lcp]-1,n++,
slope=(lcp[[n+1,2]]-lcp[[n,2]])/(lcp[[n+1,1]]-lcp[[n,1]]);
If[thek==1,
{
If[slope>0,
lcpsave=DeleteCases[lcpsave,lcp[[n+1]]];
];
}
,
{
If[slope>=0,
lcpsave=DeleteCases[lcpsave,lcp[[n+1]]];
];
}
];
];

mylines=Table[Plot[myyLine[x,lcpsave[[n,1]],lcpsave[[n,2]],lcpsave[[n+1,1]],lcpsave[[n+1,2]]],{x,0,wdegree},PlotStyle->{Dashed,Blue}],{n,1,Length[lcpsave]-1}];
newtonleg=Graphics[{Red,Thickness[0.007],Line[lcpsave]},Axes->True];
(* for fully-ramified branch in normal form, newton polygon is from first to last point *)

sg=Graphics[{Red,PointSize[0.02],Point @@{mypoints}}];
convexhull=ConvexHull[mypoints];
maxy=Max[(#[[2]]& /@ mypoints)];
maxx=Max[(#[[1]]& /@ mypoints)];
If[Length[mypoints]>2,
(* cannot draw convex hull when less than 3 points *)
pgp=PlanarGraphPlot[mypoints,convexhull,Axes->True,PlotLabel->"Newton Polygon",PlotRange->{{0,15},{0,10}}];
,
pgp=Graphics[Null];
];
{lcpsave,mypoints,pgp,sg,newtonleg,mylines,maxx,maxy}
]


3. ### theSlopeInterceptList=getSlopeIntercepts[newtonPoints]; 

Returns the slope and intercept for each segment on the lower Newton leg. For example:

theFunction=(1+z)+(z+z^2)w+z^3 w^2;
{normalizedFunction,lambda,mu}=getNormalizedFunction[theFunction];
{newtonPoints, hullPoints, gHull,gPoints, gnewtonLeg, gslopeLines, maxx, maxy} =
getPolygonSetup[normalizedFunction, 1];
theSlopeIntercepts=getSlopeIntercepts[newtonPoints];



returns:
{{slope1,intercept1},{slope2,intercept2},...,{slopen,interceptn}}


#### Function code

getSlopeIntercepts[lowerLeg_]:=Module[{lnp,myeqns,myslopes,myinter,mypointlist,slopeinterceptlist},
lnp=lowerLeg;
myeqns=Table[myyLine[x,lnp[[n,1]],lnp[[n,2]],lnp[[n+1,1]],lnp[[n+1,2]]],{n,1,Length[lnp]-1}];
myslopes=Coefficient[myeqns,x];
myinter=myeqns/.x->0;
slopeinterceptlist=DeleteDuplicates[mypointlist];
slopeinterceptlist
];


4. ### coefficientMatrix=getCoefficientMatrix[function]; 

Returns a matrix of exponents and coefficients. The exponents on z may be fractional. Each row of the matrix has an exponent list and coefficient list. The rows are the coefficients of w^n starting with n=0. For example, suppose we are given: \begin{aligned} \text{theFunction}&=(1+2c^3+c^6+c^2 z^{1/3}+c^3 z^5) \\ &+(6 c^2+6 c^5+2 c z^{1/3}+3 c^2 z^5)w \\ &+(6 c+15 c^4+z^{1/3}+3 c z^5) w^2 \\ &+(2+20 c^3+z^5) w^3 \\ &+ 15 c^2 w^4 \\ &+ 6c w^5 \\ &+ w^6 \end{aligned}

then:
coefficientMatrix=getCoefficientMatrix[theFunction]


returns:
{{{0,1/3,5},{1+2 c^3+c^6,c^2,c^3}},
{{0,1/3,5},{6 c^2+6 c^5,2 c,3 c^2}},
{{0,1/3,5},{6 c+15 c^4,1,3 c}},
{{0,5},{2+20 c^3,1}},
{{0},{15 c^2}},
{{0},{6 c}},
{{0},{1}}}


#### Function code

getCoefficientMatrix[thef_]:=Module[{plist,coefficientMatrix,exponents,zlist,array,temp},
plist=CoefficientList[thef,w];
coefficientMatrix={};
For[i=1,i<=Length[plist],i++,
exponents=Exponent[plist[[i]],z,List];
zlist=Coefficient[plist[[i]],z,exponents];
coefficientMatrix=Append[coefficientMatrix,{exponents,zlist}];
];
coefficientMatrix
];


5. ### characteristicEquations=getCharacteristicEquations[function,index]; 

Returns a list of characteristic equations for a Newton polygon. If index=1, returns negative and zero slope segments. Otherwise, returns only negative slope segments. For example, suppose: $$F(z,w)=(z^4)+(2 z^2+z^4)w+(1+z^2+z^3)w^2+(z)w^3+(1/4-z/2)w^4+(-(1/2))w^5$$

then:
{lcpsave,myPoints,pgp,sg,newtonLeg,mylines,maxx,maxy}=getPolygonSetup[normalizedFunction,1]
slopeInterceptList=getSlopeIntercepts[lcpsave]
theftemp=theFunction;
theclist=CoefficientList[theftemp,w];
expList=Exponent[theclist,z,Min];
theEquations=getCharacteristicEquations[minVals,myPoints,lcpsave,1]


returns:
{{1+2c_1+c_1^2=0, c_1^2+1/2 c_1^4-1/2 c_1^5=0}


#### Function code

getCharacteristicEquations[minvals_, mypoints_, lnp_, thek_] :=
Module[{myeqns, myslopes, myinter, mypointlist, slopeinterceptlist, totalsegmentsize, pointlist, myeqnlist, myeqn},
myeqns = Table[myyLine[x, lnp[[n,1]], lnp[[n,2]], lnp[[n + 1,1]], lnp[[n + 1,2]]], {n, 1, Length[lnp] - 1}];
myslopes = Coefficient[myeqns, x]; myinter = myeqns /. x -> 0;
mypointlist = MapThread[{#1, #2} & , {-myslopes, myinter}]; slopeinterceptlist = DeleteDuplicates[mypointlist];
Print["slopeinterceptlist: ", slopeinterceptlist]; totalsegmentsize = Length[slopeinterceptlist]; pointlist = {};
For[i = 1, i <= Length[slopeinterceptlist], i++, pointlist = Append[pointlist, {}];
eqn[x_] := (-slopeinterceptlist[[i,1]])*x + slopeinterceptlist[[i,2]]; For[j = 1, j <= Length[mypoints], j++,
If[eqn[mypoints[[j,1]]] == mypoints[[j,2]], pointlist[[i]] = Append[pointlist[[i]], mypoints[[j]]]; ]; ]; ];
myeqnlist = {}; For[i = 1, i <= Length[pointlist], i++,
myeqn = Sum[minvals[[1 + pointlist[[i,j,1]]]]*x^pointlist[[i,j,1]], {j, 1, Length[pointlist[[i]]]}] == 0;
myeqnlist = Append[myeqnlist, myeqn]; ]; myeqnlist];


6. ### {slopeIntercept,characteristicEquation,rootTally}= getNewSegment[minVals, myPoints, lcpsave, theK]; 

For a Newton polygon, returns for each segment on the lower Newton leg, the slope and intercept of the segment, the associated characteristic equation, and the roots to the characteristic equation along with their multiplicities. For example, suppose: $$F(z,w)=(z^4)+(2 z^2+z^4)w+(1+z^2+z^3)w^2+(z)w^3+(1/4-z/2)w^4+(-(1/2))w^5$$

then:
theFunction=w^2+2z^2 w+z^4+z^2 w^2+z w^3+1/4 w^4+z^4 w+z^3 w^2-1/2 z w^4-1/2 w^5;
theclist=CoefficientList[theFunction,w];
expList=Exponent[theclist,z,Min];
{lcpsave,myPoints,pgp,sg,newtonLeg,mylines,maxx,maxy}=getPolygonSetup[theftemp,theK];
{slopeInterceptList, charEquations, theRootTally} =
getNewSegment[minVals, myPoints, lcpsave, 1];



returns as theSegmentData:
{{{2,4},
1+2 x+x^2==0,
{{-1,2}}},
{{0,0},
x^2+x^4/4-x^5/2==0,
{{0,2},
{1/6 (1+(217-12 Sqrt[327])^(1/3)+(217+12 Sqrt[327])^(1/3)),1},
{1/6-1/12 (1+I Sqrt[3]) (217-12 Sqrt[327])^(1/3)-1/12 (1-I Sqrt[3]) (217+12 Sqrt[327])^(1/3),1},
{1/6-1/12 (1-I Sqrt[3]) (217-12 Sqrt[327])^(1/3)-1/12 (1+I Sqrt[3]) (217+12 Sqrt[327])^(1/3),1}}}}


Where: $\{2,4\}$ is the slope and intercept of the first segment, $1+2x+x^2=0$ is the associated characteristic equation, and $\{\{-1,2\}\}$ is the root tally, that is, $-1$ with multiplicity 2. For the horizontal segment we next have $\{0,0\}$ as the slope and intercept, $x^2+x^4/4-x^5/2=0$ as its characteristic equation, and $\{0,2\}$ as a root of zero with multiplicty two and the three remaining roots with multiplicity one

#### Function code

getNewSegment[minvals_,mypoints_,lnp_,thek_]:=Module[{myeqns,myslopes,myinter,mypointlist,slopeinterceptlist,totalsegmentsize,pointlist,myeqnlist,myeqn,mynetroots,roottally,workingPrecision,rootmatrix,rootaccuracy,rootprecision,ra,eqnlistprecision,workingprecision,eqnlistPrecision,eqnlistAccuracy,minvalaccuracy,myeqnstatus,temp,rootstatuslist,rootstatus,rootprecisionarray,rootsummary},

myeqns=Table[myyLine[x,lnp[[n,1]],lnp[[n,2]],lnp[[n+1,1]],lnp[[n+1,2]]],{n,1,Length[lnp]-1}];
myslopes=Coefficient[myeqns,x];
myinter=myeqns/.x->0;
slopeinterceptlist=DeleteDuplicates[mypointlist];
Print["slopeinterceptlist: ",slopeinterceptlist];
totalsegmentsize=Length[slopeinterceptlist];

(* now have slope-intercept list for each lower segment.  Now check each point to see what points are on those lines *)
pointlist={};
For[i=1,i<=Length[slopeinterceptlist],i++,
pointlist=Append[pointlist,{}];
eqn[x_]:=-slopeinterceptlist[[i,1]] x+slopeinterceptlist[[i,2]];
For[j=1,j<=Length[mypoints],j++,
If[eqn[mypoints[[j,1]]]==mypoints[[j,2]],
pointlist[[i]]=Append[pointlist[[i]],mypoints[[j]]];
];
];
];

myeqnlist={};

For[i=1,i<=Length[pointlist],i++,
myeqn=\!$\*UnderoverscriptBox[$$\[Sum]$, $j = 1$, $Length[pointlist[\([$$i$$]$]]$$]$minvals[$$[$$1 + pointlist[\([$$i, j, 1$$]$]$$$]$]\ x^$(pointlist[$$[$$i, j, 1$$]$])$$\)\);
If[Resultant[myeqn,D[myeqn,x],x]==0,
Print[Style["Multiple roots detected",Blue]];
];

myeqnlist=Append[myeqnlist,myeqn];
];
myeqnstatus={};

Print["Solving for roots"];
mynetroots=x/.Map[Solve[#1==0,x]&,myeqnlist];
Print["Roots solved"];
If[Length[Cases[mynetroots,_Root,\[Infinity]]]>0,
{
Print[Style["Cannot find exact roots to characteristic equation",Red,Bold]];
Print["Char equ: ",myeqnlist];
Print["Roots: ",N[mynetroots]];
Abort[];
}
];

roottally=DeleteCases[#,{0,_}]& /@ (Map[Tally[#1]&,mynetroots]);
{slopeinterceptlist,myeqnlist,roottally}

];


2. ### doPuiseux[thef_, theK_, totalIter_, theLambda_, dFlag_,rootFlag_]

Computes the fractional power series for the function $w(z)$ given implicitly as $$f(z,w)=a_0(z)+a_1(z)w+a_2(z)w^2+\cdots+a_n(z)w^n=0$$ Please refer to the section on doPuiseux for a description of this function.

3. ### getDigits[{vals}_, accuracy_]

Returns an integer representation of a floating-point number with accuracy number of digits.

getDigits[vals_, n_] :=
Module[{net, theDigits, i, theChoppedVals, theComponents, myList,
listAccuracy},
If[Length[vals] == 0,
myList = {vals};
,
myList = vals;
];
listAccuracy = Accuracy[myList];

myList = SetAccuracy[myList, listAccuracy - 2];
net = Length[myList];

theComponents = {Re[#], Im[#]} & /@ myList;
theChoppedVals = Chop[theComponents, 10^(-n)];
theDigits = RealDigits[theChoppedVals, 10, n];
For[i = 1, i <= net, i++,
If[theChoppedVals[[i, 1]] < 0,
theDigits[[i, 1]] = Prepend[theDigits[[i, 1]], -1];
,
theDigits[[i, 1]] = Prepend[theDigits[[i, 1]], 1];
];
If[theChoppedVals[[i, 2]] < 0,
theDigits[[i, 2]] = Prepend[theDigits[[i, 2]], -1];
,
theDigits[[i, 2]] = Prepend[theDigits[[i, 2]], 1];
];
];
theDigits
];
`