NAG CL Interface
c02afc (poly_complex)

Note: this function is deprecated. Replaced by c02aac.

Old Code

Complex a[n+1], z[n];
Nag_Boolean scale;
...
scale = Nag_True;
nag_zeros_poly_complex(n, a, scale, z, &fail);

New Code

Complex a[n+1], z[n];
Integer itmax;
Nag_Root_Polish polish;
double berr[n], cond[n];
Integer conv[n];
...
itmax = 30;
polish = Nag_Root_Polish_Simple;
nag_zeros_poly_complex_fpml(a, n, itmax, polish, z, berr, cond, conv, &fail);

that the roots may be returned in a different order in array z.

Keyword Search:

NAG Library Manual, Mark 28.6

Interfaces: FL CL CPP AD PY MB

NAG CL Interface Introduction

C02 (Zeros) Chapter Contents

C02 (Zeros) Chapter Introduction

c02af: FL CL CPP AD PY MB

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

▸▿ 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

Settings help

CL Name Style:
Short (a00aac)
Long (impl_details)
Full (nag_info_impl_details)

1 Purpose

c02afc finds all the roots of a complex polynomial equation, using a variant of Laguerre's method.

2 Specification

#include <nag.h>

void	c02afc (Integer n, const Complex a[], Nag_Boolean scale, Complex z[], NagError *fail)

The function may be called by the names: c02afc, nag_zeros_poly_complex or nag_zeros_complex_poly.

3 Description

c02afc attempts to find all the roots of the

n

th degree complex polynomial equation

P (z) = a_{0} z^{n} + a_{1} z^{n - 1} + a_{2} z^{n - 2} + \dots + a_{n - 1} z + a_{n} = 0 .

The roots are located using a modified form of Laguerre's method, originally proposed by Smith (1967).

The method of Laguerre (see Wilkinson (1965)) can be described by the iterative scheme

