A practical guide to understanding and using the Newton Polygon Algorithm
$$ \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)} $$Background:
Singlevalued functions, such as $f(z)=\sin(z)$, can be represented as Taylor series with integer exponents: $$ \sin(z)=z\frac{z^3}{6}+\frac{z^5}{120}+\cdots. $$ A multivalued function, such as $\sqrt{z(z1)(z+1})$, can be represented as a fractional power series: $$ \sqrt{z(z1)(z+1)}=z^{1/2}1/2 z^{5/2}1/8 z^{9/2}+\cdots $$ These fractional power series are called Puiseux series.
Algebraic functions, the subject of this section, are also multivalued.
In previous sections, we computed Puiseux expansions of algebraic functions numerically using Laurent's Theorem. The coefficients of the expansions were necessarily approximate and the powers of $z$ in the expansion were often fractions. The Newton Polygon method ^{1} is an iterative method for computing the exact Puiseux expansion of algebraic functions. In this and the second part of the section, we describe the method and work through several examples.
The following discussion is adapted from Kung and Traub ^{2}. And rather then working with the algorithm immediately, it is more productive to first study briefly, some of the steps used in the method. And although the following may seem quite complicated at first, the interested reader is cautioned not to attempt to understand every detail at this point or get discouraged but rather just briefly review the following and then use it as a reference when implementing the method in the next section. Reviewing this part now, it is hoped the reader will follow and understand much more easily, Part 2 and the actual method in practice.
There are multiple ways of implementing the Newton polygon algorithm. In this and succeeding pages, we review the simplest approach.
Given the expression for $w(z)$ written implicitly as: $$ f(z,w)=a_0(z)+a_1(z)w+a_2(z)w^2+\cdots+a_n(z)w^n=0, $$ our objective is to find a power series expansion for $w(z)$ with center at the origin. Before we review how this is accomplished, it is helpful to first read briefly the following sections and then refer back to this page while we are implementing the method.

