G Beylkin

ON THE REPRESENTATION OF OPERATORS IN BASES OF COMPACTLY SUPPORTED WAVELETS

Keywords: wavelets, differential operators, Hilbert transform, fractional derivatives, pseudodifferential operators, shift operators, numerical algorithms

This paper describes exact and explicit representations of the differential operators, d n /dx n , n = 1, 2, • • •, in orthonormal bases of compactly supported wavelets as well as the representations of the Hilbert transform and fractional derivatives. The method of computing these representations is directly applicable to multidimensional convolution operators.

Also, sparse representations of shift operators in orthonormal bases of compactly supported wavelets are discussed and a fast algorithm requiring O(N log N ) operations for computing the wavelet coefficients of all N circulant shifts of a vector of the length N = 2 n is constructed. As an example of an application of this algorithm, it is shown that the storage requirements of the fast algorithm for applying the standard form of a pseudodifferential operator to a vector (see [G. Beylkin, R. R. Coifman, and V. Rokhlin, Comm. Pure. Appl. Math., 44 (1991), pp. 141-183]) may be reduced from O(N ) to O(log 2 N ) significant entries.

Introduction.

In [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF] Daubechies introduced compactly supported wavelets which proved to be very useful in numerical analysis [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF]. In this paper we find exact and explicit representations of several basic operators (derivatives, Hilbert transform, shifts, etc.) in orthonormal bases of compactly supported wavelets. We also present an O(N log N ) algorithm for computing the wavelet coefficients of all N circulant shifts of a vector of the length N = 2 n .

Throughout this paper we only compute the nonstandard forms of operators since it is a simple matter to obtain a standard form from the nonstandard form [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF]. Meyer [START_REF] Meyer | Le calcul scientifique, les ondelettes et filtres miroirs en quadrature[END_REF], following [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF], considered several examples of nonstandard forms of basic operators from a general point of view. It is possible, however, to compute the nonstandard forms of many important operators explicitly.

First, we explicitly compute the nonstandard form of the operator d/dx. The set of coefficients that defines all nonzero entries of the nonstandard form appears as the solution to a system of linear algebraic equations. This system, in turn, arises as a consequence of the recursive definition of the wavelet bases. The operator d n /dx n is treated similarly to d/dx.

The computation of the nonstandard forms of many other operators reduces to solving a simple system of linear algebraic equations. Among such operators are fractional derivatives, Hilbert and Riesz transforms, and other operators for which analytic expressions are available. For convolution operators, there are significant simplifications in computing the nonstandard form since the vanishing moments of the autocorrelation function of the scaling function simplify the quadrature formulas. Moreover, by solving a system of linear algebraic equations combined with the asymptotics of wavelet coefficients, we arrive at an effective method for computing the nonstandard form of convolution operators. As examples, we compute the nonstandard forms of the Hilbert transform and fractional derivatives. The generalization of this method for multidimensional convolution operators is straightforward. Second, we compute the nonstandard form of the shift operator. This operator is important in practical applications of wavelets because the wavelet coefficients are not shift invariant. Since the nonstandard and standard forms of this operator are sparse and easy to compute, knowing these representations "compensates" for the lack of shift invariance. The wavelet expansion of shifts of vectors or of matrices may be obtained by applying the shift operator directly to the coefficients of the original expansion. The coefficients for the shift operators may be stored in advance and used as needed.

It is clear, however, that the particular manner in which sparseness of the shift operator may be exploited depends on the application and may be less straightforward than is indicated above. We present an example of such an application in numerical analysis. Observing that there are only N log 2 N distinct wavelet coefficients in the decomposition of all N circulant shifts of a vector of the length N = 2 n , we construct an O(N log N ) algorithm for computing all of these coefficients. Using this algorithm, we show that the storage requirements of the fast algorithm for applying the standard form of a pseudodifferential operator to a vector [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF] may be reduced from O(N log N ) to O(log 2 N ) significant entries.

2. Compactly supported wavelets. In this section, we briefly review the orthonormal bases of compactly supported wavelets and set our notation. For the details we refer to [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF].

The orthonormal basis of compactly supported wavelets of L 2 (R) is formed by the dilation and translation of a single function ψ(x),

ψ j,k (x) = 2 -j/2 ψ(2 -j x -k), (2.1)
where j, k ∈ Z. The function ψ(x) has a companion, the scaling function ϕ(x), and these functions satisfy the following relations:

ϕ(x) = √ 2 L-1 k=0 h k ϕ(2x -k), (2.2) ψ(x) = √ 2 L-1 k=0 g k ϕ(2x -k), (2.3)
where

g k = (-1) k h L-k-1 , k = 0, • • • , L -1, (2.4) and +∞ -∞ ϕ(x)dx = 1. (2.5)
In addition, the function ψ has M vanishing moments

+∞ -∞ ψ(x)x m dx = 0, m = 0, • • • , M -1. (2.6)
The number L of coefficients in (2.2) and (2.3) is related to the number of vanishing moments M , and for the wavelets in [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF], L = 2M . If additional conditions are imposed (see [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF] for an example), then the relation might be different, but L is always even.

The wavelet basis induces a multiresolution analysis on L 2 (R) [START_REF]Ondelettes et functions splines[END_REF], [START_REF] Mallat | Multiresolution approximation and wavelets[END_REF], i.e., the decomposition of the Hilbert space L 2 (R) into a chain of closed subspaces

• • • ⊂ V 2 ⊂ V 1 ⊂ V 0 ⊂ V -1 ⊂ V -2 ⊂ • • • (2.7) such that j∈Z V j = {0}, j∈Z V j = L 2 (R). (2.8)
By defining W j as an orthogonal complement of V j in V j-1 ,

V j-1 = V j ⊕ W j , (2.9) the space L 2 (R) is represented as a direct sum L 2 (R) = j∈Z W j . (2.10)
On each fixed scale j, the wavelets {ψ j,k (x)} k∈Z form an orthonormal basis of W j and the functions {ϕ j,k (x) = 2 -j/2 ϕ(2 -j x -k)} k∈Z form an orthonormal basis of V j .

The coefficients