L (z_{k}) = z_{k + 1} - z_{k} = \frac{- n P (z_{k})}{P^{'} (z_{k}) \pm \sqrt{H (z_{k})}},

where

H (z_{k}) = (n - 1) [(n - 1) {(P^{'} (z_{k}))}^{2} - n P (z_{k}) P^{''} (z_{k})]

, and

z_{0}

is specified.

The sign in the denominator is chosen so that the modulus of the Laguerre step at

z_{k}

, viz.

| L (z_{k}) |

, is as small as possible. The method can be shown to be cubically convergent for isolated roots (real or complex) and linearly convergent for multiple roots.

The function generates a sequence of iterates

z_{1}, z_{2}, z_{3}, \dots,

such that

| P (z_{k + 1}) | < | P (z_{k}) |

and ensures that

z_{k + 1} + L (z_{k + 1})

‘roughly’ lies inside a circular region of radius

| F |

about

z_{k}

known to contain a zero of

P (z)

; that is,

| L (z_{k + 1}) | \leq | F |

, where

F

denotes the Fejér bound (see Marden (1966)) at the point

z_{k}

. Following Smith (1967),

F

is taken to be

\min (B, 1.445 n R)

, where

B

is an upper bound for the magnitude of the smallest zero given by

B = 1.0001 \times \min (\sqrt{n} L (z_{k}), | r_{1} |, {| a_{n} / a_{0} |}^{1 / n}),

r_{1}

is the zero

X

of smaller magnitude of the quadratic equation

(P^{''} (z_{k}) / (2 n (n - 1))) X^{2} + (P^{'} (z_{k}) / n) X + \frac{1}{2} P (z_{k}) = 0

and the Cauchy lower bound

R

for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation

| a_{0} | z^{n} + | a_{1} | z^{n - 1} + | a_{2} | z^{n - 2} + \dots + | a_{n - 1} | z - | a_{n} | = 0 .

Starting from the origin, successive iterates are generated according to the rule

z_{k + 1} = z_{k} + L (z_{k})

for

k = 1, 2, 3, \dots

and

L (z_{k})

is ‘adjusted’ so that

| P (z_{k + 1}) | < | P (z_{k}) |

and

| L (z_{k + 1}) | \leq | F |

. The iterative procedure terminates if

P (z_{k + 1})

is smaller in absolute value than the bound on the rounding error in

P (z_{k + 1})

and the current iterate

z_{p} = z_{k + 1}

is taken to be a zero of

P (z)

. The deflated polynomial

\tilde{P} (z) = P (z) / (z - z_{p})

of degree

n - 1

is then formed, and the above procedure is repeated on the deflated polynomial until

n < 3

, whereupon the remaining roots are obtained via the ‘standard’ closed formulae for a linear (

n = 1

) or quadratic (

n = 2

) equation.

4 References

Marden M (1966) Geometry of polynomials Mathematical Surveys 3 American Mathematical Society, Providence, RI

Smith B T (1967) ZERPOL: a zero finding algorithm for polynomials using Laguerre's method Technical Report Department of Computer Science, University of Toronto, Canada

Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford

5 Arguments

1: $n$ – Integer Input: On entry: the degree of the polynomial, $n$ .

Constraint: $n \geq 1$ .
2: $a [n + 1]$ – const Complex Input: On entry: $a [i] . re$ and $a [i] . im$ must contain the real and imaginary parts of $a_{i}$ (i.e., the coefficient of $z^{n - i}$ ), for $i = 0, 1, \dots, n$ .

Constraint: $a [0] . re \neq 0.0$ or $a [0] . im \neq 0.0$ .
3: $scale$ – Nag_Boolean Input: On entry: indicates whether or not the polynomial is to be scaled. The recommended value is Nag_TRUE. See Section 9 for advice on when it may be preferable to set $scale = Nag_FALSE$ and for a description of the scaling strategy.
4: $z [n]$ – Complex Output: On exit: the real and imaginary parts of the roots are stored in $z [i] . re$ and $z [i] . im$ respectively, for $i = 0, 1, \dots, n - 1$ .
5: $fail$ – NagError * Input/Output: The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $⟨ value ⟩$ had an illegal value.
NE_COMPLEX_ZERO: On entry, the complex variable $a [0]$ has zero real and imaginary parts.
NE_INT_ARG_LT: On entry, $n = ⟨ value ⟩$ .
Constraint: $n \geq 1$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
NE_POLY_NOT_CONV: The iterative procedure has failed to converge. This error is very unlikely to occur. If it does, please contact NAG immediately, as some basic assumption for the arithmetic has been violated.
NE_POLY_OVFLOW: The function cannot evaluate $P (z)$ near some of its zeros without overflow. Please contact NAG immediately.
NE_POLY_UNFLOW: The function cannot evaluate $P (z)$ near some of its zeros without underflow. Please contact NAG immediately.

7 Accuracy

All roots are evaluated as accurately as possible, but because of the inherent nature of the problem complete accuracy cannot be guaranteed.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.

c02afc is not threaded in any implementation.

9 Further Comments

scale = Nag_TRUE

, then a scaling factor for the coefficients is chosen as a power of the base

b

of the machine so that the largest coefficient in magnitude approaches

thresh = b^{e_{\max} - p}

. You should note that no scaling is performed if the largest coefficient in magnitude exceeds

thresh

, even if

scale = Nag_TRUE

. (

b

e_{\max}

and

p

are defined in Chapter X02.)

However, with

scale = Nag_TRUE

, overflow may be encountered when the input coefficients

a_{0}, a_{1}, a_{2}, \dots, a_{n}

vary widely in magnitude, particularly on those machines for which

b^{4 p}

overflows. In such cases, scale should be set to Nag_FALSE and the coefficients scaled so that the largest coefficient in magnitude does not exceed

b^{e_{\max} - 2 p}

Even so, the scaling strategy used in c02afc is sometimes insufficient to avoid overflow and/or underflow conditions. In such cases, you are recommended to scale the independent variable

(z)

so that the disparity between the largest and smallest coefficient in magnitude is reduced. That is, use the function to locate the zeros of the polynomial

d \times P (c z)

for some suitable values of

c

and

d

. For example, if the original polynomial was

P (z) = 2^{−100} i + 2^{100} z^{20}

, then choosing

c = 2^{−10}

and

d = 2^{100}

, for instance, would yield the scaled polynomial

i + z^{20}

, which is well-behaved relative to overflow and underflow and has zeros which are

2^{10}

times those of

P (z)

If the function fails with NE_POLY_NOT_CONV, NE_POLY_UNFLOW or NE_POLY_OVFLOW, then the real and imaginary parts of any roots obtained before the failure occurred are stored in z in the reverse order in which they were found. More precisely,

z [n - 1] . re

and

z [n - 1] . im

contain the real and imaginary parts of the 1st root found,

z [n - 2] . re

and

z [n - 2] . im

contain the real and imaginary parts of the 2nd root found, and so on. The real and imaginary parts of any roots not found will be set to a large negative number, specifically

- 1.0 / (\sqrt{2.0} \times nag_real_safe_small_number)

10 Example

To find the roots of the polynomial

a_{0} z^{5} + a_{1} z^{4} + a_{2} z^{3} + a_{3} z^{2} + a_{4} z + a_{5} = 0

, where

a_{0} = (5.0 + 6.0 i)

a_{1} = (30.0 + 20.0 i)

a_{2} = - (0.2 + 6.0 i)

a_{3} = (50.0 + 100000.0 i)

a_{4} = - (2.0 - 40.0 i)

and

a_{5} = (10.0 + 1.0 i)

c02af: FL CL CPP AD PY MB

NAG CL Interfacec02afc (poly_​complex)

Old Code

New Code

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG CL Interface
c02afc (poly_complex)