Algorithm 916: Faddeyeva & Voigt Evaluation

• Algorithm 916 is a computational routine for evaluating the Faddeyeva function and its Voigt components with high precision (up to 13–15 digits) using series expansions and asymptotic methods.
• It balances accuracy and efficiency through tunable parameters, precomputed exponentials, and domain-specific shortcuts, and is implemented in MATLAB, Fortran, and Scilab.
• Benchmarks show that Algorithm 916 outperforms earlier methods in speed and precision, making it essential for plasma physics, radiative transfer, and high-throughput spectral synthesis.

Algorithm 916 is a robust, highly accurate, and tunable computational algorithm for evaluating the Faddeyeva function $w(z)$ and its real (Voigt function $V(x,y)$) and imaginary ($L(x,y)$) components. Originally implemented as a MATLAB routine, it was designed to achieve both high precision (up to 13–15 significant digits) and efficient computation across broad domains in the complex $z = x + iy$ plane, particularly for applications in plasma physics, radiative transfer, and spectroscopy. Algorithm 916 introduced systematic control over the accuracy vs. efficiency trade-off and has subsequently been optimized, extended to Fortran and Scilab, and incorporated into high-throughput frameworks such as GPU-based HELIOS-K.

1. Mathematical Foundations

The Faddeyeva function is defined as

$w(z) = e^{-z^2} \operatorname{erfc}(-iz)$

where $z = x + iy$ and $\operatorname{erfc}$ denotes the complementary error function. An equivalent form is

$$w(z) = e^{-z^2}[1 + \operatorname{erf}(iz)] = V(x, y) + i L(x, y)\quad (y > 0),$$

with $V(x, y)$ and $L(x, y)$ denoting the real and imaginary Voigt functions. The Voigt profile, common in spectroscopy and astrophysics, is the real part of $w(z)$ and is expressed as

$$V(x, y) = \frac{y}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{(x-t)^2 + y^2}\, dt.$$

For efficient computation, the algorithm employs the Salzer expansion for $e^{-t^2}$: $$e^{-t^2} \approx \frac{a}{\sqrt{\pi}} \sum_{n=-\infty}^{\infty} e^{-a^2 n^2} \cosh(2 a n t),$$ where $0 < a \leq 1$ tunes the relative error $E_{\rm rel} \sim 2 e^{-\pi^2/a^2}$.

2. Series Expansions and Integral Representations

Algorithm 916 rewrites the expressions for $\Re[w(z)]$ and $\Im[w(z)]$ as organized series using the above expansion and integral representations for the error function. The key reformulation involves expressing $w(z)$ as a sum of five one-dimensional series ($\Sigma_1 \ldots \Sigma_5$), avoiding direct computation of double sums and leveraging symmetries such as

$$w(-z) = 2 e^{-z^2} - w(z), \qquad \Re[w(-x, y)] = \Re[w(x, y)], \qquad \Im[w(-x, y)] = -\Im[w(x, y)].$$

This rearrangement is critical for uniform convergence and stable evaluation even near problematic regions of the complex plane.

For large $|z|$, Algorithm 916 switches to an asymptotic continued-fraction expansion (after Faddeeva–Terent’ev and Gautschi) to preserve both speed and accuracy: $$w(z) \approx \frac{i}{\sqrt{\pi}} \frac{1}{z - \frac{1/2}{z - \frac{1}{z - \frac{3/2}{z}}}}$$ which ensures relative error $< 10^{-13}$ for $|z| \geq 3.8 \times 10^4$ (Zaghloul, 2015).

3. Algorithmic Implementation and Parameter Control

The main MATLAB function signature is

w = faddeyeva(z) w = faddeyeva(z, tiny)

where z is a complex array and tiny is the user-specified relative error threshold, defaulting to $\sim 1.4 \times 10^{-17}$ (machine precision in IEEE double-precision). Internally, tiny is mapped to the expansion parameter $a$ via

$$\text{tiny} \approx 2 e^{-\pi^2 / a^2}, \quad 0.5 \leq a \leq 1,$$

thus directly linking accuracy and convergence rate.

For performance, the algorithm:

• Precomputes exponentials and updates them iteratively within the summation loop.
• Makes domain-specific shortcuts: direct use of erfcx for purely imaginary arguments, rapid switching to continued fractions for large $|z|$, and logic to avoid catastrophic cancellation near the real axis.
• Truncates the series when new terms fall below tiny or machine epsilon.

A pseudocode kernel captures this structure:

function w = faddeyeva(z, tiny) % initialization... for each (x, y > 0) in z % handle special cases for n = 1 to n_max % update terms for Sigma1...Sigma5 if (new_term/acc_sum < max(tiny, eps)), break, end end % combine Sigma's w(x, y) = Re + i*Im % apply symmetry if necessary end end

(Zaghloul et al., 2011, Zaghloul, 2015)

4. Accuracy vs Efficiency: Tunable Trade-Offs

Algorithm 916's distinctive feature is its tunable accuracy/speed trade-off using the tiny (or sdgts in later versions) parameter. Users can select:

Desired Accuracy (digits)	Input Parameter	Algorithmic Consequence	Typical Runtime (MATLAB, 2.8M points)
13–15	tiny ≈ 1e−17, sdgts=13	Full cycles, tight cutoff	1.55 s (v2)
~10	tiny = 1e−8, sdgts=8	Fewer cycles, earlier cutoff	0.76 s
~6	tiny = 1e−4, sdgts=4	Minimal cycles, w4 fallback	0.49 s

Forreduced precision ($\sim 10^{-4}$), the algorithm invokes a reformed version of Humlíček’s “w4” routine, correcting previous accuracy failures on the real axis (Zaghloul, 2015).

5. Benchmarks and Comparative Performance

Extensive benchmarks have demonstrated Algorithm 916's superiority across both accuracy and speed. For 2.8 million evaluations over $x \in [-200,200]$, $y \in [10^{-20},10^4]$:

• Full-precision (tiny ≈ 1.4e−17): Relative error $\sim 10^{-13}$, runtime 1.55 s (v2), versus 4.46 s (original), 107.4 s for Poppe & Wijers (Algorithm 680), 23.2 s for Humlíček (Zaghloul et al., 2011, Zaghloul, 2015).
• Low-accuracy mode (sdgts=4): Runtime 0.49 s, still avoiding the pathological failures of prior algorithms (e.g., negative Voigt values or total loss of precision near the real axis).

In all tested domains, Algorithm 916 and its v2 implementation are free of catastrophic inaccuracies documented for Hui et al., Humlíček, Weideman, and Letchworth & Benner.

6. Extensions, Implementations, and High-Throughput Applications

Subsequent work provided major efficiency improvements, fixed-cycle summation (removing dynamic convergence checks), precomputation of all necessary exponentials, logical domain partitioning, and optimized asymptotic thresholds. The “v2” codebase is available in MATLAB, Scilab, and Fortran 2008. The Fortran module supports both single and double precision, exposes optional partial derivatives, and handles large-scale evaluations with fine-grained control over significant digits via the sdgts argument (range 4–13).

Algorithm 916 is foundational in the HELIOS-K GPU-accelerated opacity calculator (Grimm et al., 2015). In this context, it is fused with Gauss–Hermite quadrature for the efficient batch evaluation of Voigt profiles on modern CUDA hardware. The regional switching (Algorithm 916 for $r^2 = a^2 + u^2 < 100$, GH3 or GH1 quadrature for larger $r^2$) ensures both optimal efficiency and uniform high-precision for exoplanet and radiative transfer models.

Code / Mode	Max Rel. Error	MATLAB Time (s, 2.8M pts)
Faddeyeva(z,13) (v2)	~1e-13	1.55
Faddeyeva(z,4) (v2)	~1e-4	0.49
Poppe & Wijers [1990]		107.4
Humlíček 1982		~10–23

(Zaghloul, 2015)

7. Practical Usage and Recommendations

For MATLAB/Scilab, the interface allows flexible selection of precision:

w = faddeyeva(z) % maximum accuracy (13–15 digits) w = faddeyeva(z,1e-8) % medium precision (~10 digits) w = faddeyeva(z,1e-4) % fast, low precision (~6 digits)

(Zaghloul et al., 2011)

For Fortran:

call Faddeyeva_v2_rk(z, sdgts, w, dVdx, dVdy, Stat)

with optional computation of partial derivatives.

Recommended settings:

• Use sdgts=13 for high-accuracy spectral synthesis, radiative transfer, or reference calculations (Zaghloul, 2015).
• For high-throughput environments (e.g., opacities, atmospheric retrieval): sdgts≈8 (10−8) balances speed and accuracy.
• For rapid, low-precision tasks where $10^{-4}$ error is acceptable: sdgts=4 activates the safe, compensated low-precision branch.
• Selection of parameters can be guided by published timing and error tables to optimize performance for domain-specific needs.

Algorithm 916's design, error analysis, and performance data position it as the state-of-the-art routine for practical and research-grade evaluation of the Faddeyeva and Voigt functions in scientific computing (Zaghloul et al., 2011, Grimm et al., 2015, Zaghloul, 2015).