H = {h k } k=L-1 k=0 and G = {g k } k=L-1 k=0 in (2.
2) and (2.3) are quadrature mirror filters. Once the filter H has been chosen, it completely determines the functions ψ and ϕ. Let us define the 2π-periodic function

m 0 (ξ) = 2 -1/2 k=L-1 k=0 h k e ikξ (2.11)
where {h k } k=L-1 k=0 are the coefficients of the filter H. The function m 0 (ξ) satisfies the equation

|m 0 (ξ)| 2 + |m 0 (ξ + π)| 2 = 1. (2.12)
The following lemma characterizes trigonometric polynomial solutions of (2.12) which correspond to the orthonormal bases of compactly supported wavelets with vanishing moments.

Lemma 1 (Daubechies [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF]). Any trigonometric polynomial solution m 0 (ξ) of (2.12) is of the form (2.13) where M ≥ 1 is the number of vanishing moments, and where Q is a polynomial such that

m 0 (ξ) = 1 2 (1 + e iξ ) M Q(e iξ ),
|Q(e iξ )| 2 = P (sin 2 1 2 ξ) + sin 2M ( 1 2 ξ) R( 1 2 cos ξ), (2.14)
where

P (y) = k=M -1 k=0 M -1 + k k y k , (2.15)
and R is an odd polynomial such that 1) . (2.17)

