Summation Paths in Clenshaw-Curtis Quadrature

Two topics concerning the use of Clenshaw-Curtis quadrature within the Bayesian automatic adaptive quadrature approach to the numerical solution of Riemann integrals are considered. First, it is found that the efficient floating point computation of the coefficients of the Chebyshev series expansion of the integrand is to be done within a mathematical structure consisting of the union of coefficient families ordered into complete binary trees. Second, the scrutiny of the decay rates of the involved even and odd rank Chebyshev expansion coefficients with the increase of their rank labels enables the definition of Bayesian decision paths for the advancement to the numerical output.


Introduction
In this paper we address two topics concerning the use of Clenshaw-Curtis (CC) quadrature [1][2][3] within the Bayesian automatic adaptive quadrature (BAAQ) [4] to the numerical solution of Riemann integrals.The first topic deals with the efficient floating point computation of the Chebyshev expansion coefficients, while the second deals with the inference of Bayesian diagnostics addressing the integrand conditioning based on the knowledge of the computed set of the Chebyshev coefficients.
The derivation of CC quadrature sums of prescribed algebraic degree of accuracy n uses Chebyshev series expansions of the integrand function at quadrature knots spanned by the extremal points of the Chebyshev polynomial of first kind T n (y), y ∈ [−1 , 1].For the sake of consistency, essentials of the CC quadrature are summarized in section 2.
The most time consuming part of the CC quadrature sum evaluation stems from the computation of the coefficients of the involved Chebyshev series expansions of the integrand.Efficient algorithms for particular n values have been reported long time ago [2,5] by taking advantage of the fast Fourier transform (FFT) approach.The analysis done in section 3 unveils the occurrence of a general mathematical structure consisting of families of coefficients ordered into complete binary trees.
The reliability of a CC quadrature sum output heavily depends on the decay rates of the involved even and odd rank Chebyshev expansion coefficients with the increase of their rank labels.The scrutiny of this matter in section 4 allows the definition of Bayesian decision paths for the advancement toward the solution.
The paper ends with conclusions in section 5.

The Clenshaw-Curtis quadrature
The computation of the Riemann definite integral by CC quadrature interpolates the reduced integrand at the set of (n + 1) CC quadrature knots by the truncated Chebyshev series expansion where the double prime shows that the first and the last terms of the sum are halved.The derivation of the CC quadrature sum is done substituting f (x) in (1) by its approximation (4) and solving analytically the resulting integral.A weight-function-g(x)-dependent quadrature sum is obtained (see, e.g., [1-3, 5, 6]).
The most computer intensive task within the CC quadrature comes from the computation of the coefficients b n k in (4).The straightforward solution expresses the column vector B = (n/2)(b n 0 , • • • , b n n ) by the matrix product B = T Φ, where T is a symmetric matrix which collects together the values of the Chebyshev polynomials at the CC knots (3), T k j = T jk = T k (y n j ), k, j = 0, 1, • • • , n, while Φ is the column vector of the computed integrand values, Φ = (φ(y n 0 ), • • • , φ(y n n )) .Since each line of the matrix T runs over the set (3), the occurrence of vanishing and ±1 elements enables a T matrix irreducible block decomposition which results in a significantly reduced computational cost of the numerical evaluation of the coefficient set B. Its derivation heavily depends on the properties of the set (3) of the CC quadrature knots.
In this paper we have adopted the original convention [1] of labelling the CC quadrature knots in the increasing order of the arguments of the cosine function.As a consequence, the elements in the set (3) monotonically decrease with the label j from y n 0 = 1 to y n n = −1 and are symmetrically distributed around the origin of [−1, 1] : y n n− j = y n j .The CC quadrature knots are characterized by inheritance, y 2n 2 j = y n j .Moreover, under even n, the middle knot vanishes, y n n/2 = 0.These properties show that, for n = 2 m , which is assumed henceforth, there are 2 m−1 − 1 nontrivial knots ( 3) inside the open interval (0, 1) and these can be ordered within a knot complete binary tree of root y (2)  1 = y (m) 2 m−2 = √ 2/2 and depth m − 2. To simplify the notation, here and in what follows, (m) stays for 2 m .At the l-th depth level inside the tree there are 2 l knots, The binary tree and its heap ordering key at given m will be referred to as H m .The occurrence of special arguments in (5) allows the computation of all the elements of H m to machine accuracy using the square root function and the recurrence.Indeed, for the root tree we have (Note that the root is its own "genetically related sibling".)

