Orthogonal Laurent polynomials and quadratures on the unit circle and the real half-line.

Abstract. The purpose of this paper is the computation of quadrature formulas based on Laurent polynomials in two particular situations: the Real Half-Line and the Unit Circle. Comparative results and a connection with the split Levinson algorithm are established. Illustrative numerical examples are approximate integrals of the form

[[integral].sup.1.sub.-1] (f(x) / [(x + [lambda])sup.r] [omega](x)dx, r = 1,2,3, ...

with f (x) a continuous function on [-1, 1], [omega](x) [greater than or equal to] 0 a weight function on this interval and [lambda] [member of] R such that [absolute value of [lambda]] > 1 is required. Here the classical Gaussian quadrature is an extremely slow procedure.

Key words. orthogonal Laurent polynomials, L-Gaussian quadrature, Szego quadrature, three-term recurrence relations, split Levinson algorithm, numerical quadrature.

AMS subject classifications. 41A55, 33C45, 65D30.

1. Orthogonal Laurent Polynomials. The important role played by the theory of orthogonal polynomials in the construction of an n-point quadrature rule

[I.sub.n] (f) = [n.summation over (j=1)][A.sub.j]f([x.sub.j]) (1.1)

to approximate integrals of the form

[I.sub.[mu]] (f) = [[integral].sup.b.sub.a] f (x) d[mu](x), -[integral] [less than or equal to] a < b [less than or equal to] + [integral] (1.2)

([mu] being a positive measure) so that [I.sub.n] (f) exactly integrates polynomials of degree as large as possible, is well known. Indeed, let [[psi].sub.n] (x) be the n-th orthonormal polynomial with respect to the measure [mu] and let [x.sub.1], ..., [x.sub.n] be its n distinct zeros. Then positive weights [A.sub.1], ..., [A.sub.n] can be uniquely determined so that the corresponding quadrature formula [I.sub.n](f) satisfies [I.sub.n] (P) = [I.sub.[mu]] (P) for all P [member of] [[PI].sub.2n-1] ([[PI].sub.k]: space of polynomials of degree k at most, [PI]: space of all polynomials). This is the well known n-point Gauss-Christoffel (or Gaussian) formula for the measure [mu] on (a, b) (for a comprehensive survey on this topic see [11]). Furthermore, these formulas are "optimal" in the sense that there is not an n-point rule (1.1) which is in [[PI].sub.2n]. It holds

[A.sub.j] = 1 / [[summation].sup.n-1.sub.k=0] [[psi].sup.2.sub.k]([x.sub.j]), j = 1, ..., n. (1.3)

The numerical power and effectiveness of such rules when dealing with smooth integrands f (x) has been convincingly demonstrated. However, convergence can become extremely slow if the integral has singularities near the integration interval (a, b). In order to overcome this drawback, quadrature formulas exactly integrating rational functions with prescribed poles outside (a, b) have been considered as an alternative in the last years (see e.g. [12], [23] or [3]). Thus, in this paper we shall be concerned with quadrature formulas (1.1) exactly integrating rational functions with all their poles at the origin and the infinity (Laurent polynomials) in order to estimate an integral like [[integral].sub.S] f (z)d[mu](z), [mu] being a positive measure supported on S [subset] C with either S = (a, b) , 0 [less than or equal to] a < b [less than or equal to] +[infinity] or S = T = {z [member of] C : [absolute value of z] = 1} (the unit circle). Here the theory of orthogonal Laurent polynomials will represent the basic ingredient in the construction of the corresponding quadrature rules.

Thus, in general, let [mu] be a positive Borel measure in the complex plane and consider the Hilbert space [L.sub.2] ([mu]) of measurable functions [phi](z) for which [integral][[absolute value of [phi](z)].sup.2] d[mu](z) < +[infinity]. As usual, in [L.sub.2] ([mu]) we can define the inner product

< [phi], [psi]) >= < [phi], [psi] > [sub.[mu]] = [integral] [phi](z) [bar.[psi](z)]d[mu](z), [phi], [psi] [member of] [L.sub.2] ([mu]) (1.4)

Suppose that [{[[phi].sub.k](z)}.sup.n.sub.k=0] is a system of linearly independent functions in [L.sub.2]([mu]). By applying the Gram-Schmidt orthogonalization process to [{[[phi].sub.k](z)}.sup.n.sub.k=0] a new system of linearly independent functions [{[[psi].sub.k](z)}.sup.n.sub.k=0] can be obtained such that [[psi].sub.n] (z) is a linear combination of the n + 1 functions [{[[phi].sub.k](z)}.sup.n.sub.k=0] and

< [[psi].sub.n](z), [[psi].sub.m](z), [[psi].sub.m](z) >= [integral] [[psi].sub.n](z) [bar.[[psi].sub.m](z)]dz = [k.sub.n][[delta].sub.n,m], [k.sub.n] > 0.

When the process is repeated for each natural n, an essentially unique system [{[[psi].sub.n](z)}.sup.[infinity].sub.n=0] of orthogonal functions with respect to the measure [mu] is obtained. If moreover

[[parallel][[psi].sub.n](z)[parallel].sup.2] = [integral] [[absolute value of [[psi].sub.n](z)].sup.2]d[mu](z) = 1,

then [{[[psi].sub.n](z)}.sup.[infinity].sub.n=0] is called an orthonormal system. For our purposes we will concentrate on the linearly independent system of monomials both with positive and negative exponents, say ..., [z.sup.-2], [z.sup.-1], 1, z, [z.sup.2] , ... yielding the so-called Laurent polynomials. Indeed, for p and q integers such that p [less than or equal to] q, [[DELTA]sub.p,q] will denote the space of Laurent polynomials of the form L (z) = [[summation].sup.q.sub.j=p] [[alpha].sub.j] [z.sup.j] with [[alpha].sub.j] [member of] C and [DELTA] the space of all Laurent polynomials. Observe that [[DELTA].sub.0,x] = [[PI].sub.k].

Now, in order to generate a sequence of nested subspaces of Laurent polynomials similar to the sequence of subspaces [{[[PI].sub.k]}.sub.k[greater than or equal to]0], we will start from two nondecreasing sequences of nonnegative integers [{p(n)}sup.[infinity].sub.n=0] and [{q(n)}.sup.[infinity].sub.n=0] such that p(n) + q(n) = n for all n = 0,1,2, ... and set

[L.sub.n] = [[DELTA].sub.-p(n),q(n)] = Span {[z.sup.j] ; -p(n) [less than or equal to] j [less than or equal to] q(n)}.

Observe that [L.sub.0] = Span{1}, dim ([L.sub.n]) = n + 1 and that [L.sub.n] [subset] [L.sub.n+1]. Furthermore if p(n) = 0 then q(n) = n and [L.sub.n] = [[DELTA].sub.0,n] = [[PI].sub.n]. In the sequel, in order to guarantee that [[union].sup.[infinity].sub.n=0] [L.sub.n] = [DELTA], we will assume that [lim.sub.n[right arrow]]p(n) = [lim.sub.n[right arrow][infinity]] q(n) = [infinity] and further restrict ourselves to the particular case

p(n) = E [n + 1 / 2], q(n) = n - p(n) = E [n/2], n [greater than or equal to] 0, (1.5)

where, as usual, E[x] denotes the integer part of x. In other words, the sequence {p(n)} has induced the following "ordering" in [DELTA]:

[[DELTA].sub.0,0], [[DELTA].sub.-1,0], [[DELTA].sub.-1,1], [[DELTA].sub.-2,1] , ....

Thus, according to this ordering we can construct a unique sequence (up to a sign) [{[[psi].sub.n](z)}.sup.[infinity].sub.n=0] of Laurent polynomials satisfying:

1. [[psi].sub.n] [member of] [L.sub.n]\[L.sub.n-1], n = 1, 2, ....

2. < [[psi].sub.n] (z), R > [sub.[mu]] = 0 for all R [member of] [L.sub.n-1] (i.e., [[psi].sub.n] [perpendicular to] [L.sub.n-1]).

3. [[parallel][[psi].sub.n](z)[parallel].sup.2.sub.[mu]] =< [[psi].sub.n](z), [[psi].sub.n](z) > [sub.[mu]] = 1.

