On the Dynamics of Laguerre’s Iteration Method for Finding the nth Roots of Unity

Previous analyses of Laguerre’s iteration method have provided results on the behavior of this popular method when applied to the polynomials p n (z) = z n − 1, n ∈ N. In this paper, we summarize known analytical results and provide new results. In particular, we study symmetry properties of the Laguerre iteration function and clarify the dynamics of the method. We show analytically and demonstrate computationally that for each n ≥ 5 the basin of attraction to the roots is a subset of an annulus that contains the unit circle and whose Lebesgue measure shrinks to zero as n → ∞. We obtain a good estimate of the size of the bounding annulus. We show that the boundary of the basin of convergence exhibits fractal nature and quasi self-similarity. We also discuss the connectedness of the basin for large values of n. We also numerically find some short finite cycles on the boundary of the basin of convergence for n = 5, ..., 8. Finally, we demonstrate that when using the floating point arithmetic and the general formulation of the method, convergence occurs even from starting values outside of the basin of convergence due to the loss of significance during the evaluation of the iteration function.


Introduction
Laguerre's iteration method (also referred to as Laguerre's method in this paper) for approximating roots of polynomials [1] is one of the least understood methods of numerical analysis. It exhibits cubic convergence to simple roots of (complex) polynomials and linear convergence to multiple roots, thus outperforming the well-known Newton's method that exhibits quadratic convergence to simple roots [2] or even the widely used and globally convergent Jenkins-Traub method, which has the order of convergence of (at least) 1+ , where = (1 + √ 5)/2 is the golden ratio [3,4]. Although Laguerre's method is implemented in Numerical Recipes [5], it is often overlooked in designing professional software. This is, perhaps, due to the lack of complete analytical understanding of the method. However, some of the known results make it an excellent candidate in many situations. For example, it is known that the method exhibits global convergence (convergence from any initial guess) for real polynomials with real roots [3,6]. It also allows for automatic switching to the complex domain if there are no real roots; this is due to the appearance of a square root in the definition of the method (see (2) in the next section). In general, although convergence is not guaranteed for all complex starting values, the method seems to perform very well in many cases [2].
It is the goal of this paper to provide additional insights into the performance of Laguerre's method when applied to simple polynomials of the form − 1. We primarily follow the work of Ray [7] and Curry and Fiedler [8] but provide additional details and clear proofs of all results. In addition, we provide computational results that demonstrate the poor performance of the method when is large and exact arithmetic is used. This is due to the fact that the basin of convergence to the roots is contained in an annulus that shrinks towards the unit circle 1 as increases. Points in the complement of the annulus converge to the two-cycle {0, ∞}. We also show that the boundary of the basin of convergence has fractal characteristics and becomes quite interesting for large . Finally, we present some examples to demonstrate that in the floating-point arithmetic the method in its general formulation (2) eventually converges from seemingly any initial complex value due to the loss of significance. This thus ironically contributes to the practicality of the method in this case, and it remains to be seen whether this is also the case for general polynomials.
We organize the paper similarly as in [8]. In Section 2, we introduce the method, briefly summarize known results, and provide a simpler expression for the method when applied to the polynomials ( ) = − 1. In Section 3, we formulate and prove three propositions regarding the symmetry properties of the iteration function for that simplify the analysis in the following sections. Sets in the complex plane that will play a significant role in the study of the dynamics are defined in Section 4, and their algebraic characterization is provided in the same section. The dynamics of the method on the unit circle and in the neighborhood of the two-cycle {0, ∞} is studied in Sections 5.1 and 5.2, respectively. In Section 5.3 we provide proofs of convergence to the roots of unity when the initial guess is in a relevant annulus containing the unit circle. The boundary of the basin of convergence is contained in two annuli shown as the "gray areas" in Figure 1, and some relevant numerical results pertinent to the boundary are shown in Section 6. We conclude with Section 7, in which we summarize some of the open questions, and demonstrate the "convergence" of the method even from the basins of attraction of the two-cycle {0, ∞}.