0 ≤ P (y) + y M R( 1 2 -y) for 0 ≤ y ≤ 1, (2.16) and sup 0≤y≤1 P (y) + y M R( 1 2 -y) < 2 2(M -
3. The operator d/dx in wavelet bases. In this section we construct the nonstandard form of the operator d/dx. The nonstandard form [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF] is a representation of an operator T as a chain of triplets

T = {A j , B j , Γ j } j∈Z (3.1)
acting on the subspaces V j and W j ,

A j : W j → W j , (3.2)
B j : V j → W j , (3.3)
Γ j : W j → V j . (3.4)
The operators {A j , B j , Γ j } j∈Z are defined as A j = Q j T Q j , B j = Q j T P j , and Γ j = P j T Q j , where P j is the projection operator on the subspace V j and Q j = P j-1 -P j is the projection operator on the subspace W j .

The matrix elements α j il , β j il , γ j il of A j , B j , Γ j , and r j il of T j = P j T P j , i, l, j ∈ Z, for the operator d/dx are easily computed as

α j il = 2 -j ∞ -∞ ψ(2 -j x -i) ψ (2 -j x -l) 2 -j dx = 2 -j α i-l , (3.5)
β j il = 2 -j ∞ -∞ ψ(2 -j x -i) ϕ (2 -j x -l) 2 -j dx = 2 -j β i-l , (3.6) γ j il = 2 -j ∞ -∞ ϕ(2 -j x -i) ψ (2 -j x -l) 2 -j dx = 2 -j γ i-l , (3.7) and r j il = 2 -j ∞ -∞ ϕ(2 -j x -i) ϕ (2 -j x -l) 2 -j dx = 2 -j r i-l , (3.8)
where

α l = +∞ -∞ ψ(x -l) d dx ψ(x) dx, (3.9) β l = +∞ -∞ ψ(x -l) d dx ϕ(x) dx, (3.10) γ l = +∞ -∞ ϕ(x -l) d dx ψ(x) dx. (3.11) and r l = +∞ -∞ ϕ(x -l) d dx ϕ(x) dx. (3.12)
Moreover, using (2.2) and (2.3) we have

α i = 2 L-1 k=0 L-1 k =0 g k g k r 2i+k-k , (3.13) β i = 2 L-1 k=0 L-1 k =0 g k h k r 2i+k-k , (3.14) and γ i = 2 L-1 k=0 L-1 k =0 h k g k r 2i+k-k , (3.15)
and, therefore, the representation of d/dx is completely determined by r l in (3.12) or, in other words, by the representation of d/dx on the subspace V 0 .

Rewriting (3.12) in terms of φ(ξ), where

φ(ξ) = 1 √ 2π +∞ -∞ ϕ(x) e ixξ dx, (3.16) we obtain r l = +∞ -∞ (-iξ) | φ(ξ)| 2 e -ilξ dξ. (3.17)
In order to compute the coefficients r l we first note that any trigonometric polynomial m 0 (ξ) satisfying (2.12) is such that

|m 0 (ξ)| 2 = 1 2 + 1 2 L/2 k=1 a 2k-1 cos(2k -1)ξ, (3.18)
where a n are the autocorrelation coefficients of

H = {h k } k=L-1 k=0 , a n = 2 L-1-n i=0 h i h i+n , n = 1, • • • , L -1. (3.19)
The autocorrelation coefficients a n with even indices are zero,

a 2k = 0, k = 1, • • • , L/2 -1. (3.20)
To prove this assertion we compute |m 0 (ξ)| 2 using (2.11) and obtain

|m 0 (ξ)| 2 = 1 2 + 1 2 L-1 n=1
a n cos nξ, (3.21) where a n are given in (3.19). Computing |m 0 (ξ + π)| 2 , we have We prove the following: Proposition 1.

|m 0 (ξ + π)| 2 = 1 2 - 1 2 L/2 k=1 a 2k-1 cos(2k -1)ξ + 1 2 L/2-
(1) If the integrals in (3.12) or (3.17) exist, then the coefficients r l in (3.12) satisfy the following system of linear algebraic equations: where the coefficients a 2k-1 are given in (3.19).

r l = 2   r 2l + 1 2 L/2 k=1 a 2k-1 (r 2l-2k+1 + r 2l+2k-1 )   , (3.
(2) If M ≥ 2, then equations (3.24) and (3.25) have a unique solution with a finite number of nonzero r l , namely, r l = 0 for -L + 2 ≤ l ≤ L -2 and

r l = -r -l . (3.26) Remark 1. If M = 1,
then equations (3.24) and (3.25) have a unique solution but the integrals in (3.12) or (3.17) may not be absolutely convergent. Let us consider Example 3.2 of [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF], where L = 4 and

h 0 = 2 -1/2 ν(ν -1) ν 2 + 1 , h 1 = 2 -1/2 1 -ν ν 2 + 1 , h 2 = 2 -1/2 ν + 1 ν 2 + 1 , h 3 = 2 -1/2 ν(ν + 1) ν 2 + 1 ,
where ν is an arbitrary real number. We have

a 1 = 1 + 3ν 2 (ν 2 + 1) 2 , a 3 = ν 2 (ν 2 -1) (ν 2 + 1) 2 ,
and

r 1 = - (1 + ν 2 ) 2 2(3ν 4 + 1) , r 2 = ν 2 (1 -ν 2 ) 2(3ν 4 + 1)
.

The parameter ν can be chosen so that the Fourier transform φ(ξ) does not have the sufficient decay to insure the absolute convergence of the integral (3.17).

For the Haar basis (h 1 = h 2 = 2 -1/2 ), a 1 = 1 and r 1 = -1 2 , and thus, we obtain the simplest finite difference operator ( 12 , 0, -1 2 ). In this case the function ϕ is not continuous and

φ(ξ) = 1 √ 2π sin 1 2 ξ 1 2 ξ e i 1 2 ξ ,
so that the integral in (3.17) is not absolutely convergent.

Proof of Proposition 1. Using (2.2) for both ϕ(x -l) and d dx ϕ(x) in (3.12) we obtain

r i = 2 L-1 k=0 L-1 l=0 h k h l +∞ -∞ ϕ(2x -2i -k) ϕ (2x -l) 2 dx (3.27)
and hence,

r i = 2 L-1 k=0 L-1 l=0 h k h l r 2i+k-l . (3.28)
Substituting l = k -m, we rewrite (3.28) as

r i = 2 L-1 k=0 k-L+1 m=k h k h k-m r 2i+m . (3.29)
Changing the order of summation in (3.29) and using the fact that

L-1 k=0 h 2 k = 1, we arrive at r l = 2r 2l + L-1 n=1 a n (r 2l-n + r 2l+n ), l ∈ Z, (3.30)
where a n are given in (3.19). Using (3.20), we obtain (3.24) from (3.30).

In order to obtain (3.25) we use the following relation:

l=+∞ l=-∞ l m ϕ(x -l) = x m + l=m l=1 (-1) l m l M ϕ l x m-l , (3.31)
where

M ϕ l = +∞ -∞ ϕ(x) x l dx, where l = 1, • • • , m, (3.32)
are the moments of the function ϕ(x). We note that (3.31) is well known if all moments (3.32) are zero. The general statement follows simply on taking Fourier transforms and using Leibniz's rule. Using (3.12) and (3.31) with m = 1, we obtain (3.25).

If M ≥ 2, then | φ(ξ)| 2 |ξ| ≤ C(1 + |ξ|) -1-, (3.33)
where > 0, and hence, the integral in (3.17) is absolutely convergent. This assertion follows from Lemma 3.2 of [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF], where it is shown that

| φ(ξ)| ≤ C(1 + |ξ|) -M +log 2 B , (3.34) where B = sup ξ∈R |Q(e iξ )|.
Due to condition (2.17), we have log 2 B = M -1 -with some > 0. The existence of the solution of the system of equations (3.24) and (3.25) follows from the existence of the integral in (3.17). Since the scaling function ϕ has a compact support there are only a finite number of nonzero coefficients r l . The specific interval -L+2 ≤ l ≤ L-2 is obtained by the direct examination of (3.24).

Let us show now that l r l = 0. (3.35) Multiplying (3.24) by e ilξ and summing over l, we obtain

r(ξ) = 2   reven (ξ/2) + 1 2 rodd (ξ/2) L/2 k=1 a 2k-1 (e -i(2k-1)ξ/2 + e i(2k-1)ξ/2 )   , (3.
r(ξ) = r(ξ/2) + r(ξ/2 + π) + (2|m 0 (ξ/2)| 2 -1)(r(ξ/2) -r(ξ/2 + π)) . (3.42)
Finally, using (2.12) we arrive at

r(ξ) = 2( |m 0 (ξ/2)| 2 r(ξ/2) + |m 0 (ξ/2) + π)| 2 r(ξ/2) + π)). (3.43)
Setting ξ = 0 in (3.43), we obtain r(0) = 2r(0) and thus, (3.35).

Uniqueness of the solution of (3.24) and (3.25) follows from the uniqueness of the representation of d/dx. Given the solution r l of (3.24) and (3.25) we consider the operator T j defined by these coefficients on the subspace V j and apply it to a sufficiently smooth function f . Since r j l = 2 -j r l (3.8), we have

(T j f )(x) = k∈Z 2 -j l r l f j,k-l ϕ j,k (x), (3.44) where f j,k-l = 2 -j/2 +∞ -∞ f (x) ϕ(2 -j x -k + l) dx. (3.45) Rewriting (3.45) f j,k-l = 2 -j/2 +∞ -∞ f (x -2 j l) ϕ(2 -j x -k) dx, (3.46)
and expanding f (x -2 j l) in the Taylor series at the point x, we have

f j,k-l = +∞ -∞ f (x) ϕ j,k (x) dx -2 j l +∞ -∞ f (x) ϕ j,k (x) dx +2 2j l 2 2 +∞ -∞ f (x) ϕ j,k (x) dx, (3.47)
where x = x(x, x -2 j l) and |x -x| ≤ 2 j l. Substituting (3.47) in (3.44) and using (3.35) and (3.25), we obtain

(T j f )(x) = k∈Z +∞ -∞ f (x) ϕ j,k (x) dx ϕ j,k (x) +2 j k∈Z 1 2 l r l l 2 +∞ -∞ f (x) ϕ j,k (x) dx ϕ j,k (x). (3.48)
It is clear that as j → -∞, operators T j and d/dx coincide on smooth functions. Using standard arguments it is easy to prove that T -∞ = d/dx and hence, the solution to (3.24) and (3.25) is unique. The relation (3.26) follows now from (3.17).

Remark 2. We note that expressions (3.13) and (3.14) for α l and β l (γ l = -β -l ) may be simplified by changing the order of summation in (3.13) and (3.14) and by introducing the correlation coefficients 2 L-1-n i=0 and2 L-1-n i=0 g i g i+n . The expression for α l is especially simple: α l = 4r 2l -r l .

g i h i+n , 2 L-1-n i=0 h i g i+n
Examples. For examples we will use Daubechies' wavelets constructed in [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF]. First, let us compute the coefficients a 2k-1 , k = 1, • • • , M , where M is the number of vanishing moments and L = 2M . Using relation (4.22) of [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF],

|m 0 (ξ)| 2 = 1 - (2M -1)! [(M -1)!] 2 2 2M -1 ξ 0 sin 2M -1 ξ dξ, (3.49)
we find, by computing

ξ 0 sin 2M -1 ξ dξ, that |m 0 (ξ)| 2 = 1 2 + 1 2 C M M m=1 (-1) m-1 cos(2m -1)ξ (M -m)! (M + m -1)! (2m -1) , (3.50)
where

C M = (2M -1)! (M -1)! 4 M -1 2 . (3.51)
Thus, by comparing (3.50) and (3.18), we have

a 2m-1 = (-1) m-1 C M (M -m)! (M + m -1)! (2m -1)
, where m = 1,

• • • , M. (3.52)
The coefficients r e are rational numbers since they are solutions of a linear system with rational coefficients (a 2m-1 in (3.52) are rational by construction). We note that the coefficients r l are the same for all Daubechies' wavelets with a fixed number of vanishing moments M , notwithstanding the fact that there are several different wavelet bases for a given M (depending on the choice of the roots of polynomials in the construction described in [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF]).

Solving the equations of Proposition 1, we present the results for Daubechies' wavelets with M = 2, 3, 4, 5, 6.

1. M = 2.

a 1 = 9 8 , a 3 = - 1 8 ,
and

r 1 = - 2 3 , r 2 = 1 12 .
The coefficients (-1/12, 2/3, 0, -2/3, 1/12) of this example coincide with one of the standard choices of coefficients for numerical differentiation.

2. M = 3. Iterative algorithm for computing the coefficients r l . We also use an iterative algorithm as a way of solving equations (3.24) and (3.25). We start with r -1 = 0.5 and r 1 = -0.5, and iterate using (3.24) to recompute r l . Using (3.43), it is easy to verify that (3.25) and (3.26) are satisfied due to the choice of initialization. Table 1 was computed using this algorithm for Daubechies' wavelets with M = 5, 6, 7, 8, 9. It only displays the coefficients {r l } L-2 l=1 since r -l = -r l and r 0 = 0. 4. The operator d n /dx n in wavelet bases. Again, as in the case with the operator d/dx, the nonstandard form of the operator d n /dx n is completely determined by its representation on the subspace V 0 , i.e., by the coefficients

a 1 =
r (n) l = +∞ -∞ ϕ(x -l) d n dx n ϕ(x) dx, l ∈ Z, (4.1)
or, alternatively,

r (n) l = +∞ -∞ (-iξ) n | φ(ξ)| 2 e -ilξ dξ (4.2)
if the integrals in (4.1) or (4.2) exist (see also Remark 3). l , l ∈ Z satisfy the following system of linear algebraic equations:

r (n) l = 2 n   r 2l + 1 2 L/2 k=1 a 2k-1 (r (n) 2l-2k+1 + r (n) 2l+2k-1 )   , (4.3) and l l n r (n) l = (-1) n n!, (4.4)
where a 2k-1 are given in (3.19).

(2) Let M ≥ (n + 1)/2, where M is the number of vanishing moments in (2.6).

If the integrals in (4.1) or (4.2) exist, then the equations (4.3) and (4.4) have a unique solution with a finite number of nonzero coefficients r and for odd n r

(n) l , namely, r (n) l = 0 for -L + 2 ≤ l ≤ L -2, such that for even n r (n) l = r (n) -l , (4.5) l l 2ñ r (n) l = 0, ñ = 1, • • • , n/2 -1, (4.6
(n) l = -r (n) -l , (4.8) l l 2ñ-1 r (n) l = 0, ñ = 1, • • • , (n -1)/2. (4.9)
The proof of Proposition 2 is completely analogous to that of Proposition 1. Remark 3. The linear system in Proposition 2 may have a unique solution whereas integrals (4.1) and (4.2) are not absolutely convergent. A case in point is the Daubechies' wavelet with M = 2. The representation of the first derivative in this basis is described in the previous section. Equations (4.3) and (4.4) do not have a solution for the second derivative n = 2. However, the system of equations (4.3) and (4.4) has a solution for the third derivative n = 3. We have

a 1 = 9 8 , a 3 = - 1 8 , and
r -2 = - 1 2 , r -1 = 1, r 0 = 0, r 1 = -1, r 2 = 1 2 .
The set of coefficients (-1/2, 1, 0, -1, 1/2) is one of the standard choices of finite difference coefficients for the third derivative. We note that among the wavelets with L = 4, the wavelets with two vanishing moments M = 2 do not have the best Hölder exponent (see [START_REF] Daubechies | Two-scale difference equations, II. Local regularity, infinite products of matrices and fractals[END_REF]), but the representation of the third derivative exists only if the number of vanishing moments M = 2.

Remark 4. Let us derive an equation generalizing to (3.43) for d n /dx n directly from (4.2). We rewrite (4.2) as

r (n) l = 2π 0 k∈Z | φ(ξ + 2πk)| 2 (-i) n (ξ + 2πk
) n e -ilξ dξ.

(M 0 f )(ξ) = |m 0 (ξ/2)| 2 f (ξ/2) + |m 0 (ξ/2) + π)| 2 f (ξ/2) + π), (4.15)
we rewrite (4.14) as