Binary tree structures of the Chebyshev coefficients
Within the variable degree nonadaptive CC quadrature, the information acquired at the polynomial degree 2 m−1 (m ≥ 2) is preserved if this degree is doubled to 2 m .The easiest way of inheriting this information makes use of the recurrence relationships [3], which start with the initial conditions defined at the integer reduced knots, Irreducible expressions of the corrections C (m) q are established by complete induction, Here δ s+1 and σ s+1 define iterated odd respectively even linear combinations at the (s + 1)-th height level of H m , with the initial conditions coming from the newly added integrand values under polynomial degree doubling.The equation (10 which collects together the set of corrections generated at the (s + 1)-th height level of H m .Let Mathematical Modeling and Computational Physics 2015 Then the set of corrections (10) can be written in matrix form Therefore, to get the coefficients b (m) in (4), it is necessary and sufficient to do (4 m−1 − 1)/3 multiplications instead of those (2 m + 1) 2 required by the straightforward matrix multiplication B = T Φ.
The vanishing correction C (m) 2 m−1 = 0 which ends the recurrence chain in (9) results in the immediate consequence that the expressions of the Chebyshev coefficients within the subset are free of multiplications, These expressions serve as the roots of m + 2 distinct coefficient complete binary trees such that: 1 0 .The subsets of coefficients are leaves, at the depth levels m − (s + 1), of the m − 1 complete binary trees starting with the roots 2 0 .The remaining three coefficients B (m) 0 , B (m) 2 m , and B (m) 2 m−1 are root leaves in the corresponding binary trees.
The occurrence of several families of coefficient complete binary trees reduces the efficiency of the fast Chebyshev transform implementation of the CC quadrature as compared to the logarithmic decrease of the number of multiplications within the FFT algorithms.The reduction factor of the number of multiplications asymptotically tends to the constant value of 9.As a consequence, the most interesting CC quadrature algorithms may be written down at values m = 2, 3, 4, and 5 with characteristic reduction factors of 25, 13.5, 10.7, and 9.7, respectively.
For the BAAQ approach to the solution of the Riemann integrals the most important CC quadrature algorithm is obtained at the polynomial degree of n = 32 (m = 5).The corresponding algebraic degree of accuracy (d = 32) exceeds by one unit that of the Gauss-Kronrod 10-21 quadrature rule [2].
The property of the interpolatory polynomial (4) of being close to the minimax polynomial approximation provides a sensitive tool for the early detection of poorly conditioned integrands f (x).

Bayesian inferences from the Chebyshev coefficients
The computation of the Chebyshev expansion coefficients b (m)  k in (4) has been implemented in two codes: one using the forward recurrence (7)-(9) and one using the backward recurrence within the complete binary tree structure defined above.
The round-off errors coming from the double precision floating point computations of the integrand values got magnified, especially for the higher rank expansion coefficient, by the differences involved in the occurring linear combinations of the integrand values.Due to the different summation paths, the outputs of the two codes were not identical, except for the lower rank coefficient staying at or near the highest absolute magnitudes.Coefficient value denormalizations to the highest floating point encountered characteristics evidenced agreement within the corresponding coefficient pairs, except for the lowest significant RAM digits.
A more interesting feature evidenced by this scrutiny was the possibility to easily get a Bayesian diagnostic on the integrand conditioning over the current integration range pointing either to the occurrence of an easily solvable integral or of major difficulties the clarification of which needs separate (and immediate) integrand profile analysis [4].
(ii) The coefficient complete binary tree structure occurring at a given m ≥ 2 inherits the complete binary tree structure defined at the Chebyshev polynomial degree 2 m−1 .This recursion property of the families of complete binary trees of the Chebyshev coefficients enables the reduction of the total number of multiplications asked for their computation from (2 m + 1) 2 to (4 m − 3m − 1)/9.
(iii) As a consequence of the parity properties of the Chebyshev polynomials of the first kind, there are two independent of each other subsets of the Chebyshev expansion coefficients: the subset of 2 m−1 odd label coefficients and the subset of 2 m−1 + 1 even label coefficients which carry at output the distinct information coming from the odd and even parts of the integrand respectively.
(iv) The comparison of the outputs of the two practical implementations of the computation of the Chebyshev expansion coefficients, done at n = 2 5 = 32, unveiled two main features.
First, if the floating point approximations of the computed Chebyshev series expansion coefficients are denormalized to the highest characteristics occurring in the union of the even and odd rank coefficients, then the forward and backward code outputs were different from each other by no more than the least significant RAM digits.
Second, the reliability of the emerging CC quadrature sum output is directly related to the maximum denormalizing value observed.A low value of this quantity points to a poorly conditioned integrand over the current integration range (hence to the need to proceed immediately to the assessment of the difficult spots inside the integrand profile), while a large value points to a well conditioned integrand skipping the need of the integrand profile analysis.
each pair of "genetically related descendants", EPJ Web of Conferences of arguments α and π/2 − α respectively, are obtained from the values of cos 2α (the reference father of argument 2α) and of sin 2α (its "genetically related sibling" of argument π/2 − 2α), as follows:

Figure 1 .
Figure 1.Patterns of variation of the absolute magnitudes of the Chebyshev expansion coefficients within the even and odd rank subsets versus the coefficient labels.The data on the left figure were derived for the function f 1 (x) over the specified subranges.The data on the right figure were derived for the function f 2 (x) at p = 5 (marked as 'p05' in the file names), at fixed x 0 = −1 (not marked), and at the specified four ω values.The file notations start with the specification of the rank of the Chebyshev subset: 'e' (for even) and 'o' (for odd).Three typical integrand conditioning diagnostics are apparent: (1) Cases (a): well-conditioned, fast converging.(2) Cases (b): well-conditioned, hopefully converging.(3) Cases (c) and (d): ill-conditioned -integrand profile analysis requested to set precise diagnostic.