Algorithm of polynomial factorization and its implementation in Maple Текст научной статьи по специальности «Математика»

MSC 47A68 DOI: 10.14529/mmp180408

ALGORITHM OF POLYNOMIAL FACTORIZATION AND ITS IMPLEMENTATION IN MAPLE

V.M. Adukov, South Ural State University, Chelyabinsk, Russian Federation, [email protected]

In the work we propose an algorithm for a Wiener-Hopf factorization of scalar polynomials. The algorithm based on notions of indices and essential polynomials allows to find the factorization factors of the polynomial with the guaranteed accuracy. The method uses computations with finite Toeplitz matrices and permits to obtain coefficients of both factorization factors simultaneously. Computation aspects of the algorithm are considered. An a priory estimate for the condition number of the used Toeplitz matrices is found. Formulas for computation of the Laurent coefficients with the given accuracy for functions that analytical and non-vanishing in an annular neighborhood of the unit circle are obtained. Stability of the factorization factors is studied. Upper bounds for the accuracy of the factorization factors are established. All estimates are effective. The proposed algorithm is implemented in Maple computer system as module "PolynomialFactorization". Numerical experiments with the module show a good agreement with the theoretical studies.

Keywords: Wiener-Hopf factorization; polynomial factorization; Toeplitz matrices.

Introduction

Mathematical modelling of wave diffraction, problems of dynamic elasticity and fracture mechanics, and geophysical problems is reduced to so-called Wiener-Hopf factorization problem for matrix functions [1-3]. The factorization of matrix functions is also a powerful tool used in various areas of mathematics [4,5]. However for the matrix case there is no constrictive solution of the factorization problem in a general setting. Moreover, it is difficult to develop approximate methods for the problem since the matrix factorization is unstable. For this reason it is very important to find cases when the problem can be effectively or explicitly solved. The current state of this problem is presented in [6]. The first stage of an explicit method for solving of the factorization problem is the factorization of scalar functions. In turn, this problem can be reduced to the polynomial factorization [7].

Let p(z) = p0 + p1z + • • • + pvzv be a complex polynomial of degree v > 1 and p0 = 0,pv = ^^ume p(z) = 0 on the unit circle T, hence p(z) = 0 on a closed circular annulus K := {z E C : r < \z\ < R} for some 0 < r < 1 < R < to. By k denote the number of zeros of p(z) inside the unit circle. Let j j = 1,..., v, be the zeros of p(z) and

0 < \&\< ... < \Ck \ <r< 1 <R< \£K+I\ < ■■■ <\£v\. Here the zeros are counted according to their multiplicity. Denote

