# NAG FL InterfaceC06 (Sum)Summation of Series

Settings help

FL Name Style:

FL Specification Language:

## 1Scope of the Chapter

This chapter is concerned with the following tasks:
1. (a)calculating the discrete Fourier transform of a sequence of real or complex data values;
2. (b)calculating the discrete convolution or the discrete correlation of two sequences of real or complex data values using discrete Fourier transforms;
3. (c)calculating the inverse Laplace transform of a user-supplied subroutine;
4. (d)calculating the fast Gauss transform approximation to the discrete Gauss transform;
5. (e)direct summation of orthogonal series;
6. (f)acceleration of convergence of a sequence of real values.

## 2Background to the Problems

### 2.1Discrete Fourier Transforms

#### 2.1.1Complex transforms

Most of the routines in this chapter calculate the finite discrete Fourier transform (DFT) of a sequence of $n$ complex numbers ${z}_{\mathit{j}}$, for $\mathit{j}=0,1,\dots ,n-1$. The direct transform is defined by
 $z^k = 1n ∑ j=0 n-1 zj exp(-i 2πjk n )$ (1)
for $k=0,1,\dots ,n-1$. Note that equation (1) makes sense for all integral $k$ and with this extension ${\stackrel{^}{z}}_{k}$ is periodic with period $n$, i.e., ${\stackrel{^}{z}}_{k}={\stackrel{^}{z}}_{k±n}$, and in particular ${\stackrel{^}{z}}_{-k}={\stackrel{^}{z}}_{n-k}$. Note also that the scale-factor of $\frac{1}{\sqrt{n}}$ may be omitted in the definition of the DFT, and replaced by $\frac{1}{n}$ in the definition of the inverse.
If we write ${z}_{j}={x}_{j}+i{y}_{j}$ and ${\stackrel{^}{z}}_{k}={a}_{k}+i{b}_{k}$, then the definition of ${\stackrel{^}{z}}_{k}$ may be written in terms of sines and cosines as
 $ak = 1n ∑ j=0 n-1 (xjcos( 2πjk n )+yjsin( 2πjk n ))$
 $bk = 1n ∑ j= 0 n- 1 (yjcos( 2πjk n )-xjsin( 2πjk n )) .$
The original data values ${z}_{j}$ may conversely be recovered from the transform ${\stackrel{^}{z}}_{k}$ by an inverse discrete Fourier transform:
 $zj = 1n ∑ k=0 n-1 z^k exp(+i 2πjk n )$ (2)
for $j=0,1,\dots ,n-1$. If we take the complex conjugate of (2), we find that the sequence ${\overline{z}}_{j}$ is the DFT of the sequence ${\overline{\stackrel{^}{z}}}_{k}$. Hence the inverse DFT of the sequence ${\stackrel{^}{z}}_{k}$ may be obtained by taking the complex conjugates of the ${\stackrel{^}{z}}_{k}$; performing a DFT, and taking the complex conjugates of the result. (Note that the terms forward transform and backward transform are also used to mean the direct and inverse transforms respectively.)
The definition (1) of a one-dimensional transform can easily be extended to multidimensional transforms. For example, in two dimensions we have
 $z^ k1 k2 = 1 n1 n2 ∑ j1=0 n1-1 ∑ j2=0 n2-1 z j1 j2 exp(-i 2 π j1 k1 n1 ) exp(-i 2 π j2 k2 n2 ) .$ (3)
Note:  definitions of the discrete Fourier transform vary. Sometimes (2) is used as the definition of the DFT, and (1) as the definition of the inverse.

#### 2.1.2Real transforms

If the original sequence is purely real valued, i.e., ${z}_{j}={x}_{j}$, then
 $z^k = ak + i bk = 1n ∑ j=0 n-1 xj exp(-i 2πjk n )$
and ${\stackrel{^}{z}}_{n-k}$ is the complex conjugate of ${\stackrel{^}{z}}_{k}$. Thus the DFT of a real sequence is a particular type of complex sequence, called a Hermitian sequence, or half-complex or conjugate symmetric, with the properties
 $a n-k = ak b n-k = -bk b0 = 0$