In this case, [{[[psi].sub.n](z)}.sup.[infinity].sub.n=0] will be called an orthonormal Laurent polynomial sequence corresponding to the measure p and the ordering induced by p(n) = E [n+1/2] which has been the most extensively treated case in the literature. For other choices of p(n), see, e.g., [27], [4] or [30]. As already said, in the rest of the paper we will concentrate on measures [mu] supported not on the whole complex plane C but either on intervals (a, b) of the positive real half-line or on the unit circle T. In the two following sections the differences and similarities for sequences of orthogonal Laurent polynomials concerning both situations will be displayed. In this respect and in order to fix a uniquely determined sequence of orthogonal Laurent polynomials, we will use the concept of leading coefficient associated with a Laurent polynomial (see [7]). Thus, for a given Laurent polynomial R there exists a unique natural number n such that R [member of] [L.sub.n] and R [not member of] [L.sub.n-1]; n is called the "L-degree" of R for the ordering established above. Observe that if P(z) is a polynomial of degree n, then its L-degree is just 2n. Indeed, P(z) [member of] [[DELTA].sub.-n,n] = [L.sub.2n] and P(z) [not member of] [L.sub.2n-1]. Setting R(z) = [[summation].sup.q(n).sub.j=-p(n)] [[alpha].sub.j][z.sup.j], one defines the leading coefficient of R as [[alpha].sub.q(n)] if n is even or [[alpha].sub.-p(n)] if n is odd. Clearly, for an orthogonal sequence [{[[phi].sub.n](z)}.sup.[infinity].sub.n=0] of Laurent polynomials, [[phi].sub.n] (z) has a non-zero leading co-efficient for each n, which will be sometimes assumed positive. When equal to one, [[phi].sub.n] (z) is said monic. Clearly, an orthonormal sequence ([[psi].sub.n](z)} of Laurent polynomials with positive leading coefficient for each n is uniquely determined.

The paper is organized as follows: In Section 2, n-point quadrature formulas for a measure [mu] on S = (a, b), 0 [less than or equal to] a < b [less than or equal to] [infinity], exactly integrating Laurent polynomials belonging to certain subspaces of dimension 2n are characterized. A three-term recurrence relation for the sequence of orthogonal Laurent polynomials and a Christoffel-Darboux formula are also mentioned. The same problem is studied in Section 3 for the unit circle T = {z [member of] C : [absolute value of z] = 1} and unlike the situation in Section 2, it is proved that there can not exist an n-point quadrature formula with nodes on T that is exact in an appropriate subspace of [DELTA] of dimension 2n. The corresponding formulas with highest degree of exactness in this situation (2n - 1) are characterized and a three-term recurrence relation for the sequence of monic orthogonal Laurent polynomials on the unit circle is proved. In Section 4 we establish that the monic orthogonal Laurent polynomials are closely related to the first and second singular predictor polynomials computed recursively for the split Levinson algorithm. Finally, in Section 5 illustrative numerical examples are given.

2. The Real Half-Line: L-Gaussian Quadratures. In this section we will assume that the measure p is supported on the interval (a, b) where 0 [less than or equal to] a < b < +[infinity]. Thus, we are interested in approximating

[I.sub.[mu]](f) = [[integral].sup.b.sub.a] f(x)d[mu](x) (2.1)

by means of an n-point quadrature rule

[I.sub.n] (f) = [n.summation over (j=1)] [A.sub.j] f ([x.sub.j]) (2.2)

where the nodes (which are assumed distinct and in (a, b)) and weights are to be determined by imposing

[I.sub.[mu]] (R) = [I.sub.n] (R) (2.3)

for any Laurent polynomial R with L-degree as large as possible, i.e., R [member of] [L.sub.N], where N = N(n) and as large as possible. For this purpose, proceed in a way similar to the polynomial case. Thus, since [L.sub.n-1], is a Tchebyshev space of dimension n on any subinterval of (a, b), given n distinct nodes [x.sub.1], [x.sub.2], ..., [x.sub.n] on (a, b), the weights [A.sub.1], ..., [A.sub.n] can be determined so that

[I.sub.[mu]] (R) = [I.sub.n] (R) = [n.summation over (j=1)] [A.sub.j] f([x.sub.j]), [for all]R [member of] [L.sub.n-1]. (2.4)

Furthermore, it also holds that if [L.sub.n-1](f, x) represents the unique Laurent polynomial in [L.sub.n-1] interpolating f at the nodes [{[x.sub.j]}.sup.n.sub.j=1], then

[I.sub.n] (f) = [I.sub.[mu]] ([L.sub.n-1](f,x)). (2.5)

For this reason [I.sub.n] (f) given by (2.4) sometimes will be called of interpolatory type.

Now, in order to increase the "L-degree" of exactness we have the following result.

THEOREM 2.1. Let [I.sub.n] (f) = [[summation].sup.n.sub.j=1] [A.sub.j] f ([x.sub.j]) be an n-point quadrature formula such that [x.sub.i] [not equal to] 0 for i = 1, ..., n. Then, [I.sub.n] (f) = [I.sub.[mu]] (f) for all f [member of] [L.sub.n+r] with r [greater than or equal to] 0 if and only if

1. [I.sub.n] (f) is of interpolatory type

2.

(2.6) < [R.sub.n] (x), h (x) >= 0, [for all]h(x) [member of] [L.sub.r] (i.e., [R.sub.n](x) [perpendicular to] [L.sub.r]),

where [R.sub.n](x) = [Q.sub.n](x)/[x.sup.p(n)] = [[PI].sup.n.sub.j=1](x-[x.sub.j]) / [x.sup.p(n)] [member of] [L.sub.n].

Proof. Similar to the proof of Theorem 3.1 (in the next section). See also [20] for the polynomial situation with r = n - 1.

It should be taken into account that since we are dealing with real-valued functions, we do not need complex conjugates in the inner product. On the other hand and concerning the integer r, one actually has 0 [less than or equal to] r [less than or equal to] n - 1, since it is very easy to check that there does not exist an n-point quadrature formula which is exact in [L.sub.N] with N [greater than or equal to] 2n. Thus, we are interested in constructing quadrature formulas, like (2.2), exact in [L.sub.n+r] with 0 [less than or equal to] r [less than or equal to] n - 1. From (2.6) it follows that

[R.sub.n](z) = [n.summation over (j=r+1)] [[alpha].sub.j][[psi].sub.j](z), (2.7)

where [{[[psi].sub.n](z)}.sup.[infinity].sub.n=0] is the corresponding sequence of orthonormal Laurent polynomials for the measure [mu] on (a, b). Now, the following question immediately arises: Is it possible to choose the parameters [[alpha]sub.j] in (2.7) so that [R.sub.n](z) has exactly n distinct zeros in (a, b) ?. A positive response can be given when the largest domain of exactness is required, i.e. r = n-1.

THEOREM 2.2. Let [[psi].sub.n] (x) denote the n-th orthonormal Laurent polynomial with respect to the measure [mu] on (a, b). Then,

1. [[psi].sub.n] (x) has exactly n distinct zeros on (a, b).

2. The zeros of [[psi].sub.n] (x) and [[psi].sub.n+1] (x) interlace.

3. Let [x.sub.1], ..., [x.sub.n] be the zeros of [[psi].sub.n] (x). Then there exist positive weights [A.sub.1], ..., [A.sub.n], such that

[I.sub.n] (f) = [n.summation over (j=1)] [A.sub.j]f ([x.sub.j]) = [I.sub.[mu]] (f) = [[integral].sup.b.sub.a] f(x)d[mu](x), [for all]f [member of] [L.sub.2n-1]. (2.8)

Proof. See [7], [19] or [16].

REMARK 2.3. The quadrature formulas given by (2.8) were earlier introduced in [19] (see also [18]). Through this last paper a new area of mathematics was opened leading to a strong moment problem and related topics (continued fractions, two-point Pade approximants, ...) and where orthogonal Laurent polynomials made their appearance. These quadrature rules which we shall refer to as L-Gaussian formulas, have been considered by different authors in the last years (see [7], [4] or [5]). For alternative approaches see [26] and [22]. As for numerical experiments see [26], [2] and [15].

In the rest of the section we will survey the most relevant features of the L-Gaussian formulas in connection with the orthogonal Laurent polynomials. The closest parallelism with classical Gaussian formulas and orthogonal polynomials will be specially emphasized. We start with a three-term recurrence relation (see [26]):

THEOREM 2.4. Let [{[R.sub.[n](z)]}.sup.[infinity].sub.n=0] be the sequence of orthogonal Laurent polynomials normalized as follows. Set [R.sub.n] (z) = [B.sub.n](z) / [z.sup.p(n)] where [B.sub.n] (z) is a monk polynomial of degree n. Then

[R.sub.n+1] (z) = (z - [[beta].sub.n+1])[z.sup.d(n)] [R.sub.n] (z) - [[alpha].sub.n+1] [R.sub.n-1](z) , n [greater than or equal to] 0, (2.9)

with [R.sub.-1] (z) [equivalent to] 0, [R.sub.0] (z) [equivalent to] 1, where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[[alpha].sub.1] = [[mu].sub.0], [[beta].sub.1] = [[mu].sub.0] / [[mu].sub.-1], [[alpha].sub.n+1] = [[rho].sub.n] / [[rho].sub.n-1], [[beta].sub.n+1] = -[[alpha].sub.n+1] [[sigma].sub.n-1] / [[sigma].sub.n], n [greater than or equal to] 1.

Here,

