## 5 Numeric Integration
### Integrating on a Double-exponential Grid  
Popular approaches to calculate numeric Fourier transforms include straightforward fast Fourier transform, and [PERSON]'s simple \"Filon-trapezoidal\" rule [35]. Both methods evaluate the Fourier integrand on an equidistant grid \(t_{k}=k\Delta t\). The Filon rule optimizes the weight of the grid points.  
In our application, especially for small \(\beta\), the decay of \(g(t)\) extends over several decades. To limit the number of grid points that must be taken into account, it is customary to use a decimation algorithm. A more efficient and perhaps even simpler alternative is the double-exponential transformation. It was first proposed by [PERSON] and [PERSON] in 1974 for the efficient evaluation of integrals with end-point  
Figure 2: Frequency limits for the series expansions for a required accuracy \(\delta=2\times 10^{-16}\) and a machine precision \(\epsilon=1\times 10^{-19}\). Orange points represent \(\omega_{\rm L}\), green points \(\omega_{\rm H}\). Black lines show the piecewise fits that are hardcoded in libkww to decide whether a series expansion is tried or whether numeric integration is used from the outset.  
singularities [36, 37]. Afterwards, it was adapted to oscillatory functions by [PERSON] and [PERSON] [38, 39]. The key idea is to choose grid points \(t_{k}\) close to the zeros of the sine function in Equation (38).  
A double-exponential transformation is a monotonous function \(\phi(x)\) that satisfies  
\[\phi(x\to-\infty)\to 0 \tag{41}\] \[\phi^{\prime}(x\to-\infty)\to 0\text{ double exponentially}\] (42) \[\phi(x\to+\infty)\to x\text{ double exponentially} \tag{43}\]  
This transformation shall now be applied to the time variable in the Fourier integral Equation (38):  
\[t=\frac{\pi}{\omega}\phi(k-\
u) \tag{44}\]  
yielding  
\[Y=\frac{\pi}{\omega}\int_{-\infty}^{\infty}\!\!\!\mathrm{d}k\,\phi^{\prime}(k- \
u)\sin(\pi(\phi(k-\
u)+\
u))g\left(\frac{\pi\phi(k-\
u)}{\omega}\right) \tag{45}\]  
It is convenient to abbreviate  
\[a_{k} := \pi\phi(k-\
u) \tag{46}\] \[b_{k} := \phi^{\prime}(k-\
u)\sin(\pi(\phi(k-\
u)+\
u))\] (47) \[\tilde{Y} := \frac{\omega}{\pi}Y \tag{48}\]  
The integral Equation (45) shall now be approximated by a sum, using the trapezoidal rule with stepwidth 1:  
\[\tilde{Y}=\sum_{k=-\infty}^{+\infty}b_{k}\,g\Big{(}\frac{a_{k}}{\omega}\Big{)} +\Delta_{\mathrm{di}}\tilde{Y} \tag{49}\]  
where the last term is the discretization error, to be discussed below (Section 5.4). As a second approximation, we truncate the summation at \(\pm N\),  
\[\tilde{Y}=\sum_{k=-N}^{+N}b_{k}\,g\Big{(}\frac{a_{k}}{\omega}\Big{)}+\Delta_{ \mathrm{di}}\tilde{Y}+\Delta_{\mathrm{tr}}\tilde{Y} \tag{50}\]  
where the new term is the truncation error, also discussed below. This sum is used in libkww to compute the KWW function at intermediate frequencies. Since \(a_{k}\) and \(b_{k}\) do not depend on \(\beta\) and \(\omega\), they must be generated only once, which greatly accelerates repeated evaluations of Equation (50).  
In practice, (50) can be well approximated with relatively small \(N\). For \(k\to-\infty\), condition (42) ensures that \(b_{k}\) goes double exponentially to \(0\). For \(k\to+\infty\), condition (43) makes the argument of the sine function in Equation (47) tend towards \(\pi k\). If \(N\) is integer, then all \(k\) are integer as well, and the sine can be expanded around \(\sin(\pi k)=0\). In consequence, \(|b_{k}|\) goes double exponentially to \(0\).