Laguerre's Iteration Method
In this section, we provide the basic details of Laguerre's method, mention known results, and apply the method to the polynomials ( ) = − 1 with ≥ 2 and ∈ C. We will denote by 1/2 the set of the two solutions { , − } ⊂ C such that 2 = (unless, of course, = 0, in which case = 0). We will use the notation √ for the principal square root of ; that is, if = with > 0 and − < ≤ , then √ = √ /2 is the principal value for the square root with argument in the interval (− /2, /2]. Laguerre's iteration method for complex polynomials ( ) of degree ≥ 2 is defined as [1,3] where ( ) denotes the Laguerre iteration function [9] given by where and where the signs are chosen so as to maximize the modulus of the denominators.

Known Results.
It is known that Laguerre's method exhibits cubic convergence to a simple root and linear convergence to a multiple root [2,3]. It is also known that for a real polynomial with real roots, the method converges to a root from any initial guess 0 ∈ R [3]. A particular feature of interest is that even if the initial guess is a real number, convergence to a complex root can occur due to the square root in the denominator of (2). In many cases, the method seems to converge to a root from any initial guess in the complex plane, although this is not the case in general [7]. For example, consider the polynomial ( ) = − 1 with ≥ 3, for which both the first and the second derivative vanish at = 0, and (0) is undefined (in what follows, we will consider the extended complex planeĈ = C∪{∞} so that (0) = ∞ and (∞) = 0, and {0, ∞} forms a two-cycle of the method). It is also known [10] that if ( ) is a polynomial of degree and ∈ C, then there exists a root * of such that | − * | ≤ √ | − ( )|.
The Laguerre iteration function (2) is sometimes claimed to be invariant under Möbius transformations [2,11], although the correct, weaker statement is given and proved in [7]. For the classes of quadratics and cubics of the form ( ) = 2 + and ( ) = ( −1) ( 2 + + ), respectively, with , ∈ C, the Laguerre iteration function (2) can be shown to not have any free critical points [11], and a generalization to all complex quadratics and cubics is suggested based on the invariance under Möbius transformations.
where again the sign is chosen to maximize the modulus of the denominator. Using the principal square root of , which will result in an expression with a nonnegative real part, we can rewrite ( ) as ( ) given by We note that the roots of are exactly the fixed points of , and it is easy to check that the derivative of vanishes at the roots, so they are attracting fixed points, and each has an open neighborhood contained in its immediate basin of attraction.
In the cases with ≥ 3, the method has a two-cycle consisting of 0 and ∞ in the extended complex planeĈ. Other than this two-cycle, the method is globally convergent for = 3, 4 [7,8].
The behavior of Laguerre's method is quite different in the cases with ≥ 5 and is the subject of our interest. In the following sections, we closely follow the analysis of Curry and Fiedler [8], which in turn is based on the work of Ray [7]. We provide additional insights and graphical illustrations for some of the results. While the main focus will be on the cases with ≥ 5, if a result applies more generally, we will state so.
For future reference, we note that for > 0, we have [8] ( )

Symmetry Properties of the Laguerre Iteration Function
When applied to the polynomials ( ) = − 1, ≥ 2, the Laguerre iteration function (5) has several symmetry properties that simplify the analysis of the method in the extended complex planeĈ. In particular, the function possesses anfold rotational symmetry around the origin, a symmetry with respect to the real axis, and also an inversion symmetry with respect to the unit circle. In Propositions 1-3 we provide the precise statements. We will use the slightly imprecise notation of [8] and denote by 0 and 1 any of the following angles for = 0, . . . , − 1: In addition, we define the rays Θ 0 = { 0 ∈ C : > 0} , Θ 1 = { 1 ∈ C : > 0} . (8) Note that the roots of lie on the rays Θ 0 , while the rays Θ 1 divide the complex plane into congruent sectors bisected by the rays Θ 0 .
The following proposition implies that it suffices to study the behavior of in the sector { ∈ C : − / < arg ≤ / }, that is, between two consecutive rays Θ 1 , and the rest follows by rotational symmetry. This is a special case of the invariance of the method with respect to certain Möbius transformations [2,7]. Proof. Let denote the rotation by an angle = 2 / about the origin; that is, ( ) = . Since ( ( )) = , substituting into (5), we get for any ∈Ĉ and the result follows.
The following proposition implies that the behavior of in the sector { ∈ C : − / < arg < / } is symmetric with respect to the real axis. In the case when arg ( ) = 1 , Proposition 1 applies.
Finally, the Laguerre iteration function (5) also exhibits inversion symmetry with respect to the unit circle as shown in the next proposition. Proof. Consider first ∈ C with arg ̸ = 1 . Using the same conjugation properties as in the proof of Proposition 2, we have 4