[[mu].sub.k] = [[integral].sup.b.sub.a] [x.sup.k]d[mu](x), [[rho].sub.n] =< [R.sub.n], [t.sup.P(n)] >, [[sigma].sub.n] =< [R.sub.n], [t.sup.-(q(n)+1)] > .

(Recall p(n) = E [n+1/2], q(n) = E [n/2]). Furthermore, both [[alpha].sub.n] and [[beta].sub.n] are positive for n [greater than or equal to] 1.

REMARK 2.5. In [7] a Favard-type theorem for the sequence [{[R.sub.n](z)}.sup.[infinity].sub.n=0] is proved.

For the weights [A.sub.1], ..., [A.sub.n] in the n-point L-Gaussian formula, similar results to the polynomial case are also true. Since [[psi].sub.n] [member of] [L.sub.n] = [[DELTA].sub.-p(n),q(n)] and has n distinct zeros on (0, [infinity]) we can write [[psi].sub.n] (x) = [v.sub.n][x.sup.-p(n)] + ... + [u.sub.n][x.sup.q(n)] ([u.sub.n][v.sub.n] [not equal to] 0). We will also require the normalization conditions: [u.sub.n] > 0 for all n = 0, 1, 2, .... First, we need the following Christoffel-Darboux formula (see [24]):

THEOREM 2.6 (Christoffel-Darboux). Let [{[[psi].sub.n](x)}.sup.[infinity].sub.n=0] be the orthonormal sequence of Laurent polynomials. Then, setting [lambda](n) = [(-1).sup.n], one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.10)

Let [l.sub.j] [member of] [L.sub.n-1] be such that [l.sub.j] ([x.sub.k]) = [[delta].sub.j,k], where {[x.sub.1], ..., [x.sub.n]} as usual are the zeros of the n-th orthonormal Laurent polynomial [[psi].sub.n] (x). Clearly [l.sup.2.sub.j] [member of] [L.sub.2n-1], hence [A.sub.j] = [I.sub.[mu]] ([l.sup.2.sub.j]) > 0, j = 1, ..., n. On the other hand, from (2.10) with x = [x.sub.j] and letting y tend to [x.sub.j], the following theorem can be proved (see [6] for details):

THEOREM 2.7. Let [[psi].sub.n] (x) be the n-th orthonormal Laurent polynomial for the measure [mu]. Then, the weights [{[A.sub.j]}.sup.n.sub.j=1] are given by

[A.sub.j] = 1 / [[summation].sup.n-1.sub.i=0][[psi].sup.2.sub.i]([x.sub.j]), j=1, ..., n.

REMARK 2.8. As a conclusion we can say that the n-point L-Gaussian quadrature rule (1.1) is completely characterized in terms of the orthonormal Laurent polynomials [{[[psi].sub.k](z)}.sup.n.sub.k=0] which can be recursively computed by relation (2.9).

3. The Unit Circle: Szego Quadratures. In this section we shall be concerned with positive measures supported on the unit circle T which in the sequel will be denoted by [sigma] (in order to avoid a possible confusion with the measure [mu] on (0, [sigma]) considered in Section 2). Since we will deal now with complex-valued functions, complex conjugation will be again required in the inner product induced by [sigma], i.e.