M 0 r = 2 -n r. (4.16)
Thus, r is an eigenvector of the operator M 0 corresponding to the eigenvalue 2 -n and, therefore, finding the representation of the derivatives in the wavelet basis is equivalent to finding trigonometric polynomial solutions of (4.16) and vice versa. (The operator M 0 is also introduced in [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF] and [START_REF] Lawton | Necessary and sufficient conditions for constructing orthonormal wavelet bases[END_REF], where the problem (4.16) with eigenvalue 1 is considered.) Remark 5. While theoretically it is well understood that the derivative operators (or, more generally, operators with homogeneous symbols) have an explicit diagonal preconditioner in wavelet bases, the numerical evidence illustrating this fact is of interest, since it represents one of the advantages of computing in the wavelet bases.

If an operator has a null space (the actual null space or a null space for a given accuracy), then by the condition number we understand the ratio of the largest singular value to the smallest singular value above the threshold of accuracy. Thus, we include the situation where the operator may be preconditioned only on a subspace. We note that the preconditioning described here addresses the problem of ill conditioning due only to the unfavorable homogeneity of the symbol and does not affect ill conditioning due to other causes.

For periodized derivative operators the bound on the condition number depends only on the particular choice of the wavelet basis. After applying such a preconditioner, the condition number κ p of the operator is uniformly bounded with respect to the size of the matrix. We recall that the condition number controls the rate of convergence of a number of iterative algorithms; for example, the number of iterations of the conjugate gradient method is O( √ κ p ). Thus, this remark implies a completely new outlook on a number of numerical methods, a topic we will address elsewhere.