and, if $n$ is even, ${b}_{n/2}=0$.
Thus a Hermitian sequence of $n$ complex data values can be represented by only $n$, rather than $2n$, independent real values. This can obviously lead to economies in storage, with two schemes being used in this chapter. In the first (deprecated) scheme, which will be referred to as the real storage format for Hermitian sequences, the real parts ${a}_{k}$ for $0\le k\le n/2$ are stored in normal order in the first $n/2+1$ locations of an array x of length $n$; the corresponding nonzero imaginary parts are stored in reverse order in the remaining locations of x. To clarify, if x is declared with bounds $\left(0:n-1\right)$ in your calling subroutine, the following two tables illustrate the storage of the real and imaginary parts of ${\stackrel{^}{z}}_{k}$ for the two cases: $n$ even and $n$ odd.
If $n$ is even then the sequence has two purely real elements and is stored as follows:
Index of x $0$ $1$ $2$ $\dots$ $\mathbit{n}/2$ $\dots$ $\mathbit{n}-2$ $\mathbit{n}-1$
Sequence ${a}_{0}$ ${a}_{1}+i{b}_{1}$ ${a}_{2}+i{b}_{2}$ $\dots$ ${a}_{n/2}$ $\dots$ ${a}_{2}-i{b}_{2}$ ${a}_{1}-i{b}_{1}$
Stored values ${a}_{0}$ ${a}_{1}$ ${a}_{2}$ $\dots$ ${a}_{n/2}$ $\dots$ ${b}_{2}$ ${b}_{1}$
 $x(k) = ak , for ​ k= 0, 1, …, n/2 , and x(n-k) = bk , for ​ k= 1, 2, …, n/2-1 .$
If $n$ is odd then the sequence has one purely real element and, letting $n=2s+1$, is stored as follows:
Index of x $0$ $1$ $2$ $\dots$ $\mathbit{s}$ $\mathbit{s}+1$ $\dots$ $\mathbit{n}-2$ $\mathbit{n}-1$
Sequence ${a}_{0}$ ${a}_{1}+i{b}_{1}$ ${a}_{2}+i{b}_{2}$ $\dots$ ${a}_{s}+i{b}_{s}$ ${a}_{s}-i{b}_{s}$ $\dots$ ${a}_{2}-i{b}_{2}$ ${a}_{1}-i{b}_{1}$
Stored values ${a}_{0}$ ${a}_{1}$ ${a}_{2}$ $\dots$ ${a}_{s}$ ${b}_{s}$ $\dots$ ${b}_{2}$ ${b}_{1}$
 $x(k) = ak , for ​ k= 0, 1, …, s , and x(n-k) = bk , for ​ k= 1, 2, …, s .$
The second (recommended) storage scheme, referred to in this chapter as the complex storage format for Hermitian sequences, stores the real and imaginary parts ${a}_{k},{b}_{k}$, for $0\le k\le n/2$, in consecutive locations of an array x of length $n+2$. If x is declared with bounds $\left(0:n+1\right)$ in your calling subroutine, the following two tables illustrate the storage of the real and imaginary parts of ${\stackrel{^}{z}}_{k}$ for the two cases: $n$ even and $n$ odd.
If $n$ is even then the sequence has two purely real elements and is stored as follows:
Index of x $0$ $1$ $2$ $3$ $\dots$ $\mathbit{n}-2$ $\mathbit{n}-1$ $\mathbit{n}$ $\mathbit{n}+1$
Stored values ${a}_{0}$ ${b}_{0}=0$ ${a}_{1}$ ${b}_{1}$ $\dots$ ${a}_{n/2-1}$ ${b}_{n/2-1}$ ${a}_{n/2}$ ${b}_{n/2}=0$
 $x(2×k) = ak , for ​ k= 0, 1, …, n/2 , and x(2×k+1) = bk , for ​ k= 0, 1, …, n/2 .$