< f, g >[sub.[sigma]] = [[integral].sup.[pi].sub.-[pi]] f ([e.sup.i[theta]]) [bar.g([e.sup.i[theta])]d[sigma]([theta]), f, g [member of] [L.sup.[sigma].sub.2] (T) (3.1)

We are interested in approximating the integral

[I.sub.[sigma]](f) = [[integral].sub.T] f(z) d[sigma](z) = [[integral].sup.[pi].sub.-[pi]] f([e.sup.i[theta]])d[sigma]([theta]) (3.2)

by an n-point quadrature rule with nodes on T,

[I.sub.n] (f) = [n.summation over (j=1)] [[lambda].sub.j]f([z.sub.j]), [z.sub.j] [not equal to] [z.sub.k] (j [not equal to] k), [z.sub.j] [member of] T, [for all]j = 1, ..., n (3.3)

Thus, starting from n distinct nodes [z.sub.1], ..., [z.sub.n] on T, and since [L.sub.n-1] is also a Tchebyshev space of dimension n on T, weights [[lambda].sub.1] , ..., [[lambda].sub.n] can be determined so that

[I.sub.[sigma]](f) = [I.sub.n](f), [for all]f [member of] [L.sub.n-1]. (3.4)

Here it also holds that [I.sub.n] (f), characterized by (3.4), can be represented as

[I.sub.n] (f) = [I.sub.[sigma]] ([R.sub.n-1](f,x)), (3.5)

where [R.sub.n-1](f, x) is the Laurent polynomial in [L.sub.n-1] which interpolates f at the nodes [{[z.sub.j]}.sup.n.sub.j=1]. However, the first difference with the real half-line appears when an appropriate selection of the nodes [{[z.sub.j]}.sup.n.sub.j=1] is required in order to increase the dimension of the subspace of Laurent polynomials where exactness of the quadrature formula takes place. Because of the inner product (3.1), we need the following subspaces for r [greater than or equal to] 0,

[L.sub.r*] = {f [member of] [DELTA], [f.sub.*] [member of] [L.sub.r]} = [[DELTA].sub.-q(r),p(r)],

where for f [member of] [DELTA]: [f.sub.*] (z) = [bar.f(1/z)] ("substar-conjugation"). Under these conditions, orthogonality immediately arises, as follows (compare with Theorem 2.1).

THEOREM 3.1. Let [I.sub.n] (f) = [[summation].sup.n.sub.j=1] [lambda]j f ([z.sub.j]) be an n-point quadrature formula such that [z.sub.j] [not equal to] 0 for j = 1, ..., n. Then [I.sub.n] (f) is exact in [L.sub.n] [L.sub.r*], r [greater than or equal to] 0, if and only if

1. [I.sub.n] (f) is exact in [L.sub.n-1].

2.

< [R.sub.n](z), g(z) > [sub.[sigma]] = 0, [for all]g [member of] [L.sub.r] (i.e., [R.sub.n] [perpendicular to] [L.sub.r]), (3.6)

where [R.sub.n](z) = [Q.sub.n](z) / [z.sup.p(n)] = [[PI].sup.n.sub.j=1](z-[z.sub.j]) / [z.sup.p(n)] [member of] [L.sub.n].

Proof.- "[right arrow]" 1.- Trivial. 2.- Recall that [L.sub.n][L.sub.r*] = [[DELTA].sub.-(p(n)+q(r)),(q(n)+p(r))] and [L.sub.r] = [[DELTA].sub.-p(r),q(r)]. Thus, < [R.sub.n](z),g(z) > [sub.[sigma]] = 0 for all g(z) [member of] [L.sub.r] is equivalent to < [R.sub.n](z), [z.sup.j] > [sub.[sigma]] = 0 for all -p (r) [less than or equal to] j [less than or equal to] q (r). Now,

< [R.sub.n](z), [z.sup.j] > [sub.[sigma]] = [[integral].sup.[pi].sub.-[pi]] [R.sub.n](z)[z.sup.-j] [sigma] ([theta]) d[theta] = [[integral].sup.[pi].sub.-[pi]] [Q.sub.n] (z) / [z.sup.p(n)+j] [sigma]([theta])d[theta].

For j = q (r) and since [Q.sub.n](z) has exact degree n, we see that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

On the other hand, for j = -p(r) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Thus, for -p (r) [less than or equal to] j [less than or equal to] q (r), [R.sub.n](z)[z.sup.-j] [member of] [L.sub.n][L.sub.r*]. Therefore one can write

< [R.sub.n](z), [z.sup.j] > [sub.[sigma]] = [I.sub.[sigma]] ([R.sub.n](z)[z.sup.-j]) = [I.sub.n] ([R.sub.n](z) [z.sup.-j]) = 0

since [R.sub.n] ([z.sub.j]) = 0.

"[left arrow]" Let [z.sub.1], ..., [z.sub.n] be the zeros of [R.sub.n] (z) [member of] [L.sub.n] so that (3.6) holds. Then a quadrature formula based upon these nodes, such as [I.sub.n] (f) = [[summation].sub.n.sub.j=1] [[lambda].sub.j] f ([z.sub.j]), can be determined so that it is exact in [L.sub.n-1]. Setting [L.sup.(j)] (z) [member of] [L.sub.n-1] such that [L.sup.(j)] ([z.sub.k]) = [[delta].sub.jk], then [[lambda].sub.j] = [I.sub.[sigma]] ([L.sup.(j)]) for all j = 1, ..., n.

Take L (z) [member of] [L.sub.n][L.sub.r*] and define A (z) = L (z) - [[summation].sup.n.sub.j=1] L([z.sub.j])[L.sup.(j)] (z). Thus, A [member of] [L.sub.n][L.sub.r*], A([z.sub.j]) = 0 for all j = 1, ..., n, and we can write

A(z) = [Q.sub.n](z)S(z) / [z.sup.p(n)+q(r)], S (z) [member of] [[PI].sub.r].

Observe that A (z) = [R.sub.n](z)H(z) with H (z) = S(z)/[z.sup.q(r)]. Hence, setting T(z) = [H.sub.*](z) [member of] [L.sub.r], we have

[I.sub.[sigma]] (A(z)) = [I.sub.[sigma]] ([R.sub.n] (z) [T.sub.*] (z)) =< [R.sub.n] (z), T (z) > [sub.[sigma]] = 0,

which yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Consider now the family of monic Szego polynomials [{[[rho].sub.n](z)}.sup.[infinity].sub.n=0] orthogonal with respect to the inner product (3.1), i.e., for n = 0,1, ... then [[rho].sub.n] (z) satisfies

< [[rho].sub.n] (z), [z.sup.m] > [sub.[sigma]] = 0 (0 [less than or equal to] m [less than or equal to] n - 1) (3.7)

and

< [[rho].sub.n](z), [[rho].sub.n](z) > [sub.[sigma]] = < [[rho].sub.n](z), [z.sup.n] > [sub.[sigma]] = [[DELTA].sub.n] = [[DELTA].sub.n-1], (3.8)

where [[DELTA].sub.n] is the n-th Toeplitz determinant of the trigonometric moments of [sigma]. Set [[rho].sup.*.sub.n](z) = [z.sup.n][[rho].sub.n*] (z) (reciprocal polynomial). Then, from the above orthogonality conditions both for [{[[rho].sub.n](z)}.sup.[infinity].sub.n=0] and [{[[phi].sub.n](z)}.sup.[infinity].sub.n=0] (recall that [[phi].sub.n](z) denotes the n-th orthogonal Laurentpolynomial with respect to the measure [sigma]) the following relation between [[phi].sub.n](z) and [[rho].sub.n](z) is deduced:

[[phi].sub.2n](z) = 1/[z.sup.n] [[rho].sub.2n](z), [[phi].sub.2n+1] (z) = 1 / [z.sup.z+1] [[rho].sup.*.sub.2n+1](z). (3.9)

On the other hand, it is well known (see [1]) that the zeros of the polynomial [[rho].sub.n](z) are located in D = {z [member of] C : [absolute value of z] < 1} and consequently, the zeros of [[rho].sup.*.sub.n](z) in E = {z [member of] C: [absolute value of] > 1}. From this fact and (3.9) we obtain the following result

LEMMA 3.2. The zeros of [[phi].sub.2n] (z) and [[phi].sub.2n+1] (z) are located in D and E respectively.

REMARK 3.3. Actually r in Theorem 3.1 should be taken such that 0 [less than or equal to] r [less than or equal to] n - 1. Indeed, assume r [greater than or equal to] n, i.e., r = n + k with k [greater than or equal to] 0. Then, [L.sub.n] [L.sub.r*] = [[DELTA].sub.-(n+q(k)),(n+P(k))] since for p(n) = E [n+1/2] and q(n) = n -p(n) = E [n/2] it holds that p(r + s) = p(r) + p(s) and q(r + s) = q(r) + q(s). Thus, if [Q.sub.n](z) = [[PI].sup.n.sub.j=1] (z - [z.sub.j]), then

[[absolute value of [Q.sub.n](z)].sup.2] = [Q.sub.n](z)[bar.[Q.sub.n](z)] [member of] [[DELTA].sub.-n,n] [subset] [L.sub.n][L.sub.r*] (z = [e.sup.i[theta]]).

Hence, 0 < [[integral].sup.[pi].sub.-[pi]] [[absolute value of [Q.sub.n](z)].sup.2]d[sigma]([theta]) and [I.sub.n]([[absolute value of [Q.sub.n](z)].sup.2]) = 0.

Let us next see what happens with the largest reachable domain of exactness, i.e., if r = n - 1. Unlike the situation on the half-line we have the following negative result.

THEOREM 3.4. There can not exist an n-point quadrature formula with nodes on T to be exact in [L.sub.n] [L.sub.(n-1)*]. Proof.- Suppose that there exists such a formula:

In (f) = [n.summation over (j=1)][[lambda].sub.j]f([z.sub.j]) = [I.sub.[sigma]](f)

where [absolute value of [z.sub.i]] = 1 for i = 1, ..., n. Define [R.sub.n] (z) = (z-[z.sub.1])...(z-[z.sub.n]) / [z.sup.p(n)] [member of] [L.sub.n]. Then by Theorem 3.1, [R.sub.n] [perpendicular to] [L.sub.n-1].

Thus, if [{[[phi].sub.k](z)}.sup.[infinity].sub.k=0] is a sequence of orthogonal Laurent polynomials, then [R.sub.n](z) = [[lambda].sub.n][[phi].sub.n](z)([[lambda].sub.n] [not equal to] 0). So, the nodes [z.sub.1], ..., [z.sub.n] are the zeros of [[phi].sub.n] (z), which cannot be located on T by Lemma 3.2 and a contradiction arises.

Now, the next step would be to consider the case r = n - 2. More precisely, one might wonder is it possible to construct an n-point quadrature formula with distinct nodes {[z.sub.j]} on T to be exact in [L.sub.n][L.sub.(n-2)*]? By Theorem 3.1 if [R.sub.n](z) = [[PI].sup.n.sub.j=1] (z-[z.sub.j])/[z.sup.p(n)] then [R.sub.n](z) [perpendicular to] [L.sub.n-2]. Hence,

[R.sub.n] (z) = [alpha][[phi].sub.n] (z) + [beta][[phi].sub.n-1] (z), [alpha], [beta] [member of] C, (3.10)

where [{[[phi].sub.k](z)}.sup.[infinity].sub.k=0] o is a sequence of orthogonal Laurent polynomials. So, from (3.10) the above question can be reformulated as follows: Is it possible to find [alpha] and [beta] in (3.10) so that [R.sub.n] (z) has its zeros on T?

THEOREM 3.5. Parameters [alpha] and [beta] in (3.10) can be conveniently chosen so that [R.sub.n] (z) = [alpha][[phi].sub.n] (z) + [beta][[phi].sub.n-1] (z) has exactly n distinct zeros on T.

Proof. Setting [R.sub.n](z) = [N.sub.n](z)/[z.sup.p(n)], with [N.sub.n](z) being a polynomial of degree n at most depending on [alpha] and [beta]. These parameters should be chosen so that [N.sub.n](z) satisfies the following:

1. [N.sub.n] (z) has exact degree n.

2. [N.sup.*.sub.n] (z) = [[lambda].sub.n] [N.sub.n] (z) with [[lambda].sub.n] [not equal to] 0.

3. < [N.sub.n](z), [z.sup.k] > [sub.[sigma]] = 0, 1 [less than or equal to] k [less than or equal to] n - 1, < [N.sub.n](z), 1 > [sigma] [not equal to] 0, < [N.sub.n](z), [z.sup.n] > [sub.[sigma]] [not equal to] 0.

Now by using Theorem 6.2 in [17] the proof follows.

Hence, from Theorems 3.1 and 3.5 we obtain the the following result

THEOREM 3.6. Let {[z.sub.1] ..., [z.sub.n]} be the n distinct zeros of [R.sub.n](z) = [alpha][[phi].sub.n] (z) [beta][[phi].sub.n-1](z) as given in Theorem 3.5. Then, there exist positive numbers [[lambda].sub.1], ..., [[lambda].sub.n] such that

[I.sub.n] (f) = [n.summation over (j=1)](zj) = [I.sub.[sigma]](f), [for all]f [member of] [L.sub.n][L.sub.(n-2)*].

Proof. It only remains to prove that the weights [[lambda].sub.1] , ..., [[lambda].sub.n] are positive. Indeed, for 1 [less than or equal to] j [less than or equal to] n, let [L.sup.(j)](z) [member of] [L.sub.(n-1)] be such that [L.sub.j]([z.sub.k]) = [[delta].sub.j,k] for 1 [less than or equal to] j, k [less than or equal to] n. Then, for z [member of] T, [[absolute value of [L.sub.j] (z)].sup.2] [member of] [L.sub.(n-1)][L.sub.(n-1)*] = [L.sub.n] [L.sub.(n-2)*]. as can be checked easily. Now, the proof immediately follows, since

0 < [I.sub.[sigma]]([[absolute value of [L.sub.j]].sup.2]) = [n.summation over (k=1)] [[lambda].sub.k] [absolute value of [L.sub.j]([z.sub.k])].sup.2] = [[lambda].sub.j], j=1, .... n.

So, from Theorems 3.1 and 3.6 quadrature formulas with nodes on T which are exact in [L.sub.n][L.sub.(n-2)*] can be essentially characterized by means of the Laurent polynomials of the form