We present here two tables illustrating such preconditioning applied to the standard form of the second derivative (see [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF] on how to compute the standard form from the nonstandard form). In the following examples the standard form of the periodized second derivative D 2 of size N × N , where N = 2 n , is preconditioned by the diagonal matrix P , D p 2 = P D 2 P where P il = δ il 2 j , 1 ≤ j ≤ n, and where j is chosen depending on i, l so that N -N/2 j-1 + 1 ≤ i, l ≤ N -N/2 j , and

P N N = 2 n .
Tables 2 and3 compare the original condition number κ of D 2 and κ p of D p 2 . 5. Convolution operators. For convolution operators, the computation of the nonstandard form is considerably simpler than in the general case [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF]. We will demonstrate that the quadrature formulas for representing kernels of convolution operators on V 0 (see, e.g., Appendix B of [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF]) are of the simplest form due to the fact that the moments of the autocorrelation function of the scaling function ϕ vanish. Moreover, by combining the asymptotics of wavelet coefficients with the system of linear algebraic equations (similar to those in previous sections), we arrive at an effective method for computing representations of convolution operators. This method is especially simple if the symbol of the operator is homogeneous of some degree.

Let us assume that the matrix t (j-1) i-l (i, l ∈ Z) represents the operator P j-1 T P j-1 on the subspace V j-1 . To compute the matrix representation of P j T P j , we have the following formula (3.26) of [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF]: which easily reduces to

t (j) l = L-1 k=0 L-1 m=0 h k h m t (j-1) 2l+k-m , (5.1)
t (j) l = t (j-1) 2l + 1 2 L/2 k=0 a 2k-1 (t (j-1) 2l-2k+1 + t (j-1) 2l+2k-1 ), (5.2)
where the coefficients a 2k-1 are given in (3.19).

We also have

t (j) l = +∞ -∞ +∞ -∞
K(x -y) ϕ j,0 (y) ϕ j,l (x) dxdy, (5.3) and by changing the order of integration, we obtain

t (j) l = 2 j +∞ -∞ K(2 j (l -y)) Φ(y) dy, (5.4)
where Φ is the autocorrelation function of the scaling function ϕ,

Φ(y) = +∞ -∞ ϕ(x) ϕ(x -y) dx. (5.5)

Let us verify that

M m Φ = 1 i ∂ ξ m | φ(ξ)| 2 ξ=0 . (5.8)
Using (5.8) and the identity φ(ξ) = φ(ξ/2)m 0 (ξ/2) (see [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF]), it is clear that (5.7) holds provided that (5.9) or (due to (2.12))

1 i ∂ ξ m |m 0 (ξ)| 2 ξ=0 = 0 for 1 ≤ m ≤ 2M -1,
1 i ∂ ξ m |m 0 (ξ + π)| 2 ξ=0 = 0 for 0 ≤ m ≤ 2M -1.
(5.10) But formula (5.10) follows from the explicit representation in (2.13).

Remark 6. Equations (5.9) and (3.21) also imply that even moments of the coefficients a 2k-1 from (3.19) vanish, namely,

k=L/2 k=1 a 2k-1 (2k -1) 2m = 0 for 1 ≤ m ≤ M -1. (5.11)
Since the moments of the function Φ vanish equation (5.4) leads to a one-point quadrature formula for computing the representation of convolution operators on the finest scale. This formula is obtained in exactly the same manner as for the special choice of the wavelet basis described in [2, eqns. (3.8)-(3.12)], where the shifted moments of the function ϕ vanish; we refer to this paper for the details.

Here we introduce a different approach for computing representations of convolution operators in the wavelet basis which consists of solving the system of linear algebraic equations (5.2) subject to asymptotic conditions. This method is especially simple if the symbol of the operator is homogeneous of some degree since in this case the operator is completely defined by its representation on V 0 . We consider two examples of such operators, the Hilbert transform and the operator of fractional differentiation (or antidifferentiation).

The Hilbert transform. We apply our method to the computation of the nonstandard form of the Hilbert transform

g(x) = (Hf )(y) = 1 π p.v. ∞ -∞
f (s) s -x ds, (5.12) where p.v. denotes a principal value at s = x.

The representation of H on V 0 is defined by the coefficients

r l = ∞ -∞ ϕ(x -l) (Hϕ)(x) dx, l ∈ Z, (5.13)
which, in turn, completely define all other coefficients of the nonstandard form. Namely, H = {A j , B j , Γ j } j∈Z , A j = A 0 , B j = B 0 , and Γ j = Γ 0 , where matrix elements α i-l , β i-l , and γ i-l of A 0 , B 0 , and Γ 0 are computed from the coefficients r l ,

α i = L-1 k=0 L-1 k =0 g k g k r 2i+k-k , (5.14)
β i = L-1 k=0 L-1 k =0 g k h k r 2i+k-k , (5.15)
and

γ i = L-1 k=0 L-1 k =0 h k g k r 2i+k-k . (5.16)
The coefficients r l , l ∈ Z in (5.13) satisfy the following system of linear algebraic equations:

r l = r 2l + 1 2
L/2 k=1 a 2k-1 (r 2l-2k+1 + r 2l+2k-1 ), (5.17) where the coefficients a 2k-1 are given in (3.19). Using (5.4), (5.6), and (5.7), we obtain the asymptotics of r l for large l,

r l = - 1 πl + O 1 l 2M .
(5.18) By rewriting (5.13) in terms of φ(ξ),