International Journal of Computational Mathematics
We can then check by direct substitution that the same result holds also when arg ( ) = 1 , and the conclusion of the proposition follows.
In summary, it is enough to study the behavior of the method in the set { : 0 ≤ ≤ 1, 0 ≤ ≤ / }, since the rest follows by the symmetries above.

Significant Sets in the Complex Plane
Consider from now on the polynomial ( ) = − 1 with ≥ 5 and the corresponding Laguerre iteration function given by (5). Following [8], we start by defining several sets in the extended complex planeĈ relevant to the study of the dynamics of Laguerre's method. We will provide relevant results, some of which are proved in [8].
It is stated in [8] that the sets in (11) below "contain all the dynamics" of (5). This is not quite true, as will be shown in Section 6, although these sets are of significance in the analysis. They divide C \ {0} into disjoint subsets and are defined as Note that, due to the use of the principal square root in (5), is continuous everywhere in C \ {0} except across the rays Θ 1 ; however, | | is continuous across Θ 1 , so the notations and are justified, since and are the boundaries of and , respectively, in C\{0}. For future reference we note that is "counter-clockwise" continuous across the rays . This is not the case in the "clockwise" direction.
As in [7,8], we next focus on the algebraic characterization of and . This will lead to a definition and study of a "characteristic function" (14) below that will allow us to determine the shapes of the boundary curves as shown in Figure 1. The figure provides an illustration of the sets in (8) and (11) for = 16; the cases with other values of are similar. In the figure, the three thick solid curves correspond to the solutions of (13) below and are, in the order of increasing distance from the origin, , 1 , and . They divide C \ {0} into four regions: , 0 , 1 , and (innermost to outermost).
The dots indicate the positions of the roots of . The thick dashed lines are the rays Θ 1 , while the thin dotted lines are the rays Θ 0 . The thin dashed circles have radii 0 < 0 < 1/ 0 < 1/ 0 as defined in (18) in Remark 5 below. Writing = , 0 < < ∞, we note that both and are characterized by the same equation: Using (6), (12) is equivalent to [7,8] ( , ) = 0, where the "characteristic function" is defined as We summarize relevant results (some stated in [8]) in the following theorem.
(2) Each of the sets in (11) corresponds to a particular sign of : (3) The boundaries and correspond to polar curves of the form = ( ) and = ( ). The function ( ) is maximized at any = 0 and minimized at any = 1 , while the function ( ) is minimized at any = 0 and maximized at any = 1 . In addition, both ( ) and ( ) are monotonic between any two consecutive angles 0 and 1 (see Figure 1).
Proof. (1) Let ≥ 5, ∈ R, and define ( ) = ( , ). Note that is a differentiable function and It is easy to show that (1) ≤ −2 2 . This implies that has at least three positive zeroes. From (12) and Proposition 3 it International Journal of Computational Mathematics 5 follows that, other than 1, the zeroes of come in reciprocal pairs, so the actual number of zeroes is an odd number greater than or equal to 3. As in [7,8], we will invoke Descartes' rule of signs. When is even, is a polynomial, so the rule can be applied directly. When is odd, we can apply it to ( ) = ( 2 ), which is a polynomial that also satisfies (16) with ( ) replaced by ( ). Hence, we focus on with the understanding that is handled exactly the same way. Note that can be expanded to contain at most 6 terms with different powers of ; hence, there are at most 5 sign changes, and has at most 5 positive zeroes. From (16) it now follows that if had 5 zeroes, two of them would have to have multiplicity greater than 1 and would have to have at least 6 positive zeroes (4 between the zeroes of and at least 2 more from the multiple roots of ). This is, however, impossible, since is another polynomial with at most 5 terms of different powers of ; hence Descartes' rule of signs implies that has at most 4 positive zeroes. Consequently, has exactly 3 simple positive zeroes as stated in the theorem.
(2) This part follows from the definition of the relevant sets in (11) and from replacing the equality in (12) by inequalities, which results in inequalities in (13) [8].
(3) Since for every ∈ R there are unique values of and , we can think of and as polar curves. To prove all of the remaining statements in this part, it is enough to consider ( ) for 0 ≤ ≤ / , since the rest follows by the symmetries discussed earlier. Note that the cosine term in (14) is the largest for = 0, so on the circle = (0) > 1, as a function of , ( (0), ) is the largest (and equal to 0) exactly when = 0 . Thus, is negative on the circle for every ̸ = 0 , and it follows that ( ) ≥ ( 0 ) for all ∈ R. Similarly, the cosine term is the smallest when = / , and by a similar argument we get ( ) ≤ ( 1 ) for all ∈ R. Finally, to prove the last assertion, we implicitly differentiate (13) with respect to and observe that / = −( / )/( / ) vanishes only when = 0 and does not exist only when = 1 , since the numerator contains a factor sin ( /2) and the denominator is positive on ( ) as the proof of part (1) implies. This concludes the proof of the theorem.