If $n$ is odd then the sequence has one purely real element and, letting $n=2s+1$, is stored as follows:
Index of x $0$ $1$ $2$ $3$ $\dots$ $\mathbit{n}-2$ $\mathbit{n}-1$ $\mathbit{n}$ $\mathbit{n}+1$
Stored values ${a}_{0}$ ${b}_{0}=0$ ${a}_{1}$ ${b}_{1}$ $\dots$ ${b}_{s-1}$ ${a}_{s}$ ${b}_{s}$ $0$
 $x(2×k) = ak , for ​ k= 0, 1, …, s , and x(2×k+1) = bk , for ​ k= 0, 1, …, s .$
Also, given a Hermitian sequence, the inverse (or backward) discrete transform produces a real sequence. That is,
 $xj = 1n (a0+2 ∑ k=1 n/2-1 (akcos( 2πjk n )-bksin( 2πjk n ))+ a n/2 )$
where ${a}_{n/2}=0$ if $n$ is odd.
For real data that is two-dimensional or higher, the symmetry in the transform persists for the leading dimension only. So, using the notation of equation (3) for the complex two-dimensional discrete transform, we have that ${\stackrel{^}{z}}_{{k}_{1}{k}_{2}}$ is the complex conjugate of ${\stackrel{^}{z}}_{\left({n}_{1}-{k}_{1}\right)\left({n}_{2}-{k}_{2}\right)}$. It is more convenient for transformed data of two or more dimensions to be stored as a complex sequence of length $\left({n}_{1}/2+1\right)×{n}_{2}×\cdots ×{n}_{d}$ where $d$ is the number of dimensions. The inverse discrete Fourier transform operating on such a complex sequence (Hermitian in the leading dimension) returns a real array of full dimension (${n}_{1}×{n}_{2}×\cdots ×{n}_{d}$).

#### 2.1.3Real symmetric transforms

In many applications the sequence ${x}_{j}$ will not only be real, but may also possess additional symmetries which we may exploit to reduce further the computing time and storage requirements. For example, if the sequence ${x}_{j}$ is odd, $\left({x}_{j}={-x}_{n-j}\right)$, then the discrete Fourier transform of ${x}_{j}$ contains only sine terms. Rather than compute the transform of an odd sequence, we define the sine transform of a real sequence by
 $x^k = 2n ∑j=1 n-1 xj sin(πjkn) ,$
which could have been computed using the Fourier transform of a real odd sequence of length $2n$. In this case the ${x}_{j}$ are arbitrary, and the symmetry only becomes apparent when the sequence is extended. Similarly we define the cosine transform of a real sequence by
 $x^k = 2n (12x0+ ∑ j=1 n-1 xjcos(πjkn)+12(-1)kxn)$
which could have been computed using the Fourier transform of a real even sequence of length $2n$.
In addition to these ‘half-wave’ symmetries described above, sequences arise in practice with ‘quarter-wave’ symmetries. We define the quarter-wave sine transform by
 $x^k = 1n ( ∑ j=1 n-1 xjsin( π j (2k-1) 2n )+12(-1) k-1xn)$
which could have been computed using the Fourier transform of a real sequence of length $4n$ of the form
 $(0,x1,…,xn,xn-1 ,…,x1,0,-x1,…,-xn, -x n-1 ,…, -x 1 ) .$
Similarly we may define the quarter-wave cosine transform by
 $x^k = 1n (12x0+∑j= 1 n- 1xjcos( π j (2k-1) 2n ))$
which could have been computed using the Fourier transform of a real sequence of length $4n$ of the form
 $(x0,x1,…, x n-1 ,0, -x n-1 ,…,-x0,-x1,…, -x n-1 ,0, x n-1 ,…,x1) .$

#### 2.1.4Fourier integral transforms

The usual application of the discrete Fourier transform is that of obtaining an approximation of the Fourier integral transform
 $F (s) = ∫ -∞ ∞ f (t) exp(-i2πst) dt$
when $f\left(t\right)$ is negligible outside some region $\left(0,c\right)$. Dividing the region into $n$ equal intervals we have
 $F (s) ≅ cn ∑ j=0 n-1 fj exp( -i 2 π s j c n )$