r l = -2 ∞ 0 | φ(ξ)| 2 sin(lξ) dξ. (5.19)
we obtain r l = -r -l and set r 0 = 0. We note that the coefficient r 0 cannot be determined from equations (5.17) and (5.18). Solving (5.17) with the asymptotic condition (5.18), we compute the coefficients r l , l = 0 with any prescribed accuracy. We note that the generalization for computing the coefficients of Riesz transforms in higher dimensions is straightforward.

Example. We compute (see Table 4) the coefficients r l of the Hilbert transform for Daubechies' wavelets with six vanishing moments with accuracy 10 -7 . The coefficients for l > 16 are obtained using asymptotics (5.18). (We note that r -l = -r l and r 0 = 0.) Fractional derivatives. We use the following definition of fractional derivatives:

(∂ α x f ) (x) = +∞ -∞ (x -y) -α-1 + Γ(-α) f (y)dy, (5.20)
where we consider α = 1, 2 • • •. If α < 0, then (5.20) defines fractional antiderivatives.

The representation of ∂ α x on V 0 is determined by the coefficients

r l = +∞ -∞ ϕ(x -l) (∂ α x ϕ) (x) dx, l ∈ Z, (5.21)
provided that this integral exists.

The nonstandard form ∂ α x = {A j , B j , Γ j } j∈Z is computed via A j = 2 -αj A 0 , B j = 2 -αj B 0 , and Γ j = 2 -αj Γ 0 , where matrix elements α i-l , β i-l , and γ i-l of A 0 , B 0 , and Γ 0 are obtained from the coefficients r l ,

α i = 2 α L-1 k=0 L-1 k =0 g k g k r 2i+k-k , (5.22)
β i = 2 α L-1 k=0 L-1 k =0 g k h k r 2i+k-k , (5.23)
and

γ i = 2 α L-1 k=0 L-1 k =0 h k g k r 2i+k-k . (5.24)
It is easy to verify that the coefficients r l satisfy the following system of linear algebraic equations:

r l = 2 α   r 2l + 1 2 L/2 k=1 a 2k-1 (r 2l-2k+1 + r 2l+2k-1 )   , (5.25)
where the coefficients a 2k-1 are given in (3.19). Using (5.4), (5.6), and (5.7), we obtain the asymptotics of r l for large l,

r l = 1 Γ(-α) 1 l 1+α + O 1 l 1+α+2M
for l > 0, (5.26) r l = 0 for l < 0. (5.27)

Example. We compute (see Table 5) the coefficients r l of the fractional derivative with α = 0.5 for Daubechies' wavelets with six vanishing moments with accuracy 10 -7 . The coefficients for r l , l > 14, or l < -7 are obtained using asymptotics

r l = - 1 2 √ π 1 l 1+ 1 2 + O 1 l 13+ 1 2
for l > 0, (5.28) r l = 0 for l < 0.

(5.29) 6. Shift operator on V 0 and fast wavelet decomposition of all circulant shifts of a vector. Let us consider a shift by one on the subspace V 0 represented by the matrix

t (0) i-j = δ i-j,1 , (6.1)
where δ is the Kronecker symbol. Using (5.1) with the a n of (3.19) we have

t (0) l = δ l,1 , t (1) l
= 1 2 a |2l-1| , • • • . (6.2)
The only nonzero coefficients t (j) l on each scale j are those with indices -L + 2 ≤ l ≤ L -2. Also, t (j) l → δ l,0 as j → ∞. As an example, the following Table 6 contains the coefficients t (j) l , j = 1, 2, • • • , 8, for the shift operator in Daubechies' wavelet basis with three vanishing moments.

We note that the shift by an integer other than one is treated similarly. However, if the absolute value of the shift is greater than L -2, then, on the first several scales j, there are nonzero coefficients t (j) l with l outside the interval |l| ≤ L -2. As j increases, all the nonzero coefficients t (j) l will have indices in the interval |l| ≤ L -1. The importance of the shift operator stems from the fact that the coefficients of wavelet transforms are not shift invariant. However, as we have just demonstrated, the nonstandard (and, therefore, the standard) forms of the shift operator are sparse and easy to compute. By applying these sparse representations directly to the wavelet coefficients, in many applications we can effectively compensate the absence of the shift invariance of the wavelet transforms. For example, if the representation of a vector in the wavelet basis is sparse, there is a corresponding reduction in the number of operations required to shift such a vector. Specifically, in image processing the shift for the shift operators can be stored in advance and used as needed. It is clear, however, that the method of using sparseness of the shift operator depends on the specific application and may be less straightforward than is indicated above. The following is an example of an application where, instead of computing shift operators, we compute all possible shifts. We describe a fast algorithm for the wavelet decomposition of all circulant shifts of a vector and then show how it may be used to reduce storage requirements of one of the algorithms of [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF].

We recall that the decomposition of a vector of length N = 2 n into a wavelet basis requires O(N ) operations. Since the coefficients are not shift invariant, the computation of the wavelet expansion of all N circulant shifts of a vector appears to require O(N 2 ) operations.

We notice, however, that there are only N log 2 (N ) distinct coefficients and present here a simple algorithm to compute the wavelet expansion of N circulant shifts of a vector in N log(N ) operations.

We recall the pyramid scheme

{s 0 k } -→ {s 1 k } -→ {s 2 k } -→ {s 3 k } • • • {d 1 k } {d 2 k } {d 3 k } • • • (6.3)
where the coefficients s 0

k for k = 1, 2, • • • , N are given, s j k = n=L-1 n=0 h n s j-1 n+2k-1 , (6.4) d j k = n=L-1 n=0
g n s j-1 n+2k-1 , (6.5) and s j k and d j k are periodic sequences with the period 2 n-j , j = 0, 1, • • • , n. In the pyramid scheme (6.3), on each scale j we compute one vector of differences {d j k } k=2 n-j k=1 and one vector of averages {s j k } k=2 n-j k=1

. Instead, let us compute on each scale j, (1 ≤ j ≤ n), 2 j vectors of differences and 2 j vectors of averages. We proceed as follows: let s j-1 k , k = 1, • • • , 2 n-j be one of the vectors of averages on the previous scale j -1 and compute