Dynamics on the Unit Circle.
In this section, we study the dynamics on the unit circle, 1 . As a consequence of Theorem 4, part (1), we have that the unit circle, 1 , is invariant under the Laguerre iteration function (see also [8]). This follows from (13) and (12) with = 1. However, more can be said about the behavior of on 1 . Proof. Due to the symmetry properties of (Propositions 1 and 2), it is enough to assume = with 0 < ≤ / and show that ( ) =̃with 0 <̃< . Since the sequence of arguments generated by the method will then be a decreasing sequence bounded below by 0, it will converge, and the corresponding sequence of points on 1 will converge to a fixed point of , that is, to a root of by the continuity of away from Θ 1 . The limiting root will have to be 1 and both assertions of the proposition will follow.

Remark 8. The regions
0 and 0 defined in Proposition 7 can be seen in Figure 1: 0 is the open ball not containing 0 bounded by the smaller gray annulus, while 0 is the region on the outside of the larger gray annulus.
Although the basin of attraction of the two-cycle {0, ∞} is significantly larger than 0 ∪ 0 (see Section 6), we can immediately extend it in the following sense. Proof. From the definitions of and in (11) it is clear that any point in {| | = 0 } ∩ or {| | = 1/ 0 } ∩ gets mapped into 0 ∪ 0 , and the claim follows.
We believe that the points 0 1 and (1/ 0 ) 1 also converge to the {0, ∞} two-cycle. Since these points belong to the set ∪ , their images under the Laguerre iteration map (5) lie on the circles with radii 1/ 0 and 0 , respectively, so it suffices to show that their arguments are different from 1 . We have not been able to find a simple proof for this statement.

Dynamics of the Basins of Convergence.
In Proposition 10 we state a result [8] that shows that the open annulus bounded by the gray annuli in Figure 1 belongs to the basin of attraction of the roots of ( ) = −1. Again, it turns out that the basin is actually larger than the annulus (see Section 6). We provide an elementary proof of the final statement of Proposition 10, since in the proof in [8] a reference is made to [6], which does not appear relevant to the proof.
We again have an extension of the above proposition, arguing as in the proof of Corollary 9.  However, the remaining points on the circles {| | = 0 } and {| | = 1/ 0 } form nontrivial, finite two-cycles [7,8]. (18) is a two-cycle for the Laguerre iteration function (5).

Proposition 12. For every
Proof. By Proposition 1, we can assume 0 = 0. From (5) we have that maps real, positive numbers to real, positive numbers, so For completeness, we state the following result that completes the dynamics on Θ 0 . Finally, the following result clearly demonstrates that Laguerre's method is not suitable for finding roots of unity for large-degree polynomials [7].