and so
 $Fk ≅ cn ∑ j= 0 n- 1 fj exp( -i 2 π jk n )$
for $k=0,1,\dots ,n-1$, where ${f}_{j}=f\left(jc/n\right)$ and ${F}_{k}=F\left(k/c\right)$.
Hence the discrete Fourier transform gives an approximation to the Fourier integral transform in the region $s=0$ to $s=n/c$.
If the function $f\left(t\right)$ is defined over some more general interval $\left(a,b\right)$, then the integral transform can still be approximated by the discrete transform provided a shift is applied to move the point $a$ to the origin.

#### 2.1.5Convolutions and correlations

One of the most important applications of the discrete Fourier transform is to the computation of the discrete convolution or correlation of two vectors $x$ and $y$ defined (as in Brigham (1974)) by
• convolution: ${z}_{k}=\sum _{j=0}^{n-1}{x}_{j}{y}_{k-j}$
• correlation: ${w}_{k}=\sum _{j=0}^{n-1}{\overline{x}}_{j}{y}_{k+j}$
(Here $x$ and $y$ are assumed to be periodic with period $n$.)
Under certain circumstances (see Brigham (1974)) these can be used as approximations to the convolution or correlation integrals defined by
 $z (s) = ∫ -∞ ∞ x (t) y (s-t) dt$
and
 $w (s) = ∫ -∞ ∞ x¯ (t) y (s+t) dt , -∞ < s < ∞ .$
For more general advice on the use of Fourier transforms, see Hamming (1962); more detailed information on the fast Fourier transform algorithm can be found in Gentleman and Sande (1966) and Brigham (1974).

#### 2.1.6Applications to solving partial differential equations (PDEs)

A further application of the fast Fourier transform, and in particular of the Fourier transforms of symmetric sequences, is in the solution of elliptic PDEs. If an equation is discretized using finite differences, then it is possible to reduce the problem of solving the resulting large system of linear equations to that of solving a number of tridiagonal systems of linear equations. This is accomplished by uncoupling the equations using Fourier transforms, where the nature of the boundary conditions determines the choice of transforms – see Section 3.3. Full details of the Fourier method for the solution of PDEs may be found in Swarztrauber (1977) and Swarztrauber (1984).

### 2.2Inverse Laplace Transforms

Let $f\left(t\right)$ be a real function of $t$, with $f\left(t\right)=0$ for $t<0$, and be piecewise continuous and of exponential order $\alpha$, i.e.,
 $|f(t)| ≤ M eαt$
for large $t$, where $\alpha$ is the minimal such exponent.
The Laplace transform of $f\left(t\right)$ is given by
 $F (s) = ∫0∞ e-st f (t) dt , t>0$
where $F\left(s\right)$ is defined for $\mathrm{Re}\left(s\right)>\alpha$.
The inverse transform is defined by the Bromwich integral
 $f (t) = 12πi ∫ a-i∞ a+i∞ est F (s) ds , t>0 .$
The integration is performed along the line $s=a$ in the complex plane, where $a>\alpha$. This is equivalent to saying that the line $s=a$ lies to the right of all singularities of $F\left(s\right)$. For this reason, the value of $\alpha$ is crucial to the correct evaluation of the inverse. It is not essential to know $\alpha$ exactly, but an upper bound must be known.
The problem of determining an inverse Laplace transform may be classified according to whether (a) $F\left(s\right)$ is known for real values only, or (b) $F\left(s\right)$ is known in functional form and can, therefore, be calculated for complex values of $s$. Problem (a) is very ill-defined and no routines are provided. Two methods are provided for problem (b).

### 2.3Fast Gauss Transform

Gauss transforms have applications in areas including statistics, machine learning, and numerical solution of the heat equation. The discrete Gauss transform (DGT), $G\left(y\right)$, evaluated at a set of target points $y\left(j\right)$, for $j=1,2,\dots ,m\in {ℝ}^{d}$, is defined as:
 $G (yj) = ∑ i=1 n qi e - ‖yj-xi‖ 2 2 / h i 2 , j=1,…,m$
