---
title: 'Algorithm 916: Faddeyeva & Voigt Evaluation'
url: https://www.emergentmind.com/topics/algorithm-916
type: topic
---

# Algorithm 916: Faddeyeva & Voigt Evaluation

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 \) [1505.06848].

## 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:
```matlab
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
```
[1106.0151, 1505.06848]

## 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 [1505.06848].

## 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 [1106.0151, 1505.06848].
- **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 [1503.03806]. 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] (dbl)    |              | ~10–23                    |

[1505.06848]

## 7. Practical Usage and Recommendations

For MATLAB/Scilab, the interface allows flexible selection of precision:
```matlab
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)
```
[1106.0151]

For Fortran:
```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 [1505.06848].
- 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 [1106.0151, 1503.03806, 1505.06848].

Source: https://www.emergentmind.com/topics/algorithm-916