The Basins of Convergence and Their Boundaries
In this section we present primarily computational results that address the structure of the basins of attraction of Laguerre's method (1) applied to ( ) = −1 with ≥ 5. All numerical results have been generated using the free software package GNU Octave [12,13]. These results show that the basins of attraction do not coincide with the sets defined in (11), and they raise additional questions that we summarize in the next section. We mentioned earlier that for = 2 it takes one iteration to get to a root of 2 from any initial guess. It is also known that for = 3 and 4 the method is globally convergent to a root of (except for when the initial guess 0 = 0) [7,8,11]. However, from the above analysis it follows that for ≥ 5 this is no longer true; more specifically, the basin of attraction for each ≥ 5 is contained in the annulus { 0 ≤ | | ≤ 1/ 0 } with 0 given in (18). Since 0 is not easily computable, we can use the upper bound 1/ 0 < ( − 1) 2/( −4) (see Theorem 4). In Figures 2 and 3 and show that the upper bound is a good estimate of 1/ 0 . Also plotted are the curves and and we see that the basins of attraction do not coincide with any of the sets of significance defined in (11). Also note how, in accordance with Theorem 4, the basin of convergence shrinks as increases. This is demonstrated in Figure 4, where the basins for = 5, 8, 12, and 16 are plotted on the same scale.
The boundary of the basin of convergence appears fractal (see, e.g., [14] for more on fractals). We demonstrate this in Note that the basin of convergence is not either connected or simply connected. In particular, it appears that the basin of attraction of the roots consists of infinitely many (quasi-) self-similar disconnected sets (zooms on the left) and infinitely many (quasi-) self-similar "holes" corresponding to basins of attraction of the two-cycle {0, ∞} (zooms on the right).   The boundaries displayed in Figure 5 appear self-similar, but they are only quasi-self-similar (see, e.g., [14,15] for more on quasi-self-similarity). We demonstrate this observation in Figure 6, where we can see slight changes of shape as we zoom in and also as we more carefully examine the shapes within each figure.
In addition, it came to us as quite a surprise; it seems that the basins of convergence as shown in color in Figures 2 and  3 are not, in general, (disregarding the "hole" in the middle) simply connected or even connected! In Figures 7 and 8 we present results with = 128 and = 1024, respectively, and several consecutive zooms into the "gray area" {1/ 0 < | | < 1/ 0 } shown in Figure 1. Note the intricate structure that becomes more prominent for larger values of . Both figures clearly demonstrate the disconnectedness of the basin of attraction of the roots of . We chose the values of = 128 and = 1024 since the "holes" become detectable with a International Journal of Computational Mathematics 13  naked eye around = 120 and we could zoom into them, and the larger value to demonstrate how much more the structure develops as increases.