where ${x}_{i}$, for $i=1,2,\dots ,n\in {ℝ}^{d}$, are the Gaussian source points, ${q}_{i}$, for $i=1,2,\dots ,n\in {ℝ}^{+}$, are the source weights and ${h}_{i}$, for $i=1,2,\dots ,n\in {ℝ}^{+}$, are the source standard deviations (alternatively source scales or source bandwidths).
The fast Gauss transform (FGT) algorithm presented in Raykar and Duraiswami (2005) approximates the DGT by using two Taylor series and clustering of the source points.

### 2.4Direct Summation of Orthogonal Series

For any series of functions ${\varphi }_{i}$ which satisfy a recurrence
 $ϕr+1 (x) + αr (x) ϕr (x) + βr (x) ϕr-1 (x) =0$
the sum
 $∑ r= 0 n ar ϕr (x)$
is given by
 $∑ r=0 n ar ϕr (x) = b0 (x) ϕ0 (x) + b1 (x) (ϕ1(x)+α0(x)ϕ0(x))$
where
 $br (x) + αr (x) br+ 1 (x) + βr+ 1 (x) br+ 2 (x) = ar b n+ 1 (x) = b n+ 2 (x) = 0 .$
This may be used to compute the sum of the series. For further reading, see Hamming (1962).

### 2.5Acceleration of Convergence

This device has applications in a large number of fields, such as summation of series, calculation of integrals with oscillatory integrands (including, for example, Hankel transforms), and root-finding. The mathematical description is as follows. Given a sequence of values $\left\{{s}_{n}\right\}$, for $\mathit{n}=m,\dots ,m+2l$, then, except in certain singular cases, arguments, $a$, ${b}_{i}$, ${c}_{i}$ may be determined such that
 $sn = a + ∑ i=1 l bi cin .$
If the sequence $\left\{{s}_{n}\right\}$ converges, then $a$ may be taken as an estimate of the limit. The method will also find a pseudo-limit of certain divergent sequences – see Shanks (1955) for details.
To use the method to sum a series, the terms ${s}_{n}$ of the sequence should be the partial sums of the series, e.g., ${s}_{n}=\sum _{k=1}^{n}{t}_{k}$, where ${t}_{k}$ is the $k$th term of the series. The algorithm can also be used to some advantage to evaluate integrals with oscillatory integrands; one approach is to write the integral (in this case over a semi-infinite interval) as
 $∫0∞ f (x) dx = ∫0a1 f (x) dx + ∫ a1 a2 f (x) dx + ∫ a2 a3 f (x) dx + ⋯$
and to consider the sequence of values
 $s1 = ∫ 0 a1 f (x) dx , s2 = ∫ 0 a2 f (x) dx = s1 + ∫ a1 a2 f (x) dx , etc.,$
where the integrals are evaluated using standard quadrature methods. In choosing the values of the ${a}_{k}$, it is worth bearing in mind that c06baf converges much more rapidly for sequences whose values oscillate about a limit. The ${a}_{k}$ should thus be chosen to be (close to) the zeros of $f\left(x\right)$, so that successive contributions to the integral are of opposite sign. As an example, consider the case where $f\left(x\right)=M\left(x\right)\mathrm{sin}x$ and $M\left(x\right)>0$: convergence will be much improved if ${a}_{k}=k\pi$ rather than ${a}_{k}=2k\pi$.

## 3Recommendations on Choice and Use of Available Routines

The fast Fourier transform algorithm ceases to be ‘fast’ if applied to values of $n$ which cannot be expressed as a product of small prime factors. All the FFT routines in this chapter are particularly efficient if the only prime factors of $n$ are $2$, $3$ or $5$.

### 3.1One-dimensional Fourier Transforms

The choice of routine is determined first of all by whether the data values constitute a real, Hermitian or general complex sequence. It is wasteful of time and storage to use an inappropriate routine.

#### 3.1.1Real and Hermitian data