s j k (0) = n=L-1 n=0
h n s j-1 n+2k-1 , (6.6)

s j k (1) = n=L-1 n=0
h n s j-1 n+2k , (6.7) and

d j k (0) = n=L-1 n=0
g n s j-1 n+2k-1 , (6.8)

d j k (1) = n=L-1 n=0
g n s j-1 n+2k . (6.9)

To compute the sum in (6.7) and (6.9), we shift by one the sequence s j-1 k in (6.6) and (6.8).

Thus, stepping from scale to scale we double the number of vectors of averages and of differences and, at the same time, halve the length of each of them. Therefore, the total number of operations in this computation is O(N log N ).

Let us organize the vectors of differences and averages as follows: on the first scale, j = 1, we set

v 1 = (d 1 k (0), d 1 k (1)) (6.10) and u 1 = (s 1 k (0), s 1 k (1)
), (6.11) where d 1 k (0), d 1 k (1), s 1 k (0), and s 1 k (1) are computed from s 0 k according to (6.6)-(6.9). On the second scale, j = 2, we set

v 2 = (d 2 k (00), d 2 k (01), d 2 k (10), d 2 k ( 11
) contains all the coefficients for odd shifts by 1, 3, • • • , 2 n -1 (see (6.9)). By the same token, the periodic sequences s 1 k (0) and s 1 k (1) contain all possible coefficients for even and odd circulant shifts of s 0 k (see (6.6) and (6.7)). Repeating this procedure on the next scale for both s 1 k (0) and s 1 k (1), we again obtain all possible coefficients for odd and even shifts which we collect in v 2 and u 2 , etc.

While the vectors v 1 , v 2 , • • • , v n contain all the coefficients, these coefficients are not organized sequentially. In order to access them, we generate two tables i loc (i s , j) and i b (i s , j) in O(N log N ) operations as follows. For each shift i s , 0 ≤ i s ≤ N -1 of the vector s 0 k , k = 1, 2, • • • , N , let us write the binary expansion of i s , i s = l=n-1 l=0 l 2 l , (6.14) where l = 0, 1. For a fixed scale j, 1 ≤ j ≤ n, we compute i loc (i s , j) = where i b (i s , j) = 0 if j = n. The number i b (i s , j) points to the begining of the subvector of differences in v j . Namely, the subvector of v j has indices between i b (i s , j) + 1 and i b (i s , j) + 2 n-j . Within this subvector (which is treated as a periodic vector with the period 2 n-j ) the number i loc (i s , j) points to the first element.

For all scales j, 1 ≤ j ≤ n, and shifts i s , 0 ≤ i s ≤ N -1, we compute two tables in (6.15) and (6.16). These tables give us the direct access to the coefficients in vectors v 1 , v 2 , • • • , v n for a constant cost per element.

We now briefly describe one of the applications of the algorithm for the fast wavelet decomposition of all circulant shifts of a vector in numerical analysis. The algorithms of [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF] are designed to evaluate the Calderon-Zygmund or a pseudodifferential operator T with kernel K(x, y), g(x) = +∞ -∞ K(x, y) f (y) dy (6.17) by constructing (for any fixed accuracy) its sparse nonstandard or standard form and thereby, reducing the cost of applying it to a function.

Let us rewrite (6.17 If the operator T is a convolution, then K(x, x -z) = K(z) is a function of z only. The nonstandard form of a convolution requires at most O(log N ) of storage (see the previous section), while the standard form of [START_REF] Beylkin | Fast wavelet transform and numerical algorithms I[END_REF] will contain O(N ) or O(N log N ) significant entries even for a convolution. Alternatively, the standard form of K(x, x -z) = K(z) in variables x and z for the convolution operators contains no more than O(log N ) significant entries for any fixed accuracy, since the kernel depends on one variable only.

If we now construct the standard form of K(x, x-z) in variables x and z for pseudodifferential operators (not necessarily convolutions), we obtain "super"-compression of the operator. Indeed, if these operators are represented in the form (6.18), then the dependence of the kernel K(x, •) on x is smooth and the number of significant entries in the standard form is of O(log 2 N ).

The apparent difficulty in computing via (6.18) is that it is necessary to compute the wavelet decomposition of f (x -z) for every x and thus, it appears to require O(N 2 ) operations. The algorithm of this section accomplishes this task in O(N log N ) operations. Therefore, the cost of evaluating (6.18) does not exceed O(N log N ). The advantage of such an algorithm is the reduced storage of O(log 2 N ) significant entries for the standard form, which is an important consideration for the multidimensional operators. We note that the extension of the algorithm for the fast wavelet decomposition of all circulant shifts of a vector to the multidimensional case is straightforward.

  Schlumberger-Doll Research, Old Quarry Road, Ridgefield, Connecticut 06877. Present address, Program in Applied Mathematics, University of Colorado at Boulder, Boulder, Colorado 80309-0526.

  ξ/2) = r(ξ/2) -r(ξ/2 + π), (3.41) and using (3.18), we obtain from(3.36)

  ξ + 2πk)| 2 (-i) n (ξ + 2πk) n , ) = m 0 (ξ/2) φ(ξ/2) (4.13) into the right-hand side of (4.11), and summing separately over even and odd indices in (4.11), we arrive at r(ξ) = 2 n ( |m 0 (ξ/2)| 2 r(ξ/2) + |m 0 (ξ/2) + π)| 2 r(ξ/2) + π)). (4.14) By considering the operator M 0 defined on 2π-periodic functions,

  m Φ(y)dy = 0 for 1 ≤ m ≤ 2M -1. (5.7) Clearly, we have

  x -z) f (x -z) dz.(6.18)