pi(z) := (z - • • • (z - iH), p2(z) := (z - Ck+1 ) • • • (z - iv). (1)

The representation

p(z )= p-(z)zKp+(z), \z\ = 1, (2)

where p-(z) = p+(z) = p2(z), is called the Wiener-Hopf factorization, and

p(z )= pi(z)p2(z) (3)

p(z) p(z)

factors as output data of the problem. The naive method of the polynomial factorization is to use formulas (1). However it is well known that roots of a polynomial are in general not well-condition functions of its coefficients (see, e.g., [8]), and coefficients of a polynomial are also not well-condition functions of its roots [9]. The latter means that, in general, we can not solve numerically the polynomial factorization problem by the naive way.

Nevertheless there exist numerical methods for solving of the problem. The basic works in this direction are cited in [10].

In this work we propose approach that is close to Algorithm 3 of D.A. Bini and

pi (z)

p2(z) simultaneously. In addition, we find effective upper bounds for the accuracy of the factorization factors. Note that our technique can be extended to the factorization problem for analytic functions. In this case we can obtain all coefficients of the polynomial factor and a required number of Taylor coefficients of the analytic factor.

1. Preliminaries

Throughout this paper, ||x|| means the Holder 1-norm ||x|| = |xi| + • • • + \xk\, where x = (x1;... ,xk)T E Ck. For a matrix A E Cixk we always use the induced norm ||A|| = maxi<i<^l=1 \Aij \.

Respectively, the norm of a polynomial p(z) = p0 + p1z + • • • + pvzv is the norm of the vector (p0,p1,... ,pv)T. For p(z) we will also apply the maximum norm ||p||C = maxxgf |p(z)| on the unit circle T.

The norms || • || and || • ||C are equivalent. CI early, ||p||C < ||p||. Since for p(z) it is fulfilled

2n

v 1 r , ___

the equality '^^k|2 = 2n P(e^)|2d^> we have ||p|| < Vv + 1 ||p|2 < Vv + 1 HpHc. k=o n 0

Thus,

Mc <|p|^vV + T Mc. (4)

In order to study stability of the factorization problem, we will need estimates for the norm of inverses of some Toeplitz matrices. Such estimates will be obtain in terms of ||p1| • ||p2||, where p1(z), p2(z) are the factorization factors of p(z). To get effective

| p 1 | | p 2 | | p|

Let q(z) = q1(z)q2(z), where q1(z) q2(z) are arbitrary monic complex polynomials and v = degq. It is obvious that ||q|| < ||q1| ||q2|. In the work [11], D.W. Boyd proved the following inequality ||q1||c ||q2|c < $v||q|c^ere 8 = e2G/n = 1, 7916228120695934247... and G is Catalan's constant. The inequality is asymptoticaly sharp as v ^

Taking into account (4), in our case we obtain

HpH < ||P1|| M < 8VV(k + 1)(v - k +1) HpI (5)

However, the exponential factor 8V can overestimate the upper bound. In some special cases we can obtain more precise estimates. For example, this can be done for a so-called

spectral factorization of polynomials. In order to take into account a specific character of p(z)

\\p\\ < \\pi\\ \\p2W < So\\p\\, (6)

instead of (5). Here 1 < So < Sv ^J(k + 1)(v - k + 1).

2. Basic Tools

Let M, N be integers, M < N, and cj = (cM,cM+l;... ,cN) a nonzero sequence of complex numbers. In this section we introduce notions of indices and essential polynomials for the sequence cThese notions were given in more general setting in the paper [12]. Here we will consider the scalar case only. The proofs of all statements of this section can be found in [12].

Let us form the family of Toeplitz matrices Tk(cj) = (ci-j) i=k,k+1,...,N, M < k < N,

j=0,l,...,k-M

and study the sequence of the spaces ker Tk (cj). Further it is more convenient to deal not with vectors Q = (q0,ql,... ,qk-M )T E ker Tk but with their generating polynomials Q(z) = q0 + ql z + • • • + qk-M zk-M. We will use the spaces Nk of the generating polynomials instead of the spaces ker Tk. The generating function Y^N=M ckzk of the sequenee cM will be denoted by cj (z).

Let us introduce a linear functional ct by the formula: a{zj} = c-j, -N < j < -M. The functional is defined on the space of rational functions of the form Q(z) = J2~=-N qjzj. Besides this algebraic definition of a we will use the following analytic definition

ct{Q(z)} = 2^1 t-lcN(t)Q(t) dt. (7)

r

Here r is any closed contour around the point z = 0.

Denote by Nk (M < k < N) the space of polynomials Q(z) with the formal degree k-M

a[z-iQ(z)} = 0,i = k,k + 1,...,N. (8)

It is easily seen that Nk is the space of generating polynomials of vectors in ker Tk. For convenience, we put NM-l = 0 and denote by Nn+1 the (N - M + 2)-dimensional space of all polynomials with the formal degree N - M + 1. If necessary, the more detailed notation Nk(cM) instead of Nk is used.

Let dk be the dimension of the space Nk and Ak = dk - dk-l (M < k < N + 1). The following inequalities are crucial for the further considerations: 0 = AM < AM+l < ... < an < an+i = 2.

It follows from the inequalities that there exist integers ¡il < fi2 such that

am = ... = Am =0,

Am+i = ... = A№ =1, (9)

AM2+i = ... = an+i = 2. If the second row in these relations is absent, we assume ¡il = fi2.

Definition 1. The integers ¡1, ¡2 will be called indices of the sequence c^.

It is easily seen that Nk and zNk are subspaces of Nk+1l M — 1 < k < Let hk+1 be the dimension of any complement Hk+1 of the subspace Nk + zNk in the whole space

Nk+1.

From (9) we see that hk+1 = 0 iff k = ¡j (j = 1, 2) hk+1 = 1 if ¡1 < ¡2, aM hk+1 = 2 for ¡1 = ¡2- Therefore, Nk+1 = Nk + zNk, for k = j and Nk+1 = (Nk + zNk) © Hk+1 for k = ¡j

Definition 2. Let ¡1 = ¡2. Any polynomials Q1(z), Q2(z) that form, a basis for the two-dimensional space N^1+1 are called the essential polynomials of the sequence cN corresponding to the index ¡1 = ¡2.

If ¡1 < ¡2, then any polynomial Qj (z) that is a basis for an one-dimensional complement +1 is said to be the essential polynomial of the sequence corresponding to the index ¡j, j = 1, 2.

3. Main Results

In this section we propose a method for solving the problem of the polynomial factorization in terms of indices and essential polynomials of some sequence. The proofs of the theorems of the section can be found in preprint [13].

Let p(z) = p0 + p1 z + • • • + pvzv be a polynomial of degree v > 1, p0 = 0, p(z) = 0 on the unit circle T, hence p(z) = 0 on a closed annulus K := {z E C : r < |z| < R} for some 0 < r < 1 < R < to. Let k = indTp(z) be the number of zeros of the polynomial inside the unit circle. We can assume that 0 < k < v, otherwise the factorizations are trivial. Put n0 = max{K, v — k}.

We will find the factorization of p(z) in the form (3): p(z) = p1(z)p2(z), where degp1 = k, degp2 = v — k, the polynomial p1(z) is monic, and all zeros of p1(z) (respectively p2(z)) lie in the domain {z E C : |z| < r} (respectively in {z E C : |z| > R}). This means that

p(z) = p-(z) zKp+(z), z E T,

is the Wiener-Hopf factorization of p(t) ^^^^^^^^^ by the condition p-(x>) = 1. Here p-(z) = ^zKr, P+(z) = p2(z). Let p-1(z) = Z-oo ckzk be the Laurent series for analytic function p-1(z) in ^te annulus K,

ck = J t-k-1p-1(t) dt, k E Z, r < p < R. (10)

|z|=P

Form the sequence c-n-K = (c-n-K,..., c-K,..., cn-K) for the given n > n0. It is easily to verify that the sequence is non-zero.

Theorem 1. For any n > n0 the integers —k, —k are the indices, and Q1(z) = zn-K+1p1(z), Q2(z) = p2(z) are essential polynomials of the sequence c-n-K-

From the theorem it follows that there exist the essential polynomials Q1(z),Q2(z) of the sequence c-n-K such that the following additional properties are fulfilled:

1. deg Q1(z) = n +1 Q1(0) = 0, Qm+1 = 1.

2. degQ2(z) <n +1 Q2(0) = 0 ^{zKQ2(z)} = 1.

Vice versa, if Qi(z),Q2(z) £ N-K+1 and satisfy conditions 1-2, then Qi(z),Q2(z) are the essential polynomials of the sequence C——-K.

Definition 3. Let n > n0. Polynomials Qi(z),Q2(z) £ N-K+1(c—K) satisfying conditions 1-2 will be called the factorization essential polynomials of C——-K.

Theorem 2. The factorization essential polynomials are uniquely determined by conditions 1-2. Let n > n0 + 1, suppose that R1(z), R2 (z) are any essential polynomials of the sequence C-—-K. Then the factorization essential polynomials of C——-K can be found by the formulas

Qi(z ) = -

Ri(z) R*(z) Rl,0 R2,0

Q2 (z) = -

^0

Rl(z) R2(z)

Ri R2

(11)

Here Rj,0 = Rj (0), Rj,n+l is the coefficient of zn+l in the polynomial Rj (z), j = 1, 2, and

Ri 0 R2 0 Rl R2 ,n+l

= 0, ao = a < z

{

Rl(z) R2(z ) Rl R2 ,n+l

}

= 0.

Now we can construct the Wiener-Hopf factorization of a polynomial with the help of the factorization essential polynomials.

Theorem 3. Let n > n0 and let Qi(z),Q2(z) be the factorization essential polynomials of the sequence Then the Wiener-Hopf factorization of the polynomial p(z) can be

constructed by the formula p(z) = p-(z)zKp+(z), where

P-(z) = z-n Ql(z), p+(z) = Q2(z).

(12)

The following theorem about explicit formulas for the factors of the Wiener-Hopf factorization of a polynomial p(z) is the main result of the section.

Theorem 4. The matrices T-K (o—K) are invertible for all n > n0.

Let n > n0 + 1. Denote by a = (al;..., an)T and ß = (ß0,..., ßn)T the solutions of the systems

T-K (o—l) a = -(c-n-K)T, T-K (o-nn--l) ß = el,

(13)

respectively. Here e 1 = (1,0,..., 0)T.

Then a 1 = ... = an-K = 0 ßo = 0 ßv-K+1 = • • • = ßn = 0, and the factors from the

p(z)

p-(z) = z-K(an-K+i + • • • + anzK 1 + zK), p+ (z) = f3o + Piz + • • • + f3—zV-K. 4. Some Computational Aspects of the Factorization Problem

In this section we consider some computation problems that arise in a numerical implementation of the results of Section 3.0ur aim is to obtain the factors p1(z), p2(z) with the guaranteed accuracy. The proofs of the statements of the section can be found in preprint [13].

x

4.1. An a Priori Estimate of the Condition Number for the Factorization Problem

Theorem 4 shows that solving of the factorization problem is equivalent to solving of linear systems with the invertible matrix T-k(c---k), n > n0 := max{K, v — k}.

Here we obtain an upper bound for the condition number k(T-K(c---K)) = \\T-k(c---k)|| WTZ^c—n—x) || in terms of the given polynomial p(z).