[R.sub.n] (z) = [alpha][[phi].sub.n] (z) + [beta][[phi].sub.n-1](z), (3.11)

which by analogy with the polynomial situation (see e.g. [1]) could be called "para-orthogonal Laurent polynomials". On the other hand, one also has that

[L.sub.n][L.sub.(n-2)*] = [L.sub.(n-1)][L.sub.(n-1)*] = [[DELTA].sub.-(n-1),(n-1)].

Quadrature formulas with nodes on T to be exact in [[DELTA].sub.-(n-1),(n-1)] were earlier introduced by Jones, et al. ([17]) in connection with the solution of the trigonometric moment problem (see also [29] and [14]). Such quadratures were named "Szego formulas" and all studied in a series of recent papers: [13], [25], [8] [9].

As a consequence, making [R.sub.n] (z) = [alpha][[phi].sub.n] (z) + [beta][[phi].sub.n-1](z) = [N.sub.n](z)/ [z.sup.p(n)] by Theorem 6.1 in [17] one can write

[N.sub.n] (z) = [[lambda].sub.n] [[[rho].sub.n] (z) + [[tau].sub.n][[rho].sup.*.sub.n](z)], [[lambda].sub.n] [not equal to] 0, [absolute value of [[tau].sub.n]] = 1.

In [17] polynomials of the form [[rho].sub.n] (z) + [tau][[rho].sup.*.sub.n] (z) are called "para-orthogonal". In the next section we will concentrate on the particular cases [[tau].sub.n] = [+ or -]1.

In the rest of this section we will continue to emphasize the analogies and differences between the real half-line and the unit circle. The following theorem yields an alternative expression for the weights [[lambda].sub.j] (compare with Theorem 2.7).

THEOREM 3.7. Let [{[[??].sub.k]}.sup.[infinity].sub.k=0] be the sequence of orthonormal Laurent polynomials and set [I.sub.n] (f) = [[summation].sup.n.sub.j=1] [[lambda].sub.j] f ([z.sub.j]) = [I.sub.[sigma]] (f) for all f [member of] [L.sub.n] [L.sub.(n-2)*]. Then

[[lambda].sub.j] = 1 / [[summation].sup.n-1.sub.k=0] [[absolute value of [[??].sub.k]([z.sub.j])].sup.2].

Proof. Direct consequence of Proposition 3.4 in [13] and (3.9).

REMARK 3.8. Observe that unlike the real half-line situation now the nodes [{[z.sub.j]}.sup.n.sub.j=1] in the quadrature formula are not the zeros of [[??].sub.n](z)

THEOREM 3.9 (Three-term recurrence relation). Let [{[[delta].sub.n]}.sup.[infinity].sub.n=1] be the sequence of Schur parameters (or reflection coefficients) for the measure [sigma],, i.e., [[delta].sub.n] = [[rho].sub.n](0), where [[rho].sub.n](z) is the n-th monic Szego polynomial with respect to the measure [sigma], and let [{[[phi].sub.n](z)}.sup.[infinity].sub.n=0] 0 be the sequence of monic orthogonal Laurent polynomials. Then, it holds that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.12)

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof.- From the recurrence relations satisfied by the Szego polynomials (see [28]):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.13)

we deduce the initial conditions [[phi].sub.0](z) = [[rho].sub.0](z) [equivalent to] 1, [[phi].sub.1](z) = 1/z [[rho].sup.*.sub.1](z) = 1/z (1 + [[bar.[delta]].sub.1]z) = [[bar.[delta]].sub.1] + 1/z (since [[rho].sub.1] z) = [[rho].sub.1](0) + z = [[delta].sub.1] + z) and the relations

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.14)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.15)

Now (3.12) follows from (3.14) and (3.15).

REMARK 3.10. For an alternative proof of Theorem 3.9 based upon certain continued fractions, see [30]. On the other hand, if the trigonometric moments [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are real, then [[delta].sub.k] [member of] R and (3.12) can be written as

[[phi].sub.n](z) = ([[delta].sub.n] + [[delta].sub.n-1][z.sup.(-1)n)] [[phi].sub.n-1](Z) + (1 - [[absolute value of [[delta].sub.n-1]].sup.2])[z.sup.(-1)n] [[phi].sub.n-2](Z), n [greater than or equal to] 2. (3.16)

EXAMPLE 3.11. Consider the biparametric family of measures

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In this case it is known (see [28]) that

[[delta].sub.2n] = [alpha] + [beta] + 1/2n + [alpha] + [beta] + 1 [member of] R, [[delta].sub.2n+1] = [alpha] - [beta] / 2n + [alpha] + [beta] + 2 [member of] R.

When [beta] = -1/2, the measure and the Schur parameters are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and in the particular case [alpha] = [beta] = -1/2, we have d[sigma]([theta]) = d[theta] (Lebesgue measure) and [[delta].sub.n] = 0 [for all]n = 1, 2,.... Now from Theorem 3.9, we deduce

[[phi].sub.n](Z) = [Z.sup.(-1)n][[phi].sub.n-2](Z)), [[phi].sub.0](Z) [equivalent to] 1, [[phi].sub.1](Z) [equivalent to] 1/z.

So, [[phi].sub.2n](Z) = [z.sup.n] and [[phi].sub.2n+1](Z) = 1/[z.sup.n+1] for all n = 0, 1, 2,.... Since

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

the orthonormal sequence [{[[??].sub.n]}.sup.[infinity].sub.n=0] is given by

[[??].sub.2n](z) = [z.sup.n]/[square root of 2[pi]], [[??].sub.2n+1](z) = 1/[square root of 2[pi]][z.sup.n+1], n = 0,1,2, ...

and from (3.11) we obtain:

1. If n is even:

[R.sub.n](z) = [alpha][[phi].sub.n](z) + [beta][[phi].sub.n-1](z) = [alpha][z.sup.n] + [beta]/[z.sup.n/2] = 0 [??] [alpha][z.sup.n] + [beta] = 0, a [not equal to] 0, [absolute value of [beta]/[alpha]] = 1.

2. If n is odd:

[R.sub.n](z) = [alpha][[phi].sub.n](z) + [beta][[phi].sub.n-1](z) = [beta][z.sup.n] + [alpha]/[z.sup.n+1] = [beta] [not equal to] 0, [absolute value of [alpha]/[beta]] = 1.

Thus, the nodes [{[z.sub.j]}.sup.n.sub.j=1] of any n-point Szego formula for the Lebesgue measure are the roots of [z.sup.n] + [tau] = 0 , [absolute value of [tau]] = 1 and the weights [{[[lambda].sub.j]}.sup.n.sub.j=1], are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(compare with [20, pp. 73-74]).

4. A connection with the split Levinson algorithm. In general, when one needs to compute an n-point Szego formula for a measure [sigma] on T, we start from the usual known information about [sigma], i.e., its trigonometric moments [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], k [member of] Z, so that, as already seen through Section 3, we essentially need to calculate the (monic) Szego polynomials [[rho].sub.0], ..., [[rho].sub.n]. As a rule, these polynomials are not explicitly known (unless for some particular measures [sigma]). Hence they should be computed from the relations (3.13). Thus, the well known Levinson algorithm arises [21]. Assume now (as will be done in the rest of the section) that the trigonometric moments [[mu].sub.k] are real (a very common situation as illustrated in Section 5). In this case, Levinson's algorithm is redundant, in the sense that more operations than the required ones are carried out. To overcome this drawback a modification of this algorithm was studied giving rise to the so-called "split Levinson algorithm" (see [10]). Here, certain polynomials closely related to Szego polynomials are computed so that the number of required operations is halved. Let us next see how these polynomials arise in the context of the construction of Szego quadrature formulas. Indeed, in order to compute the nodes of an n-point Szego formula and since we are dealing with real moments, it seems natural to handle polynomials (whose zeros provide us with the required nodes) with real coefficients. In other words, the parameters [alpha] and [beta] in (3.11) will be chosen so that

[R.sub.n](z) = [alpha][[phi].sub.n](z) + [beta][[phi].sub.n-1](z) = [N.sub.n](z)/[z.sup.P(n)] = [[lambda].sub.n]([[rho].sub.n](z) + [[tau].sub.n][[rho].sup.*.sub.n](z))/[z.sup.P(n)],

where [[lambda].sub.n] [not equal to] 0, [[tau].sub.n] [member of] R and [absolute value of [tau].sub.n] = 1.

Set first [[tau].sub.n] = 1. Then by (3.13) it follows (now [[delta].sub.n] [member of] R)

[N.sub.n](z) = [[lambda].sub.n] (1 + [[delta].sub.n]) [[[rho].sup.*.sub.n-1](z) + z[[rho].sub.n-1](z)].

Thus, taking [[lambda].sub.n] = 1/1+[[delta].sub.n], one has

[N.sub.n](z) = [[rho].sup.*.sub.n-1](z) + Z[[rho].sub.n-1](z) = [P.sub.n](z). (4.1)