c06paf transforms a single sequence of real data onto (and in-place) a representation of the transformed Hermitian sequence using the complex storage scheme described in Section 2.1.2. c06paf also performs the inverse transform using the representation of Hermitian data and transforming back to a real data sequence.
Alternatively, the two-dimensional routine c06pvf can be used (on setting the second dimension to 1) to transform a sequence of real data onto an Hermitian sequence whose first half is stored in a separate Complex array. The second half need not be stored since these are the complex conjugate of the first half in reverse order. c06pwf performs the inverse operation, transforming the the Hermitian sequence (half-)stored in a Complex array onto a separate real array.
Multiple sequences of real data can also be transformed to and from a transformed Hermitian sequence (using the complex storage scheme). c06ppf assumes that the original real data is stored such that all the corresponding ($j$-th say) elements of the different streams are stored contiguously. c06pqf assumes that the elements of each sequence are stored contiguously; this is the recommended form of storage since it normally results in more efficient computation.

#### 3.1.2Complex data

c06pcf transforms a single complex sequence in-place; it also performs the inverse transform. c06psf transforms multiple complex sequences, each stored sequentially; it also performs the inverse transform on multiple complex sequences. This routine is designed to perform several transforms in a single call, all with the same value of $n$.
If extensive use is to be made of these routines and you are concerned about efficiency, you are advised to conduct your own timing tests.

### 3.2Half- and Quarter-wave Transforms

Four routines are provided for computing fast Fourier transforms (FFTs) of real symmetric sequences. c06ref computes multiple Fourier sine transforms, c06rff computes multiple Fourier cosine transforms, c06rgf computes multiple quarter-wave Fourier sine transforms, and c06rhf computes multiple quarter-wave Fourier cosine transforms.

### 3.3Application to Elliptic Partial Differential Equations

As described in Section 2.1.6, Fourier transforms may be used in the solution of elliptic PDEs.
c06ref may be used to solve equations where the solution is specified along the boundary.
c06rff may be used to solve equations where the derivative of the solution is specified along the boundary.
c06rgf may be used to solve equations where the solution is specified on the lower boundary, and the derivative of the solution is specified on the upper boundary.
c06rhf may be used to solve equations where the derivative of the solution is specified on the lower boundary, and the solution is specified on the upper boundary.
For equations with periodic boundary conditions the full-range Fourier transforms computed by c06paf are appropriate.

### 3.4Multidimensional Fourier Transforms

The following routines compute multidimensional discrete Fourier transforms of real, Hermitian and complex data stored in complex arrays:
real Hermitian complex
$2$ dimensions c06pvf c06pwf c06puf
$3$ dimensions c06pyf c06pzf c06pxf
any number of dimensions c06pjf
The Hermitian data, either transformed from or being transformed to real data, is compacted (due to symmetry) along its first dimension when stored in Complex arrays; thus approximately half the full Hermitian data is stored.
c06puf and c06pxf should be used in preference to c06pjf for two- and three-dimensional transforms, as they are easier to use and are likely to be more efficient.
The transform of multidimensional real data is stored as a complex sequence that is Hermitian in its leading dimension. The inverse transform takes such a complex sequence and computes the real transformed sequence. Consequently, separate routines are provided for performing forward and inverse transforms.
c06pvf performs the forward two-dimensionsal transform while c06pwf performs the inverse of this transform.
c06pyf performs the forward three-dimensional transform while c06pzf performs the inverse of this transform.
The complex sequences computed by c06pvf and c06pyf contain roughly half of the Fourier coefficients; the remainder can be reconstructed by conjugation of those computed. For example, the Fourier coefficients of the two-dimensional transform ${\stackrel{^}{z}}_{\left({n}_{1}-{k}_{1}\right){k}_{2}}$ are the complex conjugate of ${\stackrel{^}{z}}_{{k}_{1}{k}_{2}}$ for ${k}_{1}=0,1,\dots ,{n}_{1}/2$, and ${k}_{2}=0,1,\dots ,{n}_{2}-1$.

### 3.5Convolution and Correlation

c06fkf computes either the discrete convolution or the discrete correlation of two real vectors.
c06pkf computes either the discrete convolution or the discrete correlation of two complex vectors.

