## 1 Introduction  
The C library libkww (please see the Supplementary Material: libkww source archive) provides functions for computing the Laplace-Fourier transform of the stretched or compressed exponential function \(\exp(-t^{\beta})\). It improves upon the previous work [1, 2] in several respects: (1) A wider \(\beta\) range is covered; (2) Results have the full accuracy of double-precision floating-point numbers; (3) The computation is very fast, thanks to two measures: (a) rigorous error bounds allow to maximally extend the low and high frequency domains where series expansions are used, and (b) the numeric integration at intermediate frequency is enormously accelerated by using a recent mathematical innovation, the Ooura-Mori double exponential transform; (4) The implementation is made available in the most portable way, namely as a C library.  
Claims (1)-(3) require some explanation. [PERSON] _et al._ published tables for values of \(\beta\) between 0.01 and 2 [1]. However, for \(\beta\lesssim 0.3\), those tables only cover an asymptotic power-law regime, which renders them practically useless. This will become clear in Section 4.3.  
With respect to accuracy, one might argue that a numeric precision of \(10^{-3}\) or \(10^{-4}\) is largely sufficient for fitting spectroscopic data. However, violations of monotonicity at a level \(\delta\) can trap a fit algorithm in a haphazard local minimum unless the minimum search step of the algorithm is correspondingly set to \(\mathcal{O}(\delta)\). Since libkww returns function values with double precision, it can be smoothly integrated into existing fit routines.  
With respect to speed of calculation, one might argue that given today's computing power this is no longer a serious concern. However, the Fourier transform of the stretched exponential is often embedded in a convolution of a theoretical model with an instrumental resolution function, which in turn is embedded in a nonlinear curve fitting routine. In such a situation, accelerating the innermost loop is still advantageous.  
The implementation in libkww is targeting a IEEE 754 compliant floating-point unit. Returned function values are accurate within the relative error of the double data type,  
\[\delta:=2.2\cdot 10^{-16}\lesssim 2^{-52} \tag{1}\]  
Internally, series expansions and trapezoid sums are computed using long double variables, expecting that this translates at least to 80 bits (extended double) so that floating-point errors are not larger than  
\[\epsilon:=5.5\cdot 10^{-20}\gtrsim 2^{-64} \tag{2}\]  
The modest signal-to-noise ratio \(\delta/\epsilon\) made painstaking fine-tuning of the numeric integration unavoidable. The good side is that it will be easy to port libkww to other architectures.  
Figure 1: Functions \(Q_{\beta}(\omega)\), \(V_{\beta}(\omega)\), \(P_{\beta}(\omega)\) for selected values of \(\beta\). Colored lines have been obtained by series expansions, solid black lines by numeric integration. Dashed lines have been computed using analytic expressions for \(\beta=1,2\) (Appendix E).  
This paper, after discussing typical applications (Section 2) and introducing some notation (Section 3), describes the mathematical foundations of the implemented algorithm in full detail (Sections 4 and 5). Based on this, the application programming interface and some test programs are documented (Section 6). Exemplary results are shown in Figure 1.