Denote mK = minz&K |p(z)|, m1 = min|^|=1 |p(z)|, p = max{r, 1/R}. Then the following estimate of the condition number is fulfilled.

nv

k(T-K(c—J) < m,JA J1±lÇ^t1^ ||p||) . (14)

[m,K (1 - p) mi J

4.2. Computation of the Laurent Coefficients of Analytic Functions

To realize the factorization method proposed in Section 3 we must calculate the Laurent coefficients c-n-K ,...,c-K,... ,cn-K1 of the func tion p-i(z) for n > n0 = max{x, v — k}.

In general, the coefficients can be found only approximately. In order to do this, we will apply the method suggested by D.A. Bini and A. Böttcher (see Theorem 3.3 in [10]). For future applications we will consider more general situation than in this work. Moreover, our proof differs from the proof in the above mentioned paper.

Let f (z) be a function that analytic in the annulus K = {z G C : r < \z\ < R}, 0 <

r < 1 < R < to. By fk denote the Laurent coefficients of f (z): fk = f t-k-if (t) dt,

m |i|=p

r < p < R. For l,k G Z, I > 2, define fk(I) = jYj=o , where = e^j are the zeros of the polynomial zl — 1.

Theorem 6. Let MK = max \f (z)|, p = max{r, 1/R}, and 1 be an even positive integer. Then

\fk - fk(l)| < pl/2 (15)

(1 - Pl)

fork = -¿/2,..., 0,...,l/2.

By the theorem, in order to compute every element of the sequence fM, fM+i,...,fN with the given accuracy, we have to select an appropriate number I.

4.3. Stability of the Factors p1(z), p2(z)

Now we study the behavior of the factors p1; p2 under small perturbations of p(z). Let m1 = min |p(z)|. It is easily seen that pC(z) = 0 on T and indT p(z) = indT p(z) if

\\p — p\\ < шь Let p(z) = p1(z)p2(z) be the factorization of p (z). By Cj we denote the Laurent coefficients of p-1(z). Estimate \\p1 — pC1\ \\p2 — pC2^ia \\p — pC\\.

Theorem 7. Let n > no + 1. // \\p — p\\ < min | , } • then

¿o\\p\\ 1 + p

.. .... 4(2n +1)So \ \\Pi - Pi\\ <-

mf

and

mK 1 — p

+ 1

\\P — p\\, (16)

4(2n + 1)«

2

\\P2 — P2\ < --m2 \\p — P \\- (17)

mi

Now we consider perturbations of the polynomials p1(z), p2(z) caused by the

approximation of the sequence by where ck = j=0 g^1)^ • Le^ Pi(z),