### 3.6Inverse Laplace Transforms

Two methods are provided: Weeks' method (c06lbf) and Crump's method (c06laf). Both require the function $F\left(s\right)$ to be evaluated for complex values of $s$. If in doubt which method to use, try Weeks' method (c06lbf) first; when it is suitable, it is usually much faster.
Typically the inversion of a Laplace transform becomes harder as $t$ increases so that all numerical methods tend to have a limit on the range of $t$ for which the inverse $f\left(t\right)$ can be computed. c06laf is useful for small and moderate values of $t$.
It is often convenient or necessary to scale a problem so that $\alpha$ is close to $0$. For this purpose it is useful to remember that the inverse of $F\left(s+k\right)$ is $\mathrm{exp}\left(-kt\right)f\left(t\right)$. The method used by c06laf is not so satisfactory when $f\left(t\right)$ is close to zero, in which case a term may be added to $F\left(s\right)$, e.g., $k/s+F\left(s\right)$ has the inverse $k+f\left(t\right)$.
Singularities in the inverse function $f\left(t\right)$ generally cause numerical methods to perform less well. The positions of singularities can often be identified by examination of $F\left(s\right)$. If $F\left(s\right)$ contains a term of the form $\mathrm{exp}\left(-ks\right)/s$ then a finite discontinuity may be expected in the inverse at $t=k$. c06laf, for example, is capable of estimating a discontinuous inverse but, as the approximation used is continuous, Gibbs' phenomena (overshoots around the discontinuity) result. If possible, such singularities of $F\left(s\right)$ should be removed before computing the inverse.

### 3.7Fast Gauss Transform

The only routine available is c06saf. If the dimensionality of the data is low or the number of source and target points is small, however, it may be more efficient to evaluate the discrete Gauss transform directly.

### 3.8Direct Summation of Orthogonal Series

The only routine available is c06dcf, which sums a finite Chebyshev series
 $∑ j=0 n cj Tj (x) , ∑ j=0 n cj T 2j (x) or ∑ j=0 n cj T 2j+1 (x)$
depending on the choice of argument.

### 3.9Acceleration of Convergence

The only routine available is c06baf.

## 4Decision Trees

### Tree 1: Fourier Transform of Discrete Complex Data

 Is the data one-dimensional? Multiple vectors? Stored as rows? c06prf yes yes yes no no no Stored as columns? c06psf yes c06pcf Is the data two-dimensional? c06puf yes no Is the data three-dimensional? c06pxf yes no Transform on one dimension only? c06pff yes no Transform on all dimensions? c06pjf yes

### Tree 2: Fourier Transform of Real Data or Data in Complex Hermitian Form Resulting from the Transform of Real Data

 Quarter-wave sine (inverse) transform? c06rgf yes no Quarter-wave cosine (inverse) transform? c06rhf yes no Sine (inverse) transform? c06ref yes no Cosine (inverse) transform? c06rff yes no Is the data three-dimensional? Forward transform on real data? c06pyf yes yes no no Inverse transform on Hermitian data? c06pzf yes Is the data two-dimensional? Forward transform on real data? c06pvf yes yes no no Inverse transform on Hermitian data? c06pwf yes Is the data multi one-dimensional? Sequence elements stored by stride (number of sequences)? c06ppf yes yes no no Sequence elements stored contiguously? c06pqf yes c06paf