By [10] it can be shown that the polynomial [P.sub.n](z) represents the basic ingredient in the "symmetric version" of the Levinson-split algorithm (here it should be noted that from the context of the solution of Toeplitz systems, [[rho].sup.*.sub.n](z) is called the "predictor polynomial"). For this reason, [P.sub.n](z) given by (4.1) is sometimes called the "first singular predictor polynomial of degree n", since roughly speaking it can be interpreted as though in (3.13), [[delta].sub.n] is forced to take the value 1.

Setting [P.sub.n](z) = [[summation].sup.n.sub.j=0][P.sub.n,j][z.sup.j] in [10] it is proved that the sequence {[P.sub.n](z)} satisfies the following three-term recurrence relation,

[P.sub.n+1](z) - (1 + z)[P.sub.n](z) + [[alpha].sub.n]z[P.sub.n-1](z) = 0, n [greater than or equal to] 1, (4.2)

with initial conditions [P.sub.0](z) [equivalent to] 2, [P.sub.1](z) = 1 + z, and where [[alpha].sub.n] = [[gamma].sub.n]/[[gamma].sub.n-1] with [[gamma].sub.n] = [[summation].sup.n.sub.i=0][[mu].sub.i][P.sub.n,i] for n [greater than or equal to] 1 and [[gamma].sub.0] = [[mu].sub.0].

On the other hand, for [T.sub.n] = -1, again by (3.13) it follows that

[N.sub.n](z) = [[lambda].sub.n]([[delta].sub.n] - 1) [[[rho].sup.*.sub.n-1](z) - z[[rho].sub.n-1](z)].

Thus, with [[lambda].sub.n] = 1/[[delta].sub.n-1], one can write for n [greater than or equal to] 1,

[N.sub.n](z) = [[[rho].sup.*.sub.n-1](z) - z[[rho].sub.n-1](z) = [[??].sub.n](z). (4.3)

Now, [[??].sub.n](z) is the basis of the "antisymmetric version" of the Levinson-split algorithm and used to be called the "second singular predictor polynomial" (set [[delta].sub.n] = -1 in (3.13)). By defining [[??].sub.0](z) [equivalent to] 1, the polynomials [P.sub.n](z) for n = 1, 2, ... are given by [[??].sub.1](z) = 1 - z, [[??].sub.2](z) = 1 - [z.sup.2] and for n [greater than or equal to] 2 it holds

[[??].sub.n+1](z) - (1 + z)[[??].sub.n](z) + [[??].sub.n]Z[[??].sub.n-1](z) = 0, (4.4)

where [[??].sub.n] = [[??].sub.n]/[[??].sub.n-1] with [[??].sub.n] = [[summation].sup.n.sub.i=0][[mu].sub.i][[??].sub.n,I] and [[??].sub.n](z) = [[summation].sup.n.sub.j=0][[??].sub.n,j][z.sup.j].

In order to illustrate relations (4.2) and (4.4), let us consider the Lebesgue measure, i.e., d[sigma]([theta]) = d[theta]. Then,

[[mu].sub.k] = [integral].sup.[pi].sub.[pi] [e.sup.-ik[theta]] d[theta] = 0, k [greater than or equal to] 1, [[mu].sub.0] = 2[pi].

From (4.2) it can be easily checked that [P.sub.k](0) = [p.sub.k,0] = 1 for k [greater than or equal to] 1, which gives [[gamma].sub.k] = [[mu].sub.0] for k [greater than or equal to] 0, and hence [[alpha].sub.k] = 1 for k [greater than or equal to] 1. Then (4.2) becomes

[P.sub.n+1](z) - (1 + z)[P.sub.n](z) + Z[P.sub.n-1](z) = 0, n [greater than or equal to] 1, (4.5)

with [P.sub.0](z) [equivalent to] 2 and [P.sub.1](z) [equivalent to] 1. Since (4.5) is a linear second order finite diference equation with constant coefficients, it holds that

[P.sub.n](z) = [C.sub.1] + [C.sub.2][z.sup.n], [C.sub.1], [C.sub.2] [member of] R.

From the initial conditions, one has [C.sub.1] = [C.sub.2] = 1, so that an explicit expression for the polynomial [P.sub.n](z) is obtained, namely

[P.sub.n](Z) = [z.sup.n] + 1, n = 0, 1, 2, ... (4.6)

Similarly for the sequence {[[??].sub.n](z)} one has

[[??].sub.n+1](z) - (1 + z)[[??].sub.n](z) + z[[??].sub.n-1](z) = 0, n [greater than or equal to] 2, (4.7)

with [[??].sub.1](z) = 1 - z and [[??].sub.2](z) = 1 - [z.sup.2]. Thus, for n [greater than or equal to] 2,

[[??].sub.n](z) = 1 - [z.sup.n]. (4.8)

Clearly, (4.8) holds for n = 1.

When dealing with other measures [sigma], in general, equation (4.5) or (4.7) cannot be explicitly solved so that we need to compute recursively the polynomials [{[P.sub.n](z)}.sup.[infinity].sub.n=0] or [{[[??].sub.n](z)}.sup.[infinity].sub.n=0] (for details concerning stability, see [10]).

We conclude this section rewriting the polynomials [{[P.sub.n](z)}.sup.[infinity].sub.n=0] and [{[[??].sub.n](z)}.sup.[infinity].sub.n=0] in terms of the orthogonal Laurent polynomials [[phi].sub.n](z) and [[phi].sub.n-1](z).

THEOREM 4.1. Let [sigma] be a measure on T with real moments and let [{[[phi].sub.k]}.sup.[infinity].sub.k=0] be the sequence of monic orthogonal Laurent polynomials. Let [P.sub.n](z) and [[??].sub.n](z) be the first and second singular predictor polynomials of degree n respectively. Then

1. [P.sub.n](z) = [Z.sup.E][absolute value of n+1/2] [[[phi].sub.n](z) + (1 - [[delta].sub.n])[[phi].sub.n-1](z)]

2. [[??].sub.n](z) = [Z.sup.E][absolute value of n+1/2] [(-1).sup.n+1][[[phi].sub.n](z) - (1 + [[delta].sub.n])[[phi].sub.n-1](z)]

Proof.

1. Recall that [[phi].sub.2n](z) = 1/[z.sup.n][[rho].sub.2n](z) and [[phi].sub.2n+1](z) = 1/[z.sup.n+1][[rho].sup.*.sub.2n+1](z). Now for [alpha] and [beta] complex numbers, consider

[R.sub.n](z) = [alpha][[phi].sub.n](z) + [beta][[phi].sub.n-1](z).

Then, for n = 2m it follows that

[R.sub.2m](z) = [alpha][[rho].sub.2m](z) + [beta][[rho].sup.*.sub.2m-1](z)/[z.sup.m]

and for n = 2m + 1

[R.sub.2m+1](z) = [alpha][[rho].sup.*.sub.2m+1](z) + [beta]z[[rho].sub.2m](z)/[z.sup.m+1]

Now, making use of (3.13), we have

[R.sub.2m](z) = [alpha]z[[rho].sub.2m-1](z) + ([alpha][[delta].sub.2m] + [beta])[[rho].sup.*.sub.2m-1](z)/[z.sup.m].

Thus, taking [alpha] = 1 and [beta] = 1 - [[delta].sub.2m], it follows that

[R.sub.2m](z) = [P.sub.2m](z)/[z.sup.m],

or equivalently,

[P.sub.2m](z) = [z.sup.m][R.sub.2m](z) = [z.sup.m] [[alpha][[phi].sub.2m](z) + [beta][[phi].sub.2m-1](z)] = = [z.sup.m][[[phi].sub.2m](z) + (1 - [[delta].sub.2m])[[phi].sub.2m-1](z)], m [greater than or equal to] 1. (4.9)

Similarly for n = 2m + 1, taking [alpha] = 1 and [beta] = 1 - [[delta].sub.2m+1], one can write

[P.sub.2m+1](z) = [z.sup.m+1] [[[phi].sub.2m+1](z) + (1 - [[delta].sub.2m+1])[[phi].sub.2m](z)]. (4.10)

So, from (4.9) and (4.10) the proof follows.

2. It can be proved in a similar way

5. Numerical Examples. Suppose that the evaluation of integrals of the form

[[integral].sup.1.sub.-1]f(x)/[(x + [lambda]).sup.r] w(x) dx, r = 1,2,3, ... (5.1)

where f(x) is a continuous function on [-1,1], w(x) [greater than or equal to] 0 a weight function on this interval and [lambda] [member of] R such that [absolute value of [lambda]] > 1, is required. Two methods are to be proposed:

Method 1

Let a and b be any two positive real numbers, such that

b/a = [absolute value of [lambda]] + 1/[absolute value of [lambda]] - 1,

and introduce the change of variable x = h(t) = [absolute value of [lambda]]/[lambda] x (2t - b - a/b - a). Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5.2)