p2(z) be polynomials that define by Eq. (13) for the sequence Ic-n-x-

Theorem 8. Let n > n0 + 1, £ № an even integer such that 1 > 2(n + x), and

pm <^rjm^nr. (is)

1 — pl 2(4n + 2)So

Then

.. .... 2(4n — 2)5o \ \\Pi — Pi\\ <~

mK

2So(1 + r,nt-1, + 1 (1 — p)mK .

pl/2 1 — pl

< 2(4n + 2)S02\p\2 pl/2

\\p2—P2\ <—mK— T—p ■

From the theorem it is easy to find the estimate of 1 in order to obtain the factors p1(z) p2(z) with the given accuracy e.

Corollary 1. Let 0 < e < 1 and

a = e " '> 0 + • «4» + 2*w} .

If I is an even integer such that

l> 2ma^< n + x,--,-L } , (19)

1 I logP\ '

then \\pi — pi\\ < e, \\P2 — P2\\ < e.

5. Algorithm and Its Implementation 5.1. Algorithm

The above results can be summarized in the form of the following algorithm. Algorithm. Polynomial factorization

INPUT. The polynomial p, the parameter par giving a variant of estimate of \\pi\ • \\p2\\ via \\p\ (par =1, 2, 3), the integer n > degp, the accuracy A = 10"9 for p: \\p — p\\ < A.