Conclusions and Outstanding Questions
In the previous sections we have analyzed the behavior of Laguerre's method applied to the polynomials ( ) = − 1 in the extended complex plane and provided computational results demonstrating the interesting behavior of the method. We now have an almost complete understanding of the behavior of the method outside of the two gray areas that contain the boundary of the basin of convergence. We concluded that for initial guesses in̂0 ∪̂1 ∪ 1 the method converges to a root of (Proposition 10), and for initial guesses in 0 ∪ 0 the method converges to the two-cycle {0, ∞} (Proposition 7).
The numerical results indicate that the basin of attraction of the roots and the basin of attraction of the two-cycle share a common boundary, which should then be an invariant set under the Laguerre iteration function (5) and consist only of finite cycles and infinite orbits. We have not pursued this direction in great depth, as it would likely require extending the theory of Julia and Fatou sets [16] to functions that are not rational. Note that the Laguerre iteration function (5) is not rational even if is even due to the choice of sign in the denominator of the method. We have, however, attempted to find some short cycles, other than those given in Proposition 12, numerically in the following way. First, we used the computational software program Mathematica to generate the contour plots of Re ( ( ) − ) = 0 and Im ( ( ) − ) = 0 in the sector 0 < < / (recall the symmetries in Propositions 1-3), and used the visually discovered points of intersection as initial guesses for the command FindRoot in Mathematica applied to ( ) − . This way we have been able to find some 2-, 4-, and 6-cycles for polynomials of low degrees. In particular, it appears that 2-cycles in the sector 0 < < / only exist for ≥ 10 with 10 -16 having one such 2-cycle each; 17 -26 having two; 27 -38 having three, and so forth. Regarding 4-cycles, we found two for = 5, 6, 7; four for = 8; eight for = 9; nine for = 10; ten for = 11 and 12, and so forth. Finally, 6-cycles appear to start at = 6; we found ten of them for = 6, twelve for = 7, and twenty-three for = 8. Not surprisingly, we have not found any short oddcycles, which seems reasonable due to the expected behavior of points near getting mapped close to and vice versa. We list the found 4and 6-cycles for = 5, 6, 7, 8 in Table 1, where all numbers have been computed to 16 significant digits accuracy.
Many questions remain. What is the shape of the boundary of the basin of convergence? We see in Figures 2 and  3 that the boundary appears to track and , but it does not coincide with these sets. The boundary is fractal (Figures 5 and 6) and, moreover, has many other components in the annuli determined by 0 and 0 (Figures 7 and 8). It may be of interest to see whether a fractal dimension of the boundary has a simple dependence on . We speculate that the dimension might grow from 1 to 2 as increases from 5 to ∞, but we have not pursued this idea further.
It would also be interesting to see if other families of polynomials exhibit similar features to those observed for ( ) = − 1.
In particular, what determines the size and shape of the basins of convergence to the roots? Is it due to the symmetry of the roots that the measure of the basins of convergence tends to zero? If so, would other symmetric arrangements of the roots yield similar results? Perhaps the questions should be reversed. Are there families of polynomials for which Laguerre's method converges to a root except if starting from a set of zero measure? If so, what are they? We intend to look into some of these questions in future work.
We conclude with the following interesting observation. The fact that the method theoretically converges only in the small annulus in the neighborhood of the unit circle 1 suggests that Laguerre's method is unsuitable practically and raises a valid concern for general polynomials. On the other hand, when the method is implemented in its general formulation (2) and applied to polynomials ( ) = − 1, the resulting image of the basins of attraction may look like Figure 9, in which the basin of attraction of the roots of 7 is shown as computed on a 1000 × 1000 grid of points. Note that visually the method appears to converge from any point in the displayed squares, which is not the case when the mathematically equivalent formulation (5) is used. We allow the iterations to run to convergence (or a prescribed maximum number of iterations, 100 in this computation) and observe convergence to apparently random roots even for initial guesses in the set of theoretical convergence to the two-cycle {0, ∞} (compare to Figure 2(c), where the basin of attraction of {0, ∞} is colored white). The reason for this peculiar behavior is the loss of significance in the computation of the expression ( − 1) 2 ( ( )) 2 − ( − 1) ( ) ( ) in the denominator of (2). Note that both terms in the difference have leading terms 2 ( − 1) 2 2 −2 , and the actual difference should be equal to 2 ( − 1) 2 −2 . We therefore see that, for large | |, significant errors will occur in the computation of the square root in (2). In fact, the relative error in the computation of the square root is roughly proportional to √ (1 + | | ), where is the machine epsilon, so, for example, with = 8 and the usual 64-bit double precision, the relative error is on the order of 1 with | | as small as 100. Such errors then lead to the subsequent iterates attaining values from which convergence occurs. We also note that such loss of significance will occur for any polynomial ( ) of degree and | | large enough, since for a general polynomial of degree the difference ( − 1) 2 ( ( )) 2 − ( − 1) ( ) ( ) will have a leading term of order 2 −4 , two orders of magnitude smaller than the leading terms of ( − 1) 2 ( ( )) 2 and ( − 1) ( ) ( ). Perhaps this observation helps explain the popular notion that Laguerre's method seems to converge to a root from almost any initial guess.