Irreducible polynomials:
We will work with irreducible polynomials. Reducible polynomials such as $f(z,w)=g(z,w) h(z,w)$ simply represent two distinct algebraic functions, one for each factor. We can easily check if a particular expression is reducible using the Mathematica function If[! IrreduciblePolynomialQ[theFunction].

Convex Hulls, the Newton Polygon and the lower Newton leg:
Consider a set of points in the xy plane and a border which envelopes the points such as a rubber band stretched around the points. The rubber band will touch the outer points and the remaining points will be inside the border. The border defined by the rubber band is the convex hull of the points. For example the points $(0,9),(1,3),(2,6),(3,5),(4,1),(5,5),(6,10),(7,6),(8,1),(9,8),(10,4)$ have the convex hull shown as the red border in Figure 1 and this is exactly the border that would be created if we were to stretch a rubber band around the points.
Now consider again the general form of an algebraic function: $$ f(z,w)=a_0(z)+a_1(z)w+a_2(z)w^2+\cdots+a_n(z)w^n $$ where $a_i(z)$ is a polynomial in $z$. We now define the order of $a_i(z)$ as the lowest power of $z$ with the exception that if there is no $a_i(z)$, that is no coefficient for $w^i$, then the order of $a_i=\infty$. We write $\textbf{ord}(a_i)$ to represent the order. Then the Newton polygon of $f(z,w)$ is the convex hull of finite points $(i,\textbf{ord}(a_i))$, that is, the points $\textbf{ord}(a_i)=\infty$ are omitted. Let's take a notsosimple example to illustrate how to construct the Newton polygon for a function. Condsider: $$ \begin{aligned} f(z,w)&=(z^{14}+3 z^{15})+(2 z^{10}+5 z^{16}) w^2+(3z^{15}20 z^{16})w^3 \\ &+3z^7 w^4+4 z^8 w^6+3 z w^8+2 w^9+5 w^{10} \end{aligned} $$ and therefore: $$ \begin{aligned} \textbf{ord}(a_0)&=14 \\ \textbf{ord}(a_1)&=\infty \\ \textbf{ord}(a_2)&=10 \\ \textbf{ord}(a_3)&=15 \\ \textbf{ord}(a_4)&=7 \\ \textbf{ord}(a_5)&=\infty \\ \textbf{ord}(a_6)&=8 \\ \textbf{ord}(a_7)&=\infty \\ \textbf{ord}(a_8)&=1 \\ \textbf{ord}(a_9)&=0 \\ \textbf{ord}(a_{10})&=0 \end{aligned} $$ so that the Newton polygon for this function is the set of points $\big\{(0,14),(2,10),(3,15),(4,7),(6,8),(8,1),(9,0),(10,0)\big\}$ shown in Figure 2.
We next define the lower Newton leg as the segments with negative or zero slope such that all remaining point on the convex hull are above and to the right of the lower leg. These segments are color coded in the figure as the red, blue, brown and black legs. Segment $S_1$ has the points $\{(0,14),(2,10)\}$. And $S_2=\{(2,10),(4,6),(8,1)\}$, $S_3=\{(8,1),(9,0)\}$, $S_4=\{(9,0),(10,0)\}$.

The degenerate case:
If we have an expression for $f(z,w)$ such that each polynomial $a_i(z)$ has the form $k+p_i(z)$ where $k$ is a constant, then we obtain for its Newton polygon, a single line segment along the horizontal axis. We call this a degenerate polygon. And in general if all of the convex hull points are along a straight line we have a degenerate case. As explained below, the Newton polygon algorithm is an iterative process where the results of one step are fed into the next step and a degenerate polygon with zero slope is an acceptable lower Newton leg for the first iteration but not succeeding iterations. A degenerate polygon with negative slope would be acceptable for any iteration. An example of this case will be reviewed in the next section.

The characteristic equation for a lower Newton segment:
We now wish to focus on the coefficient on the lowest powers of $z$ for each $a_i$ belonging to each segment. In Figure 2 we have two points along segment 1, three points are on segment 2, and two points for the third and fourth segments. Now form the set of points $(c_i,i)$ where $c_i$ is the coefficient on the power of $z$ for each point of the segment and $i$ is the power of $w$. For segment $S_1$ we would then have $C_1=\{(1,0),(2,2)\}$. That is, the coefficient on the lowest power of $z$ for $a_0$ is 1 and the coefficient of lowest power of $z$ in $a_2$ is 2. Now look at $S_2$. We have the coefficients on the lowest power of $z$ for the points on this segment are $C_2=\{(2,2),(3,4),(3,8)\}$. And likewise, $C_3=\{(3,8),(2,9)\}$, $C_4=\{(2,9),(5,10)\}$. In each of these we have the sets $(c_i,e)$ where $e$ is an exponent. We can then form the characteristic equation for each of these by forming a polynomial: $$ E_i(z)=\sum_{n=1}^N c_i x^e $$ So that for $C_1=\{(1,0),(2,2)\}$ we would have $E_1(x)=1+2x^2$. And for $C_2=\{(2,2),(3,4),(3,8)\}$ we have $E_2(x)=2x^2+3 x^4+3 x^8$ and likewise, $E_3(x)=3x^8+2x^9$ and $E_4(x)=2x^9+5 x^{10}$. These polynomials, $E_i(x)$, are the characteristics equations for the lower Newton legs and during the Newton polygon analysis, we will be concerned with the nonzero roots to these polynomials and their multiplicities. 
Normalizing $f(z,w)$:

Polynomials solutions in $z$
In the case of polynomial solutions $w(z)=a_0+a_1 z^{1/d}+a_2 z^{2/d}+\cdots+a_n z^{n/d}$ where the exponents can be either integer or fractional, the Newton polygon algorithm would, after $n$ iterations, lead to a convex hull without a lower leg signaling that the solution is a polynomial. However, in the implementation described in this paper, we will use regular substitution (described below) before reaching this point. An example of this case will be given in the next section.
One possible issue with regular substitution is that generating a finite number of terms may indicate incorrectly that a solution is a polynomial solution. In the case of the solution $w(z)=a_0+a_1 z^{1/d}+a_2 z^{2/d}+\cdots+a_n z^{n/d}$, regular substitution will after $n$ terms, produce zero for each succeeding coefficient $a_i$. This is expected. However, there are some solutions which have many zero terms only to be followed by nonzero terms. And if we generated an insufficient number of terms for these solutions, the large number of zero terms would suggest the solution is finite. Consider $$ f(z,w)=(1zz^2+z^5)+(8 z^2+8 z^34 z^5)w+(2 z6 z^28 z^32 z^4)w^2+(4 z^2+4 z^3)w^3+(z^2)w^4 $$ Using regular substitution with only 20 terms would produce five nonzero terms and fifteen succeeding zero terms. But the solution is nonfinite; the twentyfirst term is nonzero as are succeeding terms. Had we only used twenty terms, we may have been led into thinking the solution was finite. One way to check this is to backsubstitute the solution into $f(z,w)$. If the results is nonzero, then the solution is not correct. And in fact backsubstituting the first 20 terms of this solution produces a nonzero result.

Polynomials in $w$
For a strict polynomial $f(z,w)=a_1+a_2w+a_3 w^2+\cdots+a_n w^n$ where the $a_i$ are constants, the solutions are simply the roots to the equation $a_1+a_2 w+a_3 w^2+\cdots+a_n w^n=0$. In this case, the algorithm leads, before we begin regular substitution, to a convex hull without a lower leg for each root signaling that we have a polynomial solution (a constant polynomial). The algorithm must identify this case and terminate successfully with each root as the solution. This case will be examined in the test cases.

Simple segments and regular substitution
The Newton Polygon algorithm is an iterative procedure; the results of one step are substituted into the next step and we construct a new Newton Polygon for the current step. If we reach a lower Newton leg consisting of a single segment extending from the vertical axis to the horizontal axis with only the end points and no other points along the segment, the characteristic equation for the segment has only simple roots as it will be of the form $a+x^n=0$. This is called a simple segment and we define all other segments as compound segments. The objective of the Newton polygon algorithm is to iterate all compound segments until we reach simple segments. When we reach a simple segment, we can begin a much simpler substitution in which we use one over the denominator $d$ of the slope of the segment for the remaining exponents $\lambda_i$ and then substitute: $$ w_n=a_n z^{1/d}+a_{n+1} z^{2/d}+a_{n+2} z^{3/d}+\cdots+a_K z^{p+d} $$ for some finite number of terms into the recursive equation described below. 
Newton Polygon Recursive Equation
During the iterative process of constructing the convex hulls, we will be concerned with computing the negative of the slope of a segment, $\lambda$ (different from the $\lambda$ above), the yintercept of the segment, $\beta$, and a coefficient, $c_i$. These are used in the Newton polygon recursive equation which is shown below. However, before we look at the recursive equation, we need to become familiar with the format of the solution we will be working with. This solution will have the form: $$ w(z)=c_1z^{\lambda_1}+c_2 z^{\lambda_1+\lambda_2}+c_3 z^{\lambda_1+\lambda_2+\lambda_3}+\cdots $$ and observe that we can write this solution as $$ \begin{aligned} &=\vdots \\ w_1&=\wvalxxxx\\ &=z^{\lambda_1}(c_1+w_2) \\ &=z^{\lambda_1}(c_1+z^{\lambda_2}(c_2+w_3)) \\ &=z^{\lambda_1}(c_1+z^{\lambda_2}\{c_2+z^{\lambda_3}(c_3+w_4)\})\\ &=z^{\lambda_1}\left(c_1+z^{\lambda_2}\{c_2+z^{\lambda_3}(c_3+z^{\lambda_4}(c_4+w_5))\}\right)\\ &=\vdots \end{aligned} $$ and where the vertical ellipses represents that we can continue the expansions. Understanding this sequence of expressions is a key in understanding the Newton Polygon algorithm. For example, if we have the expression, $$ F(z,w)\equiv F_1(z,w_1)=(z+z^2)+(z1)w_1+w_1^2 $$ and if $w_1=z^{\lambda_1}(c_1+w_2)$, we wish to know $F_1(z,w_1)$. This is easily computed as $$ \begin{aligned} F_1(z,w_1)&=(z+z^2)+(z1)z^{\lambda_1}(c_1+w_2) +(z^{\lambda_1}(c_1+w_2))^2\\ &=F_2(z,w_2). \end{aligned} $$ Notice now this expression has only one coefficient, $c_1$ which we can solve for. Once we compute $c_1$, we then make the substitution $w_2=z^{\lambda_2}(c_2+w_3)$ in $F_2(z,w_2)$, compute another convex hull for $F_2$ if necessary, and then solve for $c_2$. We can continue in this way using $w_3, w_4, \cdots$ as necessary and solve for the coefficients. The recursive equation, which will be explained in the next section in greater detail, is just a compact version of these substitutions. We let:It helps to understand the process if we explicitly write the first few iterations. Let: $$ F(z,w)\equiv F_1(z,w_1) $$ now let $w_1=z^{\lambda_1}(c_1+w_2)$ and $$ F_2(z,w_2)=z^{\beta_1}F_1\left[z,z^{\lambda_1}(c_1+w_2)\right] $$ We then solve for $c_1$. Now let $w_2=z^{\lambda_2}(c_2+w_3)$ and $$ F_3(z,w_3)=z^{\beta_2}F_2\left[z,z^{\lambda_2}(c_2+w_3)\right] $$ and we then solve for $c_2$. So that in general, we write $$ \begin{aligned} w_i&=z^{\lambda_i}(c_i+w_{i+1})\\ &=c_i z^{\lambda_i}+c_{i+1} z^{\lambda_i+\lambda_{i+1}}+\cdots \end{aligned} $$ and then define the recursive equation as $$ F_{i+1}(z,w_{i+1})=z^{\beta_i} F_{i}\left(z,z^{\lambda_i}(c_i+w_{i+1})\right) $$

The Newton polygon algorithm returns all conjugate Pusieux series for each branch
Recall in Section 2 we reviewed how to compute power expansions for algebraic functions using Laurent Theorem. The format of the resulting series was dependent upon the initial conditions of the integrations resulting in $n$ possible versions of the same power series for an $n$cycle branch differing only by the $n$roots of the coefficients. We call the set of $n$ powers series for an $n$cycle branch, the conjugate series for the branch. The Newton polygon algorithm generates all conjugate series for each branch. For example, if we encounter a $4$cycle branch in the analysis, we will likewise encounter one or more polynomials for one or more of the $c_i$ terms leading to $4$ possible cases ultimately giving rise to $4$ conjugate sereis for the branch. See Section 2 for more information about this.

An exact power series obtained by the Newton polygon algorithm is limited by the ability to solve the characteristic polynomial exactly.
The coefficients of the Puiseux series obtained by this method are obtained by finding the roots to a polynomial equation called the "characteristic equation" described in Section 6. Exact solution are only possible for equations of order four or less in the general case. For cases involving degree five or higher, we can resort to numerical methods but with the possibility of encountering numerical errors.

Expanding a function around $p\neq 0 $
In the first three pages of this web site, we considered expansions of the function $w(z)$ around the point $z_0=0$. We can compute the expansions around a point $z_0\neq 0$ by making the simple transformation $F(z,w)=f(z_0+z,w)$ and then an expansion of $F(z,w)$ around zero is an expansion of $f(z,w)$ around $z_0$. Take for example the function studied in Section 2: $f(z,w)=2w^2 z^2+9z^2+2$. If we wanted the power expansion of the function around the point $z=1$, we make the transformation $F(z,w)=2w^2(z+1)^2+9(z+1)^2+2$ and then compute the expansion of $F(z,w)$ around zero.
We will therefore translate the function by $z\to z+z_0$ using this technique anytime we wish an expansion around a point $z_0\neq 0$ so that from now on, we will focus exclusively on expansions about the origin.

Expanding a function around the point $z=\infty$
If we make the transformation $F(z,w)=z^n F(\frac{1}{z},w)$ where $n$ is the highest power of $z$ in the expression, than an expansion of $F(z,w)$ around $z_0=0$ is an expansion around $f(z,w$ around $z_0=\infty$. Take the case again of $f(z,w)=2w^2 z^2+9z^2+2$. Then $z^2f(\frac{1}{z},w)=F(z,w)=2w^2+9+2z^2$.
Normalizing $f(z,w)$ is a crucial step in the method and is a common source of failure with students learning the algorithm for the first time. Recall the power expansions shown in equation $(5)$ of Section $2$. These had a $\frac{a}{z}$ term in their expansions because the function had a pole at the origin. And given $f(z,w)=0$, the function $w(z)$ has a pole at each root of $a_n(z)=0$. In the case of a pole at the origin, $a_n(z)$ has the form $z^kp(z)$ for $k>0$ and $p(z)$ a polynomial. And recall $f(z,w)$ in Section $2$ had a $zw^{12}$ term. In order to use this particular application of Newton Polygon successfully, we must normalize the function and remove the pole at the origin for $f(z,w)$ before we can work with it. And although the following algorithm is tediouslooking, it is the key to working with expansions around poles. We first begin with a definition:
Consider the function $f(z,w)$ of degree $n$ in $w$. We need to find nonnegative integers $\mu,\lambda$ such that: $$ \begin{equation} \begin{aligned} \mu+\textbf{ord}(a_n)&=n\lambda \\ \mu+\textbf{ord}(a_i)&\geq i\lambda, i=1,\cdots, n1 \end{aligned} \end{equation} $$ and we then let $$ F(z,w)=z^{\mu}f\left(z,\frac{w}{z^{\lambda}}\right). $$
Best to work out an example using this algorithm. Consider the function $$ \begin{aligned} f(z,w)&=(1+z)+ (z+z^2) w + z^3 w^2\\ f(z,w)&=a_0(z)+ a_1(z)w + a_2(z)w^2 \end{aligned} $$ for which we obtain $\textbf{ord}(a_0)=0$, $\textbf{ord}(a_1)=1$, and $\textbf{ord}(a_2)=3$. So that we now need to find $\mu$ and $\lambda$ which obeys equation (1) and we can certainly find these values by hand. However, convertToNormalForm computes these values. Refer to Mathematica code 1.a for a description of this routine.
Using the example above, we can code in Mathematica:
convertToNormalForm[myFunction]
And we briefly mention here that we can expand the function around a point other than zero:
References:
1 Algebraic Curves, Walker, R.J.2 All Algebraic Functions can be Computed Fast, H. T. Kung and J. F. Traub
@Paris: Thanks for comments. I corrected lambda expression to mu+ord(a_n)=n lambda as required but your second suggestion is not correct: A pole of w(z) exists whenever a_n(z)=0 and not when a_n(z)!=0.
ReplyDelete