## 5Functionality Index

 Acceleration of convergence c06baf
 Convolution or Correlation,
 complex vectors c06pkf
 real vectors,
 time-saving c06fkf
 Discrete Fourier Transform,
 multidimensional,
 complex sequence,
 complex storage c06pjf
 real storage c06fjf
 multiple half- and quarter-wave transforms,
 Fourier cosine transforms,
 simple use c06rbf
 Fourier cosine transforms, simple use c06rff
 Fourier sine transforms,
 simple use c06raf
 Fourier sine transforms, simple use c06ref
 quarter-wave cosine transforms,
 simple use c06rdf
 quarter-wave cosine transforms, simple use c06rhf
 quarter-wave sine transforms,
 simple use c06rcf
 quarter-wave sine transforms, simple use c06rgf
 one-dimensional,
 multiple transforms,
 complex sequence,
 complex storage, contiguous sequence elements c06psf
 complex storage by rows c06prf
 Hermitian/real sequence,
 complex storage, contiguous sequence elements c06pqf
 complex storage, strided sequence elements c06ppf
 multi-variable,
 complex sequence,
 complex storage c06pff
 real storage c06fff
 single transforms,
 complex sequence,
 time-saving,
 complex storage c06pcf
 real storage c06fcf
 Hermitian/real sequence,
 time-saving,
 complex storage c06paf
 Hermitian sequence,
 time-saving,
 real storage c06fbf
 real sequence,
 time-saving,
 real storage c06faf
 three-dimensional,
 complex sequence,
 complex storage c06pxf
 real storage c06fxf
 Hermitian/real sequence,
 complex-to-real c06pzf
 real-to-complex c06pyf
 two-dimensional,
 complex sequence,
 complex storage c06puf
 Hermitian/real sequence,
 complex-to-real c06pwf
 real-to-complex c06pvf
 Fast Gauss Transform c06saf
 Inverse Laplace Transform,
 Crump's method c06laf
 Weeks' method,
 compute coefficients of solution c06lbf
 evaluate solution c06lcf
 Summation of Chebyshev series c06dcf

None.

## 7 Withdrawn or Deprecated Routines

The following lists all those routines that have been withdrawn since Mark 23 of the Library or are in the Library, but deprecated.
Routine Status Replacement Routine(s)
c06dbf Withdrawn at Mark 25 c06dcf
c06eaf Withdrawn at Mark 26 c06paf
c06ebf Withdrawn at Mark 26 c06paf
c06ecf Withdrawn at Mark 26 c06pcf
c06ekf Withdrawn at Mark 26 c06fkf
c06fpf To be withdrawn at Mark 28 c06ppf or c06pqf
c06fqf To be withdrawn at Mark 28 c06ppf or c06pqf
c06frf Withdrawn at Mark 26 c06psf
c06fuf Withdrawn at Mark 26 c06puf
c06gbf Withdrawn at Mark 26 No replacement required
c06gcf Withdrawn at Mark 26 No replacement required
c06gqf Withdrawn at Mark 26 No replacement required
c06gsf Withdrawn at Mark 26 No replacement required
c06haf Withdrawn at Mark 26 c06ref
c06hbf Withdrawn at Mark 26 c06rff
c06hcf Withdrawn at Mark 26 c06rgf
c06hdf Withdrawn at Mark 26 c06rhf

## 8References

Brigham E O (1974) The Fast Fourier Transform Prentice–Hall
Davies S B and Martin B (1979) Numerical inversion of the Laplace transform: A survey and comparison of methods J. Comput. Phys. 33 1–32
Fox L and Parker I B (1968) Chebyshev Polynomials in Numerical Analysis Oxford University Press
Gentleman W S and Sande G (1966) Fast Fourier transforms for fun and profit Proc. Joint Computer Conference, AFIPS 29 563–578
Hamming R W (1962) Numerical Methods for Scientists and Engineers McGraw–Hill
Raykar V C and Duraiswami R (2005) Improved Fast Gauss Transform With Variable Source Scales University of Maryland Technical Report CS-TR-4727/UMIACS-TR-2005-34
Shanks D (1955) Nonlinear transformations of divergent and slowly convergent sequences J. Math. Phys. 34 1–42
Swarztrauber P N (1977) The methods of cyclic reduction, Fourier analysis and the FACR algorithm for the discrete solution of Poisson's equation on a rectangle SIAM Rev. 19(3) 490–501
Swarztrauber P N (1984) Fast Poisson solvers Studies in Numerical Analysis (ed G H Golub) Mathematical Association of America
Swarztrauber P N (1986) Symmetric FFT's Math. Comput. 47(175) 323–346
Wynn P (1956) On a device for computing the ${e}_{m}\left({S}_{n}\right)$ transformation Math. Tables Aids Comput. 10 91–96