Table 1

 1 The coefficients {r l } l=L-2 l=1 for Daubechies' wavelets, where L = 2M and M = 5, • • • , 9. If the integrals in (4.1) or (4.2) exist, then the coefficients r

		Coefficients			Coefficients
	l	r l		l	r l
	M = 5 1	-0.82590601185015	M = 8	1	-0.88344604609097
	2	0.22882018706694		2	0.30325935147672
	3	-5.3352571932672E-02		3	-0.10636406828947
	4	7.4613963657755E-03		4	3.1290147839488E-02
	5	-2.3923582002393E-04		5	-6.9583791164537E-03
	6	-5.4047301644748E-05		6	1.0315302133757E-03
	7	-2.5241171135682E-07		7	-7.6677069083796E-05
	8	-2.6960479423517E-10		8	-2.4519921109537E-07
				9	-3.9938104563894E-08
				10	7.2079482385949E-08
	M = 6 1	-0.85013666155592		11	9.6971849256415E-10
	2	0.25855294414146		12	7.2522069166503E-13
	3	-7.2440589997659E-02		13	-1.2400785360984E-14
	4	1.4545511041994E-02		14	1.5854647516841E-19
	5	-1.5885615434757E-03
	6	4.2968915709948E-06
	7	1.2026575195723E-05	M = 9	1	-0.89531640583699
	8	4.2069120451167E-07		2	0.32031206224855
	9	-2.8996668057051E-09		3	-0.12095364936000
	10	6.9686511520083E-13		4	3.9952721886694E-02
				5	-1.0616930669821E-02
				6	2.1034028106558E-03
	M = 7 1	-0.86874391452377		7	-2.7812077649932E-04
	2	0.28296509452594		8	1.9620437763642E-05
	3	-9.0189066217795E-02		9	-4.8782468879634E-07
	4	2.2687411014648E-02		10	1.0361220591478E-07
	5	-3.8814546576295E-03		11	-1.5966864798639E-08
	6	3.3734404776409E-04		12	-8.1374108294110E-10
	7	4.2363946800701E-06		13	-5.4025197533630E-13
	8	-1.6501679210868E-06		14	-4.7814005916812E-14
	9	-2.1871130331900E-07		15	-1.6187880013009E-18
	10	4.1830548203747E-10		16	-4.8507474310747E-24
	11	-1.2035273999989E-11
	12	-6.6283900594600E-16

Table 2

 2 Condition numbers of the matrix of periodized second derivative (with and without preconditioning) in the basis of Daubechies' wavelets with three vanishing moments M = 3.

	N	κ	κp
	64	0.14545E+04	0.10792E+02
	128	0.58181E+04	0.11511E+02
	256	0.23272E+05	0.12091E+02
	512	0.93089E+05	0.12604E+02
	1024	0.37236E+06	0.13045E+02

Table 3

 3 Condition numbers of the matrix of periodized second derivative (with and without preconditioning) in the basis of Daubechies' wavelets with six vanishing moments M = 6.

	N	κ	κp
	64	0.10472E+04	0.43542E+01
	128	0.41886E+04	0.43595E+01
	256	0.16754E+05	0.43620E+01
	512	0.67018E+05	0.43633E+01
	1024	0.26807E+06	0.43640E+01

Table 4

 4 The coefficients r l , l = 1, • • • , 16 of the Hilbert transform for Daubechies' wavelet with six vanishing moments.

			Coefficients		Coefficients
	l		r l	l	r l
	M = 6	1	-0.588303698	9	-0.035367761
		2	-0.077576414	10	-0.031830988
		3	-0.128743695	11	-0.028937262
		4	-0.075063628	12	-0.026525823
		5	-0.064168018	13	-0.024485376
		6	-0.053041366	14	-0.022736420
		7	-0.045470650	15	-0.021220659
		8	-0.039788641	16	-0.019894368

Table 5

 5 The coefficients {r l } l , l = -7, • • • , 14 of the fractional derivative α = 0.5 for Daubechies' wavelet with six vanishing moments.

			Coefficients		Coefficients
		l	r l	l	r l
	M = 6	-7 -2.82831017E-06	4	-2.77955293E-02
		-6 -1.68623867E-06	5	-2.61324170E-02
		-5	4.45847796E-04	6	-1.91718816E-02
		-4 -4.34633415E-03	7	-1.52272841E-02
		-3	2.28821728E-02	8	-1.24667403E-02
		-2 -8.49883759E-02	9	-1.04479500E-02
		-1	0.27799963	10	-8.92061945E-03
		0	0.84681966	11	-7.73225246E-03
		1	-0.69847577	12	-6.78614593E-03
		2	2.36400139E-02	13	-6.01838599E-03
		3	-8.97463780E-02	14	-5.38521459E-03

Table 6

 6 The coefficients {t L+2 for Daubechies' wavelet with three vanishing moments, where L = 6 and j = 1, • • • , 8.

		(j) l } l=L-2 l=-Coefficients		Coefficients
		l	t	(j) l	l	t	(j) l
	j = 1	-4	0.		j = 5 -4	-8.3516169979703E-06
		-3	0.		-3	-4.0407157939626E-04
		-2	1.171875E-02	-2	4.1333660119562E-03
		-1	-9.765625E-02	-1	-2.1698923046642E-02
		0	0.5859375	0	0.99752855458064
		1	0.5859375	1	2.4860978555807E-02
		2	-9.765625E-02	2	-4.9328931709169E-03
		3	1.171875E-02	3	5.0836550508393E-04
		4	0.		4	1.2974760466022E-05
	j = 2	-4	0.		j = 6 -4	-4.7352138210499E-06
		-3	-1.1444091796875E-03	-3	-2.1482413927743E-04
		-2	1.6403198242188E-02	-2	2.1652627381741E-03
		-1	-1.0258483886719E-01	-1	-1.1239479930566E-02
		0	0.87089538574219	0	0.99937113652686
		1	0.26206970214844	1	1.2046257104714E-02
		2	-5.1498413085938E-02	2	-2.3712690179423E-03
		3	5.7220458984375E-03	3	2.4169452359502E-04
		4	1.3732910156250E-04	4	5.9574082627023E-06
	j = 3	-4	-1.3411045074463E-05	j = 7 -4	-2.5174703821573E-06
		-3	-1.0904073715210E-03	-3	-1.1073373558501E-04
		-2	1.2418627738953E-02	-2	1.1081638044863E-03
		-1	-6.9901347160339E-02	-1	-5.7198034904338E-03
		0	0.96389651298523	0	0.99984123346637
		1	0.11541545391083	1	5.9237906308573E-03
		2	-2.3304820060730E-02	2	-1.1605296576369E-03
		3	2.5123357772827E-03	3	1.1756409462604E-04
		4	6.7055225372314E-05	4	2.8323576983791E-06