pi p2 e e

COMPUTATION.

1. Compute the index k of p.

2. Compute the radii r, R of the circular annulus K and p = max{r, 1/R}, compute m1 = min|z|=1 |p(z)|, mK = min^K |p(z)|.

3. Find accuracy e1; e2 for p1; p2 by formulas (16), (17). Compute the theoretically guaranteed accuracy e := max{e1 ,e2}. Define d such that e < 10_d.

4. To compensate the loss of accuracy under calculations with the Toeplitz matrix T-k(c___k) we must improve accuracy of Laurent coefficients ck- To do this we estimate the condition number k of the matrix by formula (14), find d such that k < 10d, and put e := 10-d-d.

5. Find an even integer 1 satisfying inequality (19) (in the inequality put e = e).

6. Form the sequence 'c---:,^.

7. Form the Toeplitz matrix T_k+1 (c___k) and find a basis {R1; R2} of its kernel. The last operation can be done with the help of SVD.

I ~n+1

uuiy ииш1шо wiyz ) — aiz I - - - I

во + в1 z + ••• + PnZn by (11).

8. Find the factorization essential polynomials Qi(z) = a\z + • • • + anzn + zn+\ Q2(z)

z n

nz

9. Verify that the absolute values of the coefficients a\,..., an-x, (3v-K+i,..., ¡3n are less than e and delete these coefficients (see Theorem 4).

10. pi(z) := z-n+K-1Qi(z), p2(z) := Q2(z) e.

11. end

5.2. Implementation in Maple

The algorithm was implemented in Maple as the module "PolynomialFactorization". The module can work in Maple 17. To access "PolynomialFactorization", use the commands

> read("PolynomialFactorization.mpl");

> with(PolynomialFactorization);

Then it is necessary to set the value of Digits taking into account initial data. If the user needs the results of intermediate calculations, then to set

> PrintOn:= true;

The basic routines of "PolynomialFactorization" are indpoly, annulusn (annulus), condnpoly, getellpoly, LaurentCoeff, factesspoly, SolverPF.

• indpoly returns the index k of the polynomial p found by formula (12.6) from [14].

• p mK

• condnpoly returns the estimate of the condition number k by formula (14).

• getellpoly returns the integer 1 required for calculation of the Laurent coefficients with the given accuracy e.

• LaurentCoeff returns the sequence cM of the Laurent coefficients of the function 1/p.

•

• SolverPF returns the factors p1; p2, and the guaranteed accuracy e of calculation.

The polynomial factorization is realized by routine SolverPF. The following input data are passed to SolverPF:

• a polynomial p in a variable z.

• the parameter par taking the values 1, 2, 3. A value of par gives a variant of estimate of (pi) • ||p21| via ||p||- If p = Yl2=0pjzj be a real polynomial of degree 2m such that p2m-j = pj for j = 0,... ,m, p0 = 1, and all roots of p(z) have negative real parts, then par = 1 If p = Y12=0 pjzj be a complex polynomial of degree 2m such that p2m-j = pj j = 0, . . . , m p0 = 1 par = 2 p

par = 3

• an integer n > degp, where 2n + 1 is the length of the used sequence c——n—K-

• the accuracy A = 10—'q of calculation p. Here is the example to call the program.

> Digits:= 15;

> p:= l+3iz/2+z~2;

> pi, p2, epsilon := SolverPF(p, 2, 3, 10"(-10)):

"Calculation in progress ...." "End of calculation"

> pi;

1.00000000000000*z-0.500000000000004*1

> p2;

1.00000000000001*z+l.99999999999999*1

> epsilon;

0.00000755534207942884 6. Numerical Examples

In the following examples we use module "PolynomialFactorization". All computations were performed on a desktop.

p(z)

and its Wiener - Hopf factorization is actually the spectral factorization.

Example 1. Let p(z) = (z + 1/2)(z + 1/3) • • • (z + 1/12)(z + 2)(z + 3) • • • (z + 12). Taking into account the values of the coefficients of p(z), we choose the precision Digits := 20. Assume that the accuracy of the input data A is equal to 10-15.

We have v = 22, k = 11 ||p|| = 20237817600. Computations show that p := 0, 51, m1 = 3,326340 x 106, mK = 30,448076. Put n = v + 1 = 23. Here par := 1 and S0 = 1. Calculation of the theoretically guaranteed accuracy e gives the following result e = 0, 695883 x 10-5. By formula (14), we obtain the following estimate k(T-K(c———K)) < 2, 859480 x 105. It follows from this that e = 10-22 and we get I = 136.

In this example the exact output is known

pi(z) = (z + 1/2)(z + 1/3) • • • (z + 1/12), P2(z) = (z + 2)(z + 3) • • • (z + 12).

Table 1 shows the results of approximate computations of the factors p^z), p2(z). It contains coefficients pjlk1 pj2k1 absolute errors \pj\ — plk\, \pj2k — pi\ for the coefficients p\, pi, and ||p1 — p1||, ||p2 — p2||. For pk, pi the number of decimal places obtained accurately is shown.

Table 1

Coefficients pk p2k

k pk \pk- pk 1 pk \pk- pk \

0 0 2,087675e-9 479001600,00000 l,04000e-9

1 0 l,60751e-7 1007441280,00000 l,26000e-8

2 0 0,55114e-5 924118272,00000 3,78100e-8

3 0,00011 l,97436e-18 489896616,00000 5,78000e-8

4 0,00145 4,98731e-17 167310220,00000 5,46800e-8

5 0,01300 9.511 !()c-17 38759930,00000 3,62700e-8

6 0,08091 1.20571c-16 6230301,00000 2,20830e-8

7 0,34928 l,18000e-16 696333,00000 l,81875e-8

8 1,02274 8,87000e-17 53130,00000 l,72958e-8

9 1,92925 4,76000e-17 2640,00000 l,24386e-8

10 2,10321 l,40000e-17 77,00000 2,80480e-9

11 1,00000 0 0,999999 9,23655e-9

llpi - Pill 0.567 13('-5

\\p2 - P2 1 2.822 16c-7

Thus ||p1 — p11| = 0, 56743 x 10-5 < 0, 695884 x 10-5 = e, and ||p2 — p21| = 2, 82246 x 10-7 < 0, 695884 x 10-5 = e. We obtain p1(z) p2(z) with the desired accuracy.

The following example was taken from [10]. Since p(z) has real coefficients pj and pv-j = pj, ^te factorization of p(z) is also the spectral factorization.

Example 2. Let p(z) = E1==o zi + 4z^ Digits := 20 A = 10-12. Now p = 0,83,

m1 = 1, 542464 mK = 0, 062855

In this example v =10 k = 5 ||p| = 15 n = v + 1 = 11 par := 2, 50 = k + 1 = 6. For the accuracy e we obtain e = 0, 536458 x 10-4. From formula (14) it follows the following estimate k{T-x(C---k)) < 1342, 00899^. ^to yields e = 10-17 and I = 418.

Table 2

Coefficients p\, pk

k 0 1 2 3 4 5

pk 0,23193 0,20715 0,17674 0,14253 0,10685 1,00000

pk 4,31154 0,46071 0,61452 0,76203 0,89314 1,00000

The computed coefficients of the factors pi(z), p2(z) are given by Table 2. We only indicate 5 decimal places here.

In order to verify the computation correctness of p2(z), we can use the following relation between the factors pi(z) and p2(z) in the spectral factorization: p2(z) = zKpi(1/z)/pi(0). For our example we have \\p2 — zKp1(1/z) /p1 (0)|| = 7, 36 x 10"18. Moreover, the residual error is \p1p2 — p\\ = 8,1 x 10"18.

p(z)

package Random Tools. Example 3. Let

Digits := 20, A = 10"18. Calculations show that p = 0, 943396, m1 = 2, 293009, mK = 0,66281.

Table 3

Coefficients plk1 pk

k Pk Pk

0 -0, 099841 - 0,150475« -5, 090491 - 10,133912«

1 -0, 236722 + 0,118527« -14,129949 + 0, 552043«

2 -0, 385402 - 0, 732498« -4, 543939 + 4, 838437«

3 1,00000 -7, 958489 + 1, 840704«

4 -5, 515909 + 9, 645327«

5 4,196252 + 7, 320240«

6 0, 930308 + 0, 031004«

7 -0,181264 + 0, 732498«

8 1,000000

For the polynomial we have v =11 k = 3, n = v + 1 = 12 ||p|| = 42, 442968. Since p(z) is a polynomial of general type, par := 3, S0 = 3663, 225630 is found by the formula S0 = §v\j(k + 1)(v — k + 1). By this reasons, we are forced to use more accurate input data. Then the guaranteed accuracy in the output is e = 0, 254667 x 10"4.

From formulas (14) we obtain the estimate k(T-K(c———K)) < 1342, 00899. Hence e = 10-2^d I = 2096.

The computed coefficients of the factors pi(z), p2(z) are given by Table 3.The residual error is \\p 1p2 - p\\ = 4, 354070 x 10"7.

Let p1(z), p2(z) be the factorization factors of p(z) obtained by the naive method (via the roots ofp(z)). Then \\p 1 - p1\ = 1.3 x 10"10, \\p2 - p2\\ = 5, 239393 x 10"8.

Acknowledgement. The program module "PolynomialFaetorization" was designed in collaboration with VI.M. Adukov. The author is grateful to him for his invaluable help in programming.

References

1. Daniele V.G., Zieh R.S. The Wiener-Hopf Method in Electromagnetics. ISMB Series, New Jersey, SciTech Publishing, 2014. DOI: 10.1049/SBEW503E

2. Abrahams I.D. On the Application of the Wiener-Hopf Technique to Problems in Dynamic Elasticity. Wave Motion, 2002, vol. 36, pp. 311-333. DOI: 10.1016/S0165-2125(02)00027-6

3. Lawrie J.B., Abrahams I.D. A Brief Historical Perspective of the Wiener-Hopf Technique. Journal of Engineering Mathematics, 2007, vol. 59, pp. 351-358. DOI: 10.1007/sl0665-007-9195-x

4. Clancey K., Gohberg I. Factorization of Matrix Functions and Singular Integral Operators. Basel, Birkäuser, 1987.

5. Gohberg I.C., Feldman I.A. Convolution Equations and Projection Methods for Their Solution. Providence, American Mathematical Society, 1974.

6. Rogosin S., Mishuris G. Constructive Methods for Factorization of Matrix-Functions. IMA Journal of Applied Mathematics, 2016, vol. 81, no. 2, pp. 365-391. DOI: 10.1093/imamat/hxv038

7. Kisil A.V. A Constructive Method for an Approximate Solution to Scalar Wiener-Hopf Equations. Proceedings of The Royal Society A Mathematical Physical and Engineering Sciences, 2013, vol. 469, no. 2154, article ID: 20120721. DOI: 10.1098/rspa.2012.0721

8. Gautschi W. On the Condition of Algebraic Equations. Numerische Mathematik, 1973, vol. 21, no. 5, pp. 405-424. D01:10.1007/BF01436491

9. Uhlig F. Are the Coefficients of a Polynomial Well-Condition Function of Its Roots? Numerische Mathematik, 1992, vol. 61, no. 1, pp. 383-393. DOI: 10.1007/BF01385516

10. Bini D.A., Böttcher A. Polynomial Factorization through Toeplitz Matrix Computations. Linear Algebra Application, 2003, vol. 366, pp. 25-37.

11. Boyd D.W. Two Sharp Inequalities for the Norm of a Factor of a Polynomials. Mathematika, 1992, vol. 39, no. 2, pp. 341-349. DOI: 10.1112/S0025579300015072

12. Adukov V.M. Generalized Inversion of Block Toeplitz Matrices. Linear Algebra Appliction, 1998, vol. 274, pp. 85-124.

13. Adukov V.M. On Wiener-Hopf Factorization of Scalar Polynomials. 2018. Available at: http://arxiv.org/abs/1806.01646.

14. Gakhov F.D. Boundary Value Problems. N.Y., Dover, 1990.

Received July 20, 2018

УДК 519.688

DOI: 10.14529/ mmp180408

АЛГОРИТМ ПОЛИНОМИАЛЬНОЙ ФАКТОРИЗАЦИИ И ЕГО ИМПЛЕМЕНТАЦИЯ В MAPLE

В.М. Адуков, Южно-Уральский государственный университет, г. Челябинск Российская Федерация, [email protected]

В работе предложен алгоритм факторизации Винера - Хопфа скалярных многочленов. Алгоритм, основанный на понятиях индексов и существенных многочленов, позволяет найти факторизационные множители многочлена с гарантированной точностью. Метод использует вычисления с конечными теплицевыми матрицами и дает возможность получить коэффициенты обоих факторизационных факторов одновременно. Рассмотрены вычислительные аспекты алгоритма. Найдена априорная оценка числа обусловленности используемой теплицевой матрицы. Получены формулы для вычисления лорановских коэффициентов с заданной точностью для функций аналитических и не обращающихся в нуль в кольцевой окрестности единичной окружности. Изучена устойчивость факторизационных множителей. Установлены верхние границы точности вычисления факторизационных множителей. Все оценки являются эффективными. Предложенный алгоритм был реализован в компьютерной системе Maple в виде модуля «PolynomialFactorization». Численные эксперименты с модулем показали хорошее согласие с теоретическим исследованием.

Ключевые слова: факторизация Винера - Хопфа; полиномиальная факторизация; теплицевы матрицы.

Виктор Михайлович Адуков, доктор физико-математических наук, профессор, кафедра «Математический анализ и методика преподавания математики>, ЮжноУральский государственный университет (г. Челябинск, Российская Федерация), [email protected].

Поступила в редакцию 20 июля 2018 г.