where g(t) = f(h(t)) and [mu](t) = w(h(t)). The integral of the right-hand side in formula (5.2) can be approximated by a Gauss-Laurent quadrature formula.

Method 2

Set z = [e.sup.i[theta]], x = 1/2(z + [z.sup.-1]) = cos[theta] and define the symmetric weight function [sigma]([theta]) on [-[pi], [pi]] by

[sigma]([theta]) = [2.sup.r][[absolute value of [alpha]].sup.r] w(cos[theta])[absolute value of sin[theta]]/[[absolute value of z - [alpha]].sup.r] (5.3)

where [alpha] [not equal to] 0 is the root of the equation [z.sup.2] + 2[lambda]z + 1 = 0 with [absolute value of [alpha]] < 1, that is,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then

[[integral].sup.1.sub.-1]f(x)/[(x + ([lambda] + [lambda]).sup.r]w(x)dx = [[integral].sup.[pi].sub.-[pi]]g([e.sup.i[theta])[sigma]([theta])d[theta] (5.4)

with

g([e.sup.i[theta]) = 1/2f ([e.sup.i[theta]) + [e.sup.-i[theta]/2).

The integral on the right side in formula (5.4) can be approximated by an n-point Szego quadrature formula.

For the numerical experiments we will consider two particular weight functions:

Case 1: w(x) = 1/[square root of 1 - [x.sup.2]] :Tchebyshev weight function of the first kind (x [member of] (-1,1)). The weight function [mu](t) = w(h(t)) is the strong Tchebyshev weight function of the first kind on (a, b) given by

[mu](t) = 1/[square root of b - t][square root of t - a], 0 < a < b < + [infinity]

In this particular situation explicit formulas for the weights and nodes are known (see [26]):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where m = 1,2, ..., p(n) = E [n + 1 / 2],

[alpha] = [([square root of b] - [square root of a]).sup.2], [beta] = [square root of ab], [v.sub.m] = 1 + cos (2m - 1/n [pi])

and

[A.sub.k] = 2[pi]/n [x.sub.k]/[x.sub.k] + [beta], k = 1,2, ..., n, n [greater than or equal to] 1.

On the other hand, the weight function [sigma]([theta]) on [-[pi], [pi]] given by (5.3) is a rational modification of the Lebesgue measure, namely

[sigma]([theta]) = [2.sup.r][[absolute value of [alpha]].sup.r]/[[absolute value of z - [alpha]].sup.r], z = [e.sup.i[theta]],

Choosing the particular values [lambda] = 1.1, 1.01, r = 1, 2, the smooth function f(x) = [e.sup.x] and n = 10 (number of nodes) fixed we will approximate the next four integrals where exact values were computed with MATHEMATICA:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For the Tchebyshev weight function of the first kind w(x), the nodes and weights of the Gaussian quadrature formula have explicit expressions. Because of the presence of a singularity near the interval of integration, we can expect to obtain poor results using this formula, as is established in the next tables:
```Table 1

r = 1 Approximation Exact value

[lambda] = 1.01 9.36117219 ... 10.2639877 ...
[lambda] = 1.1 4.39825773 ... 4.39889820 ...

Table 2

r = 2 Approximation Exact value

[lambda] = 1.01 247.983129 ... 414.487345 ...
[lambda] = 1.1 15.0292310 ... 15.0611749 ...
```

In order to improve the slow convergence in this situation we use the two above mentioned methods:

Results for Method 1

Choosing a = 1 implies b = 21 when [lambda] = 1.1 and b = 201 when [lambda] = 1.01. The nodes and weights obtained are:
```Table 3

[lambda]=1.1 n=10
Node Weights

1.015968076 0.1140210017
1.153792215 0.1263777081
1.490163636 0.1541804078
2.180737182 0.2025926653
3.529658620 0.2733833654
5.949583872 0.3549351653
9.629771149 0.4257258654
14.09241206 0.4741381229
18.20085084 0.5019408226
20.66994082 0.5142975291

Table 4

[lambda]=1.01 n=10

Node Weights

196.7288912 0.5860819973
164.9922240 0.5786005592
113.4050389 0.5584974065
60.84072229 0.5095745959
24.34950355 0.3971049908
8.254788423 0.2312135400
3.303708313 0.1187439348
1.772408016 0.0698211241
1.218239230 0.0497179710
1.021710633 0.0422365330
```

The approximations of the above four integrals are:
```Table 5

r = 1 Approximation Exact value

[lambda] = 1.01 10.26398781 ... 10.26398779 ...
[lambda] = 1.1 4.398898202 ... 4.398898203 ...

Table 6

r = 2 Approximation Exact value

[lambda] = 1.01 414.487344 ... 414.4873459 ...
[lambda] = 1.1 15.06117499 ... 15.06117492 ...
```

Results for Method 2

Expressions for the trigonometric moments for r = 1, 2 are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

By using the split Levinson algorithm computed in FORTRAN with double precision, the following nodes and weights are obtained:
```Table 7

r = 1, [lambda] = 1.1 n = 10

Nodes Weights

0.5231956081 [+ or -] 0.8522126235i 0.4170158485
-0.9718655888 [+ or -] 0.2355361486i 3.8990461817
-0.6941445801 [+ or -] 0.7198356076i 1.5283999088
0.9425118124 [+ or -] 0.3341728347i 0.3334871136
-0.1205684668 [+ or -] 0.992705014i 0.6775681574

Table 8

r = 1, [lambda] = 1.01 n = 10

Nodes Weights

0.5077320439 [+ or -] 0.8615150443i 0.4552583228
-0.1556028488 [+ or -] 0.9878196968i 0.8023120227
0.9406267445 [+ or -] 0.3394426719i 0.3550340683
-0.7396737275 [+ or -] 0.672965658i 2.4403459272
-0.9871949772 [+ or -] 0.1595782654i 18.1061357088

Table 9

r = 1, [lambda] = 1.01 n = 10

Nodes Weights

-0.7826063935 [+ or -] 0.6225168535i 5.7286222677
-0.2523182279 [+ or -] 0.9676443106i 0.9628770791
0.9315406380 [+ or -] 0.3636372365i 0.1801413420
0.4439464841 [+ or -] 0.8960533016i 0.3067180364
-0.9823049310 [+ or -] 0.1872886076i 28.7314933198

Table 10

r = 2, [lambda] = 1.01 n = 10

Nodes Weights

0.9264937183 [+ or -] 0.376310231i 0.2056746266
-0.8669062623 [+ or -] 0.4984711951i 30.7425691490
0.4037996706 [+ or -] 0.9148474332i 0.3833193864
-0.3361546099 [+ or -] 0.9418068158i 1.6432603190
-0.9954580194 [+ or -] 0.0952015315i 1080.4916890184
```

The approximation of the above four integrals are
```Table 11

r = 1 Approximation Exact value

[lambda] = 1.01 10.26398785 .. 10.26398779 ..
[lambda] = 1.1 4.398898196 .. 4.398898203 ..

Table 12

r = 2 Approximation Exact value

[lambda] = 1.01 414.4873471 .. 414.4873459 ..
[lambda] = 1.1 15.06117499 .. 15.06117492 ..
```

From Tables 5-6 and 11-12 one can see that both methods provide us with similar numerical results. In both cases results given by Gaussian formulas are strongly improved.

Case 2: w(x) = [square root of 1 - [x.sup.2]] :Tchebyshev weight function of the second kind (x [member of] [-1,1]). The weight function [mu](t) = w(h(t)) is the strong Tchebyshev weight function of the second kind on [a, b] given by

[mu](t) = [square root of b - t][square root of t - a].

Choosing the particular values [lambda] = 1.1, 1.01, r = 1, the smooth function f(x) = [e.sup.x] and again n = 10 (number of nodes) fixed, we will approximate the next two integrals with exact values computed with MATHEMATICA:

[I.sub.1] = [[integral].sup.1.sub.-1][e.sup.x][square root of 1 - [x.sup.2]]/[absolute value of x + 1.1] dx = 1.67594127475 ...

[I.sub.2] = [[integral].sup.1.sub.-1][e.sup.x][square root of 1 - [x.sup.2]]/[absolute value of x + 1.01] dx = 2.03543204817 ...

For the Tchebyshev weight function of the second kind w(x), the nodes and weights of the Gaussian quadrature formula have also explicit expressions. As in case 1, because of the presence of a singularity near the interval of integration, we can expect to obtain poor results using this formula as displayed in the next table:
```Table 13

r = 1 Approximation Exact value

[lambda] = 1.01 2.421155521 ... 2.03543204817 ...
[lambda] = 1.1 1.608794159 ... 1.67594127475 ...
```

Again, in order to improve the slow convergence in this situation we use the two mentioned methods:

Results for Method 1

From the relation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with

g(x) = f(x)(1 - [x.sup.2]),

the mentioned explicit formulas for the nodes and weights in case 1 can be used. The approximation of the above two integrals are
```Table 14

r = 1 Approximation Exact value

[lambda] = 1.01 2.035432052 ... 2.03543204817 ...
[lambda] = 1.1 1.675941276 ... 1.67594127475 ...
```

Results for Method 2

Expressions for the trigonometric moments are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

when r = 1, and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

when r = 2. By the split Levinson algorithm computed in FORTRAN with double precision, nodes and weights are given by:
```Table 15

r = 1, [lambda] = 1.1 n = 10

Nodes Weights

0.8473567133 + [+ or -] 0.5310241053i 0.0809799064
-0.9095350917 + [+ or -] 0.4156271369i 0.4250829076
-0.5968012792 + [+ or -] 0.8023890784i 0.6749529041
0.4377179772 + [+ or -] 0.8991123247i 0.2923674821
-0.0996095352 + [+ or -] 0.9950266029i 0.5427101048
```

and,
```Table 16

r = 1, [lambda] = 1.01 n = 10

Nodes Weights

0.8432691294 + [+ or -] 0.5374915584i 0.0884264387
-0.9380253965 + [+ or -] 0.3465665239i 0.8084274344
0.6342873464 + [+ or -] 0.7730973821i 0.8785194976
0.4228259917 + [+ or -] 0.906210892i 0.3244615048
-0.1278951430 + [+ or -] 0.9917876952i 0.6277760750
```

Finally the approximation of the two integrals are displayed in:
```Table 17

r = 1 Approximation Exact value

[lambda] = 1.01 2.03543204774 ... 2.03543204817 ...
[lambda] = 1.1 1.67594127382 ... 1.67594127475 ...
```

Again, numerical results offered by both methods are quite similar.

Acknowledgements. The authors thank the referee for valuable suggestions which have contributed to improve the final form of this paper. The work of this author was partially supported by the Scientific Research Proyect BFM2001-3411 of the Spanish Ministry of Science and Technology.

REFERENCES

[1] N. I. AKHIEZER, The classical moment problem and some related questions in analysis, Hafner, New York, 1965.

[2] A. BULTHEEL, C. DIAZ-MENDOZA, P. GONZALEZ-VERA AND R. ORIVE, On the convergence of certain Gauss-type quadrature formulas for unbounded intervals, Math. Comp., 69 (2000), pp. 721-747.

[3] A. BULTHEEL, P. GONZALEZ-VERA, E. HENDRIKSEN AND O. NASTAD, Orthogonal rational functions and quadrature on the real half-line, Journal of Complexity, to appear.

[4] A. BULTHEEL, C. DIAZ-MENDOZA, P. GONZALEZ-VERA AND R. ORIVE, Quadrature on the half-line and two point Pade approximants to Stieljes functions, Part II: Convergence, J. Comput. Appl. Math., 77 (1997), pp. 53-76.

[5] A. BULTHEEL, C. DIAZ-MENDOZA, P. GONZALEZ-VERA AND R. ORIVE, Orthogonal Laurent polynomials and quadrature formulas for unbounded intervals: L Gauss-type formulas, Rocky Mountain J. Mathematics, to appear.

[6] A. BULTHEEL, P. GONZALEZ-VERA AND R. ORIVE, Quadrature on the half-line and two point Pade approximants to Stieljes functions, Part I. Algebraic aspects, J. Comput. Appl. Math., 65 (1996), pp. 57-72.

[7] L. COCHRAN AND S. CLEMENT COOPER, Orthogonal Laurent polynomials on the real line, in S. Clement Cooper, and W.J. Thron, eds., Lecture Notes in Pure and Applied Mathematics., vol. 154, Marcel Dekker, 1994, pp. 47-100.

[8] L. DARUIS, P. GONZALEZ-VERA AND O. NJASTAD., Szego quadrature formulas for certain Jacobi-type weight functions, Mathematics of Computation, 71 (2001), No. 238, pp. 683-701.

[9] L. DARUIS AND P. GONZALEZ-VERA, Interpolatory quadrature formulas on the unit circle for Chebyshev weight functions, Numer. Math., 90 (2002), pp. 641-664.

[10] P. DELSARTE, Y. GENIN, The Split Levinson Algorithm, IEEE Trans Acoust. Spreech Signal Proc., 34 (1986), pp. 470-478.

[11] W. GAUTSCHI, A survey of Gauss-Christoffel quadrature formulae, in Butzer, P.L., Feher, F., eds., E.B. Christoffel--The influence of his work in mathematics and the physical sciences, Birkhauser Verlag, Basel, 1981, pp. 72-147.

[12] W. GAUTSCHI, Gauss-type quadrature rules for rational functions, in H. Brass, G. Hammerlin, eds., Numerical Integration IV, International Series of Numerical Mathematics, vol. 112, Birkhauser, Basel, 1993, pp. 111-130.

[13] P. GONZALEZ-VERA, O. NJASTAD AND J.C. SANTOS-LEON, Some results about numerical quadrature on the unit circle, Advances in Computational Mathematics, 5 (1996), pp. 297-328.

[14] W.B. GRAGG, Positive definite Toeplitz matrices, the Arnoldi process for isometric operators, and Gaussian quadrature on the Unit Circle, J. Comput. Appl. Math., 46 (1993), pp. 183-198.

[15] P.E. GUSTAFSON AND B.A. HAGLER, Gaussian quadrature rules and numerical examples for strong extensions of mass distribution functions, J. Comput. Appl. Math., 105 (1999), pp. 317-326.

[16] W. B. JONES, O. NJASTAD AND W.J. THRON, Two-point Pade expansions for a family of analytic functions, J. Comput. Appl. Math., 9 (1983), pp. 105-123.

[17] W.B. JONES, O. NJASTAD AND W.J. THRON, Moment theory, orthogonal polynomials, quadrature, and continued fractions associated with the unit circle, Bull. London Math. Soc., 21 (1989), pp. 113-152.

[18] W.B. JONES, W.J. THRON AND H. WAADELAND, A strong Stieljes moment problem, Trans. of the AMS, 261 (1980), pp. 503-528.

[19] W.B. JONES AND W.J. THRON, Orthogonal Laurent polynomials and Gaussian quadrature, in Quantum Mechanics in Mathematics, Chemistry and Physics, K.E. Gustafson et al., eds., Plenum Publ. Comp., 1981, pp. 449-455.

[20] V.I. KRYLOV, Approximate Calculation of Integrals, The MacMillan Company, New York, 1962.

[21] N. LEVINSON, The Wiener RMS (root mean square) error criterion in filter design and prediction, J.Math. Phys., 25 (1947), pp. 261-278.

[22] G. LOPEZ-LAGOMASINO AND A. MARTINEZ-FINKELSHTEIN, Rate of convergence of two-point Pade approximants and logarithmic asymptotics of Laurent-type orthogonal polynomials, Constr. Approx., 11 (1995), pp. 255-286.

[23] G. LOPEZ-LAGOMASINO AND J. ILLAN-GONZALEZ, Quadrature formulas for unbounded intervals, Rev. Cienc. Mat., 3 (1982), pp. 29-47.

[24] O. NJASTAD AND W.J. THRON, The theory of sequences of orthogonal L-polynomials, in Pade approximants and Continued Fractions, H. Waadeland et a]., eds., Det Kongelige Norske Videnskabers Selskab Skrifter, No. 1, 1983, pp. 54-91.

[25] J.C. SANTOS, Computation of integrals over the unit circle with nearby poles, Appl. Numer. Math., 36 (2001), No. 2-3, pp. 179-195.

[26] A. SRI RANDA, Another quadrature rule of highest algebraic degree of precision, Numer. Math., 68 (1991), pp. 283-294.

[27] A. SRI RANDA, J-Fractions and strong moment problems, in W.J. Thron, ed., Analytic Theory of Continued Fractions, II, Procedings Pidochry and Aviemore 1985. Lecture Notes in Mathematics Springer-Verlag, No. 1199, 1986, pp. 269-284.

[28] G. SZEGO, Orthogonal polynomials, A.M.S., 4th ed., vol. 23, Colloquium Publications, 1975.

[29] G. SZEGO, On biorthogonal systems of trigonometric polynomials, Magyar. Tud. Akad. Kutato Int. Kozl., 8 (1963), pp. 255-273.

[30] W.J. THRON, L-polynomials orthogonal on the unit circle, in Nonlinear Methods and Rational Approximation, A. Cuyt, ed., D. Reidel Publishing Company, 1988, pp. 271-278.

* Received July 30, 2002. Accepted for publication May 10, 2003. Communicated by J. Arvesu.

RUYMAN CRUZ-BARROSO AND PABLO GONZALEZ-VERA ([dagger])

([dagger]) Department of Mathematical Analysis, La Laguna University. 38271 La Laguna. Tenerife. Spain.
COPYRIGHT 2005 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.