# Fast and Reliable Evaluation of the Distribution of Quadratic Forms of Gaussian Random Variables

Denis Arzelier, Florent Bréhard, Mioara Joldes

► To cite this version:

Denis Arzelier, Florent Bréhard, Mioara Joldes. Fast and Reliable Evaluation of the Distribution of Quadratic Forms of Gaussian Random Variables. 2026. hal-05507599

HAL Id: hal-05507599

https://hal.science/hal-05507599v1

Preprint submitted on 12 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Fast and Reliable Evaluation of the Distribution of
Quadratic Forms of Gaussian Random Variables

Denis Arzelier
LAAS-CNRS,

Université de Toulouse, CNRS,
07 Avenue Colonel Roche,
BP 54200,
Toulouse, 31031 France,
arzelier@laas.fr,

orcid: 0000-0003-0149-1369

Florent Bréhard

Univ. Lille, CNRS, Centrale Lille,
UMR 9189 CRISTAL
F-59000, Lille, France,
florent.brehard@univ-lille.fr,

orcid: 0000-0003-3909-2099

Mioara Joldes

LAAS-CNRS,

Université de Toulouse, CNRS,
07 Avenue Colonel Roche,
BP 54200,
Toulouse, 31031 France,
mmjoldes@laas.fr,

orcid: 0000-0001-7781-3745

Abstract

Quadratic forms in Gaussian random variables yield generalized noncentral chi-
square distributions which appear in many test statistics. Classical cdf methods
(e.g., characteristic-function inversion, saddlepoint approximations) can be effec-
tive but are not designed for reliable finite-precision computation and offer limited
structural or complexity insight. We study finite-precision cdf evaluation and its
binary complexity.

Building on earlier low-dimensional results, we note that relevant modified Laplace transforms of the cdf are available in closed-form and are D-finite functions. This reduces the general *d* dimensional cdf evaluation problem to evaluating a truncated holonomic power series with positive coefficients, together with explicit tail bounds. Next, we develop a reduced-rounding-error evaluation by switching from the standard alternating-sign scalar recurrence to an equivalent coupled recurrence in which all update coefficients are nonnegative. This “positivity feature” supports a sharper forward-error analysis and improves robustness of the finite-precision evaluation.

Finally, we propose an accelerated coefficient-generation algorithm. Using a
transposed multipoint-evaluation primitive, it computes collectively—via a trans-
posed Vandermonde/remainder-tree construction—the truncated contributions in-
duced by the rational factors of the holonomic equation. The remaining step, which assembles the solution by exponentiating a truncated series formed from these aggregated transforms, is carried out with FFT-based power-series exponentiation. Overall, the method achieves quasi-linear complexity in the maximum between the number of Gaussian random variables and the truncation order.

# 1 Introduction

Quadratic forms in Gaussian random variables arise throughout statistics and applied probability, e.g., as test statistics [18], in spectral analysis and optimal control [14], or space collision risk evaluation [24], [17]. In this work, we study the efficient and accurate evaluation of the cumulative distribution function (cdf)

$$F_Q(q) = \mathbb{P}[Q \le q], \quad q \ge 0,$$

of the Gaussian positive semidefinite quadratic form

$$Q = \Xi^\top A \Xi, \quad \Xi \sim \mathcal{N}(\mu, \Sigma), \quad A \succeq 0,$$

where $\Xi \in \mathbb{R}^d$ is a normal vector, with mean $\mu$ and covariance $\Sigma$.

A classical reduction expresses $Q$ as a sum of squares of $d$ independent normal variables $X_i \sim \mathcal{N}(\nu_i, \sigma_i^2)$. Writing the spectral decomposition $\Sigma^{1/2} A \Sigma^{1/2} = P \Lambda P^\top$, where $P$ is orthogonal and $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_d)$ is a diagonal matrix with entries $\lambda_i \ge 0$, define $b := P^\top \Sigma^{-1/2} \mu$ and let $Z \sim \mathcal{N}(0, I_d)$ (so that $Z_1, \dots, Z_d$ are independent standard normals), then

$$Q = \sum_{i=1}^{d} \lambda_i (Z_i + b_i)^2 = \sum_{i=1}^{d} X_i^2, \qquad (1)$$

with $\nu_i = \sqrt{\lambda_i} b_i$, and $\sigma_i^2 = \lambda_i$. Equivalently, $Q$ is a generalized noncentral chi-square variable, that is, a weighted sum of independent noncentral chi-square distributions [13, 14, 18].

## 1.1 Related works

Three main algorithmic paradigms enable fast evaluation of $F_Q(q)$.

A first major class of algorithms relies on a numerical inversion of the characteristic function of $Q$ using a particular instance of the inversion theorem [13], [9]. Although $F_Q$ has no closed-form in general, the characteristic function is explicit and can be inverted numerically [13]. Related refinements and implementations including efficient numerical integration strategies and accuracy guarantees, were developed later [9, 10, 8]. While these approaches are broadly applicable, their numerical performance depends strongly on the choice of integration contour and on the oscillatory nature of the integrand, which can be problematic in high accuracy or repeated-query settings.

A second line of work is based on saddlepoint approximations and related steepest-descent techniques [5]. These expansions originate in Daniels' work on approximating
the probability density function (pdf) or cdf of sample means of independent normal variables [6, 7]. Drawing on Watson lemma and a generalization of Laplace method [25], they typically provide accurate pdf/cdf approximations with few terms, and are especially effective for tail probabilities. Modern refinements include the Lugannani-Rice approximation, higher-order corrections [5, 25] and Barndorff-Nielsen-type formulae [3]. Saddlepoint methods have been adapted to the computation of the cdf of Q in [15]. Yet, this class of algorithms relies on numerically solving a saddlepoint equation and may require special care near singular points of the Laplace transform of the cdf.

A third family of methods is based on *convergent series expansions*; a unified presentation of several classes of quadratic forms is given in [14, 18]. Classical expansions include power series and mixture representations in terms of central or noncentral chi-square distributions. In particular, Ruben established a mixture expansion of the distribution function as a convex combination of chi-square CDFs and derived recurrence relations for the coefficients [22, 23]. Such convergent expansions can be highly effective when their coefficients can be generated in a numerically stable way, but direct power series are known to suffer from catastrophic cancellation when coefficients alternate in sign [11].

More recently, holonomic methods [21] have been proposed to compute a holonomic system of linear (ordinary or partial) differential equations for the cdf. Likewise, the *D-finiteness* of Laplace transforms of cdfs—together with their closed-form availability—yields a more refined strategy for efficient finite-precision evaluation, developed in dimensions 2 and 3 for orbital collision probability and related risk-assessment problems [24]. The approach computes preconditioned series with *nonnegative coefficients*, avoiding catastrophic cancellation and giving rigorous truncation error bounds, and was complemented in [2] by a rounding-error analysis based on majorizing-series techniques.

We build on this symbolic-numeric viewpoint recalled below and sharpen the complexity/accuracy analysis in higher dimension.

## 1.2 Problem setup

Related to equation (1), we introduce the function

$$g(\xi) := \mathbb{P}\left[\sum_{i=1}^{d} X_i^2 \le \xi q\right], \quad \xi \ge 0, \quad q \ge 0,$$

which can be written as the *d*-dimensional Gaussian integral

$$g(\xi) = \frac{1}{\prod_{i=1}^{d} \sqrt{2\pi}\sigma_i} \int \exp\left(-\sum_{i=1}^{d} \frac{(x_i - \nu_i)^2}{2\sigma_i^2}\right) dx_1 \cdots dx_d. \quad (2)$$

For a locally integrable function $h: \mathbb{R}_{\ge 0} \to \mathbb{R}$, we denote by

$$\mathcal{L}_h(\lambda) := \int_0^{+\infty} e^{-\lambda\xi} h(\xi) \, d\xi, \quad \text{Re}(\lambda) > 0, \quad (3)$$
its unilateral Laplace transform. Following [24], we define the *modified Laplace transform* of the function $h$ as:

$$ \hat{h}(\lambda) := \frac{\mathcal{L}_h(\lambda^{-1})}{\lambda^{d/2+1}}. \qquad (4) $$

In particular, the closed-form of the modified Laplace transform $\hat{g}$ of $g$ defined in (2), may be readily obtained as

$$ \hat{g}(\lambda) = \frac{\mathcal{L}_g(\lambda^{-1})}{\lambda^{d/2+1}} = \frac{\left(\prod_{i=1}^{d} \sqrt{p_i}\right) \exp\left(-\sum_{i=1}^{d} \frac{m_i}{1+p_i\lambda}\right)}{\prod_{i=1}^{d} \sqrt{1+p_i\lambda}}, \qquad (5) $$

where the parameters $p_i$ and $m_i$ are defined w.l.o.g. by

$$ p_i := \frac{q}{2\sigma_i^2}, \quad m_i := \frac{\nu_i^2}{2\sigma_i^2} \quad (1 \le i \le d), \quad p_1 = \max_{1 \le i \le d} p_i > 0. \qquad (6) $$

Hence, $\hat{g}$ is analytic for $|\lambda| < 1/p_1$ and it is straightforward to verify that $\hat{g}$ satisfies a first-order linear differential equation with polynomial coefficients. Given the parameters in (6), we consider the finite-precision evaluation of $g(1) = F_Q(q)$ and analyze its binary complexity.

**Notations** For sets $S$, $\#S$ denotes cardinality; $\mathbb{R}_{\ge 0}$ denotes the non-negative real numbers and similarly for $\mathbb{N}_{>0}$, $\mathbb{R}_{>0}^d$. For real polynomials, $\deg(\cdot)$ is the degree and $\mathbb{R}[z]_{<N}$ denotes the space of polynomials in $z$ with real coefficients of degree $< N$. For a polynomial or formal power series $P(z) = \sum_{n \ge 0} p_n z^n$, we write $P(z) \% z^N$ for truncation modulo $z^N$. For formal power series $A(\lambda) = \sum_{n \ge 0} a_n \lambda^n$ and $B(\lambda) = \sum_{n \ge 0} b_n \lambda^n$ with $b_n \ge 0$, we write $A \ll B$ if $|a_n| \le b_n$ for all $n$, and $A \gg 0$ if $a_n \ge 0$ for all $n$. The Euler Gamma function is denoted by $\Gamma$.

## 1.3 Contributions

A direct $d$-dimensional generalization of the algorithm of [24] follows from (5). Section 2 recalls the construction, including the preconditioning that yields a series with nonnegative coefficients and the resulting truncation bounds, as well as the rounding-error analysis [2]. Despite positivity of the target coefficients, the generating recurrence has mixed-sign coefficients, so finite-precision cancellations may occur and worsen as $d$ increases.

To mitigate these cancellations, Section 3 introduces an alternative recurrence that computes the same coefficients using only additions, leading to improved numerical robustness and the complexity/precision trade-off summarized by the next theorem. We also derive effective truncation and evaluation error bounds, together with an available implementation [1].

**Theorem 1.1.** Given parameters $(p_i, m_i)_{i=1}^d$, and letting $M = \sum_{i=1}^d m_i$, $g(1)$ can be approximated with relative accuracy $2^{-\tau}$ by summing $N = O(p_1+M+\sum_i \ln(p_i/p_i)+\tau)$ terms of the power series of $g$ in $O(Nd)$ arithmetic operations with a working precision of $t = \tau + \ln(N) + \ln(d) + \ln(M) + O(1)$ bits, yielding a total bit complexity of $\tilde{O}((p_1+M+\sum_i \ln(p_i/p_i)+\tau)d\tau)$ (where $\tilde{O}()$ hides polylogarithmic factors in $p_1, d, M$ and $\tau$).
Section 4 presents a second, currently mostly theoretical (but with an on-going implementation), approach that achieves a quasi-linear dependence on the target precision, at the cost of an absolute (rather than relative) error guarantee, as stated in the following theorem.

**Theorem 1.2.** Given parameters $(p_i, m_i)_{i=1}^d$, $g(1)$ can be approximated with absolute accuracy $2^{-\tau}$ by summing $N = O(p_1 + \tau)$ terms of the power series of $g$ in $\tilde{O}(N+d)$ arithmetic operations with a working precision of $t = \tilde{O}(\tau + \ln(N) + \ln(d) + \ln(M))$ bits, yielding a total bit complexity of $\tilde{O}((p_1+d+M+\tau)(\tau+\ln(M)))$ (where $\tilde{O}()$ hides polylogarithmic factors in $p_1, d$ and $\tau$).

We conclude in Section 5 with numerical experiments, including a comparison with an existing Matlab implementation.

## 2 D-finite series with nonnegative coefficients

The power series of $g$ often exhibits mixed-sign coefficients, causing cancellation in the finite-precision series evaluation. Following [11, 24], we therefore precondition $g$ by $f(\xi) = e^{p\xi}g(\xi)$, with a carefully chosen parameter $p \in \mathbb{R}$, so that the power series coefficients of $f$ are nonnegative and hence, evaluate $g(1)$ as follows.

**Proposition 2.1.** The target probability is given by

$$ \mathbb{P}[Q \le q] = g(1) = e^{-p_1} \sum_{n \ge 0} c_n, \qquad (7) $$

with $c_n \ge 0$ satisfying an order-2d linear recurrence given in (11).

**Proof sketch adapted from [11, 24]** Firstly, eq. (4) and the shift relation $\mathcal{L}_f(\lambda) = \mathcal{L}_g(\lambda - p)$, for $\mathrm{Re}(\lambda) > p$ yield

$$ \hat{f}(\lambda) = \frac{\mathcal{L}_f(\lambda^{-1})}{\lambda^{d/2+1}} = (1-p\lambda)^{-d/2-1} \hat{g}\left(\frac{\lambda}{1-p\lambda}\right), \quad |\lambda| < 1/p. $$

Let $\alpha_i = p - p_i$ for $1 \le i \le d$, and notice that substituting $\lambda$ for $\lambda/(1-p\lambda)$ in the term $1+p_i\lambda$ (cf. eq. (5)) gives $(1-\alpha_i\lambda)/(1-p\lambda)$. Furthermore, since the power series of $(1-a\lambda)^{-1/2}$ and of $\exp(-m(1-p\lambda)/(1-a\lambda))$ have nonnegative coefficients for, respectively, $a \ge 0$ and $m \ge 0$ and $0 \le a \le p$, this motivates choosing $p \ge \max_i p_i$. From (6), take $p = p_1$, so $\alpha_1 = 0$, and finally obtain:

$$ \hat{f}(\lambda) = \frac{\left(\prod_{i=1}^{d} \sqrt{p_i}\right) \exp \left[-\sum_{i=1}^{d} m_i \frac{1-p_i\lambda}{1-\alpha_i\lambda}\right]}{(1-p_1\lambda) \prod_{i=2}^{d} \sqrt{1-\alpha_i\lambda}}. \qquad (8) $$
Secondly, differentiating (8) gives

$$ \hat{f}'(\lambda) = \varphi(\lambda) \hat{f}(\lambda), \quad \hat{f}(0) = \left(\prod_{i=1}^{d} \sqrt{p_i}\right) \exp\left(-\sum_{i=1}^{d} m_i\right), \qquad (9) $$

where $\varphi(\lambda)$ is the rational function:

$$ \begin{aligned} \varphi(\lambda) &= \frac{p_1}{1 - p_1 \lambda} + \sum_{i=1}^{d} \left( \frac{\alpha_i / 2}{1 - \alpha_i \lambda} + \frac{p_i m_i}{(1 - \alpha_i \lambda)^2} \right), \\ &= m_1 p_1 + \frac{p_1}{1 - p_1 \lambda} + \sum_{i=2}^{d} \left( \frac{\alpha_i / 2}{1 - \alpha_i \lambda} + \frac{p_i m_i}{(1 - \alpha_i \lambda)^2} \right). \end{aligned} \qquad (10) $$

Note that $\varphi(\lambda) = \frac{P(\lambda)}{Q(\lambda)}$, with $P, Q \in \mathbb{R}_{\ge 0}[\lambda]$ defined by $P(\lambda) = \sum_{j=0}^{2d-1} (-1)^j P_j \lambda^j$ and $Q(\lambda) = \sum_{j=0}^{2d-1} (-1)^j Q_j \lambda^j$, with $Q_0 = 1$.

Letting $\hat{f}(\lambda) = \sum_{n \ge 0} c_n \Gamma(n + \frac{d}{2} + 1) \lambda^n$, the differential equation $Q(\lambda) \lambda \hat{f}'(\lambda) - P(\lambda) \lambda \hat{f}(\lambda) = 0$ translates into a recurrence for $c_n$,

$$ \begin{aligned} c_0 &= \frac{\hat{f}(0)}{\Gamma(\frac{d}{2} + 1)}, \\ nc_n + \sum_{i=1}^{\min(2d,n)} (-1)^i \frac{\Gamma(n - i + \frac{d}{2} + 1)}{\Gamma(n + \frac{d}{2} + 1)} [P_{i-1} + Q_i(n-i)] c_{n-i} = 0, \quad n \ge 1. \end{aligned} \qquad (11) $$

Finally, a term-by-term Laplace inversion argument (see [26, Chap. 9]), yields the convergent expansion $f(\xi) = \sum_{n \ge 0} c_n \xi^n$ and we conclude with $g(1) = e^{-p_1} f(1)$. □ Hence, from (7), $g(1)$ can be approximated by truncating the series at order $N$. The truncation error $r_N := g(1) - e^{-p_1} \sum_{n=0}^{N-1} c_n = e^{-p_1} \sum_{n=N}^{\infty} c_n$ satisfies the following:

**Proposition 2.2 (Truncation error).** Let $\delta \in (0, 1)$.

(i) (absolute truncation error) $r_N \le \delta$ for $N \ge N_1$, where

$$ N_1 = \left\lceil \max \left( 2ep_1, \frac{\ln(\delta^{-1})}{\ln 2} \right) - \frac{d}{2} \right\rceil. $$

(ii) (relative truncation error) $r_N \le \delta g(1)$ for $N \ge N_2$, where

$$ N_2 = \left\lceil \max \left( 2ep_1, \frac{\ln(\delta^{-1})}{\ln 2} + \frac{M}{\ln 2} + \frac{1}{2\ln 2} \sum_{i=1}^{d} \ln \frac{p_i}{p_i} \right) - \frac{d}{2} \right\rceil. $$

For $\delta = 2^{-\tau}$, $N_1 = O(p_1 + \tau)$ and $N_2 = O(p_1 + \tau + M + \sum_i \ln \frac{p_i}{p_i})$.
Proof. Letting $\hat{f}(\lambda) = \sum_{n \ge 0} \hat{c}_n \lambda^n$, we first show that

$$ \hat{c}_0 p_1^n \le \hat{c}_n \le p_1^{n+\frac{d}{2}}, \quad n \ge 0. \qquad (12) $$

Let $\hat{f}^*(\lambda) := (1 - p_1\lambda)\hat{f}(\lambda) = \sum_{n \ge 0} \hat{c}_n^* \lambda^n \gg 0$. Since $\hat{f}^*(\lambda) \gg \hat{c}_0$, we get $\hat{f}(\lambda) \gg \hat{c}_0/(1 - p_1\lambda)$, hence $\hat{c}_n \ge \hat{c}_0 p_1^n$. Moreover, $\hat{f}^*(1/p_1) = \sum_{n \ge 0} \hat{c}_n^* p_1^{-n} = p_1^{d/2}$ (using $p_i/(1 - \alpha_i/p_1) = p_1$), so

$$ \hat{c}_n = \sum_{k=0}^{n} p_1^{n-k} \hat{c}_k^* \le p_1^n \sum_{k \ge 0} \hat{c}_k^* p_1^{-k} = p_1^{n+\frac{d}{2}}. $$

(i) Using $\Gamma(x+1)\Gamma(y+1) \le \Gamma(x+y+1)$ and (12),

$$ \sum_{n \ge N} c_n \le \sum_{n \ge N} \frac{p_1^{n+\frac{d}{2}}}{\Gamma(n+\frac{d}{2}+1)} \le e^{p_1} \frac{p_1^{N+\frac{d}{2}}}{\Gamma(N+\frac{d}{2}+1)}. $$

By Stirling inequality, $\Gamma(x+1) \ge (x/e)^x$ for $x > 1$, hence

$$ r_N = e^{-p_1} \sum_{n \ge N} c_n \le \frac{p_1^{N+\frac{d}{2}}}{\Gamma(N+\frac{d}{2}+1)} \le \left( \frac{e p_1}{N+\frac{d}{2}} \right)^{N+\frac{d}{2}}. $$

Thus, $N \ge N_1$ implies $ep_1/(N + \frac{d}{2}) \le 1/2$ and $r_N \le \delta$.

(ii) Let $\alpha = \hat{c}_0 p_1^{-\frac{d}{2}} = \exp(-M - \frac{1}{2} \sum_{i=1}^{d} \ln \frac{p_1}{p_i})$ be the ratio between the lower and upper bounds in (12). Also, $e^{-p_1} \sum_{n=0}^{\infty} \frac{p_1}{\Gamma(n+\frac{d}{2}+1)} \le 1$ holds since $\sum_{n=0}^{\infty} \frac{x^n}{\Gamma(n+1)}$ and $\sum_{n=0}^{\infty} \frac{x^{n+\frac{1}{2}}}{\Gamma(n+\frac{3}{2})}$ are the respective power series expansions of $e^x$ and $e^x \operatorname{erf}(\sqrt{x}) \le e^x$. Therefore, we have $\alpha \le g(1) \le 1$, and (ii) follows from (i) by replacing $\delta$ with $\alpha\delta$. $\square$

In Appendix B, several refined truncation bounds are derived and were also implemented in practice.

Algorithm CDFREC evaluates the truncated series by computing the coefficients $c_n$ via the recurrence (11). In finite precision arithmetic, the computed coefficients $\tilde{c}_n$ differ from $c_n$ because of rounding in the precomputed data (the $\alpha_i$ and the derived $P_j, Q_j$) and rounding in each recurrence step. Following the majorant-series approach of [19, 2], one can bound the sum of the errors $\delta_n = c_n - \tilde{c}_n$ by noting that $\delta(\lambda) = \sum_{n \ge 0} \delta_n \Gamma(n + \frac{d}{2} + 1) \lambda^n$ satisfies an inhomogeneous first-order linear differential equation, with homogeneous term given by (9) and inhomogeneous part driven by the series of local errors injected at iteration $n$ of the loop evaluating $\tilde{c}_n$. Coefficientwise positivity of the majorants then yields an explicit uniform control of the accumulated rounding error on $f(1)$ (and hence $g(1)$).
**Algorithm 1:** CDFREC(d, $(p_i)_{i=1}^d$, $(m_i)_{i=1}^d$, N)

**input:** $d \in \mathbb{N}_{>0}$; $(p_i)_{i=1}^d$, $(m_i)_{i=1}^d \in (\mathbb{R}_{>0})^d$; $N \in \mathbb{N}_{>0}$

**output:** an approximation $s$ of $g(1)$ cf. eq. (7).

1 **for** $i = 1$ **to** $d$ **do**

2 $\varphi(\lambda) \leftarrow p_1 m_1 + \frac{p_1}{1 - p_1 \lambda} + \sum_{i=2}^{d} \left( \frac{\alpha_i / 2}{1 - \alpha_i \lambda} + \frac{p_i m_i}{(1 - \alpha_i \lambda)^2} \right)$

3 $P(\lambda) = P_0 - P_1 \lambda + \dots - P_{2d-1} \lambda^{2d-1} \leftarrow \text{numerator}(\varphi)$

4 $Q(\lambda) = 1 - Q_1 \lambda + \dots - Q_{2d-1} \lambda^{2d-1} \leftarrow \text{denominator}(\varphi)$

5 $c_0 \leftarrow \frac{\sqrt{\prod_i p_i} \exp(-\sum_i m_i)}{\Gamma(d/2 + 1)}$; $s \leftarrow c_0$

6 **for** $n = 1$ **to** $N$ **do**

7 $c_n \leftarrow 0$; $\text{fac} \leftarrow 1$

8 **for** $i = 1$ **to** $\min(2d, n)$ **do**

9 $\text{fac} \leftarrow (d/2 + n + 1 - i)$ $\text{fac}$

10 $c_n \leftarrow c_n + (-1)^{i-1} \frac{P_{i-1} + Q_i(n-i)}{\text{fac}} c_{n-i}$

11 **end**

12 $c_n \leftarrow c_n/n$; $s \leftarrow s + c_n$

13 **end**

14 **return** $s \leftarrow \exp(-p_1) s$

**Proposition 2.3 (Evaluation error for Algorithm CDFREC).** Consider radix-2, precision-t, round-to-nearest floating-point arithmetic with unbounded exponent range, unit round-off $u = 2^{-t}$ and define $\gamma_k = \frac{ku}{1-ku}$, for $ku \le 1$. Assume that the coefficients of P and Q are pre-computed with relative accuracy bounded by $\gamma_{K_{\text{coef}}}$.
Then, the total rounding error accumulated while computing $f(1)$ satisfies

$$ \sum_{n=0}^{+\infty} |\tilde{c}_n - c_n| = |\delta|(1) = (e_0 + \gamma C^{\#} + o(u)) f(1), $$

with $e_0 = (1 + \gamma_{d+7})e^{\gamma_{d-1} \sum_i^d m_i} - 1 \ge |\delta_0|/c_0$, $\gamma := \gamma_{14d+1+K_{\text{coef}}}$ accounting to local rounding errors and

$$ C^{\#} := \int_{0}^{1} W^{\#}(\tau) d\tau, $$

with $W^{\#}$ an explicit degree-2d-1 polynomial with nonnegative coefficients $w_r^{\#} = O(\frac{p_r}{r!})$.

The proof of Prop. 2.3 is a technical but straightforward generalization of [2]. It is given for completeness in Appendix C.
## 3 A new recurrence with reduced evaluation error

While $(c_n)_{n \ge 0}$ are nonnegative and $f(1) = \sum_{n \ge 0} c_n$ is a sum of positive terms, the direct recurrence (11) derived from (9) via $\varphi = P/Q$ contains alternating signs. Therefore, the finite precision computation of $c_n$ can still be affected by the cancellation phenomenon. We thus derive an equivalent recurrence whose update coefficients are all nonnegative.

### 3.1 A recurrence with positive coefficients

We introduce the following power series expansions:

$$ 
\begin{aligned}
\hat{f}(\lambda) &= \sum_{n=0}^{\infty} \hat{c}_n \lambda^n, & \frac{\hat{f}(\lambda)}{1 - \alpha_i \lambda} &= \sum_{n=0}^{\infty} \hat{a}_{i,n} \lambda^n, \quad (2 \le i \le d), \\
\frac{\hat{f}(\lambda)}{1 - p_1 \lambda} &= \sum_{n=0}^{\infty} \hat{a}_{1,n} \lambda^n, & \frac{\hat{f}(\lambda)}{(1 - \alpha_i \lambda)^2} &= \sum_{n=0}^{\infty} \hat{b}_{i,n} \lambda^n, \quad (2 \le i \le d).
\end{aligned}
$$

Using the differential equation (9) with (10) and identifying the coefficients yields, for $n \ge 1$,

$$ 
\left\{
\begin{aligned}
\hat{c}_n &= \frac{1}{n} \left( p_1 m_1 \hat{c}_{n-1} + p_1 \hat{a}_{1,n-1} + \sum_{i=2}^{d} \left( \frac{\alpha_i}{2} \hat{a}_{i,n-1} + p_i m_i \hat{b}_{i,n-1} \right) \right), \\
\hat{a}_{1,n} &= p_1 \hat{a}_{1,n-1} + \hat{c}_n, \\
\hat{a}_{i,n} &= \alpha_i \hat{a}_{i,n-1} + \hat{c}_n, \quad (2 \le i \le d), \\
\hat{b}_{i,n} &= \alpha_i \hat{b}_{i,n-1} + \hat{a}_{i,n}, \quad (2 \le i \le d).
\end{aligned}
\right.
$$

with initial conditions $\hat{a}_{1,0} = \hat{a}_{i,0} = \hat{b}_{i,0} = \hat{c}_0 = \hat{f}_0$ for $2 \le i \le d$.

To perform the inverse Laplace transform as in the previous section, we rescale the coefficients by setting $\hat{c}_n := \Gamma(n + \frac{d}{2} + 1) c_n$, $\hat{a}_{i,n} := \Gamma(n + \frac{d}{2} + 1) a_{i,n}$, and $\hat{b}_{i,n} := \Gamma(n + \frac{d}{2} + 1) b_{i,n}$, which gives the updated relations:

$$ 
\left\{
\begin{aligned}
c_n &= \frac{p_1 m_1 c_{n-1} + p_1 a_{1,n-1} + \sum_{i=2}^{d} \left( \frac{\alpha_i}{2} a_{i,n-1} + p_i m_i b_{i,n-1} \right)}{n(n+d/2)}, \\
a_{1,n} &= \frac{p_1}{n+d/2} a_{1,n-1} + c_n, \\
a_{i,n} &= \frac{\alpha_i}{n+d/2} a_{i,n-1} + c_n, \quad (2 \le i \le d), \\
b_{i,n} &= \frac{\alpha_i}{n+d/2} b_{i,n-1} + a_{i,n}, \quad (2 \le i \le d).
\end{aligned}
\right.
$$

All coefficients appearing in these recurrences are nonnegative. Consequently, when their initial conditions are nonnegative, the sequences $(c_n)$, $(a_{i,n})$, and $(b_{i,n})$ remain nonnegative, which significantly reduces cancellation effects compared with the scalar recurrence (11). The implementation of this new recurrence is given in Algorithm CDFPOSREC.
**Algorithm 2:** CDFPOSREC(d, (pᵢ)ᵢ=1, (mᵢ)ᵢ=1, N)

**input:** d ∈ ℝ>0; (pᵢ)ᵢ=1, (mᵢ)ᵢ=1 ∈ (ℝ>0)ᵈ; N ∈ ℝ>0

**output:** an approximation s of g(1) cf. eq. (7).

1 $\alpha_i \leftarrow p_1 - p_i \quad (1 \le i \le d)$; $c_0 \leftarrow \frac{\sqrt{\prod_i p_i \exp(-\sum_i m_i)}}{\Gamma(d/2 + 1)}$

2 $a_{i,0} \leftarrow c_0 \quad (1 \le i \le d)$; $b_{i,0} \leftarrow c_0 \quad (2 \le i \le d)$; $s \leftarrow c_0$

3 **for** n = 1 **to** N **do**

4 $c_n \leftarrow p_1 m_1 c_{n-1} + p_1 a_{1,n-1}$

5 **for** i = 2 **to** d **do** $c_n \leftarrow c_n + \frac{\alpha_i}{2} a_{i,n-1} + p_i m_i b_{i,n-1}$

6 $c_n \leftarrow \frac{c_n}{n(n+d/2)}$; $a_{1,n} \leftarrow \frac{p_1}{n+d/2} a_{1,n-1} + c_n$

7 **for** i = 2 **to** d **do**

8 $a_{i,n} \leftarrow \frac{\alpha_i}{n+d/2} a_{i,n-1} + c_n$

9 $b_{i,n} \leftarrow \frac{\alpha_i}{n+d/2} b_{i,n-1} + a_{i,n}$

10 **end**

11 $s \leftarrow s + c_n$

12 **end**

13 **return** $s \leftarrow \exp(-p_1) s$

## 3.2 Roundoff error analysis

**Floating-point model** We assume radix-2, precision-t, round-to-nearest floating-point (FP) arithmetic with unbounded exponent range and unit roundoff $u = 2^{-t}$. We use standard error models [12, Chap. 2.2.]: for any real $x$, the round to nearest FP number is $\text{RN}(x) = x(1 + \delta) = \frac{x}{1+\delta'}$, with $|\delta|, |\delta'| \le u$. For each basic operation $\circ \in \{+, -, \times, /\}$, the computed result is $\text{RN}(x \circ y)$. Multiplications by powers of two are assumed exact (hence $\text{RN}(x/2) = x/2$). Composite expressions are evaluated from left to right. We assume that $\exp$ and $\sqrt{}$ as well as $\Gamma(d/2 + 1)$ are faithfully rounded, i.e., their relative errors are bounded by $2u$.

**Positive-sum perturbation and θ-calculus.** Variables with a tilde denote computed quantities. We also use Higham [12]'s $\theta_k$-notation: each occurrence of $\theta_k$ denotes a (possibly different) quantity of the form $1+\theta_k = \prod_{j=1}^k (1+\delta_j)^{\pm 1}$ with $|\delta_j| \le u$. Assuming $ku < 1$, one has $|\theta_k| \le \gamma_k$ where $\gamma_k := \frac{ku}{1-ku}$. Moreover, multiplicative perturbations can be chained using the standard $\theta$-notation: for any $r, s$, $(1+\theta_r)(1+\theta_s) = 1+\theta_{r+s}$. Since all intermediate quantities in Algorithm CDFPOSREC are nonnegative, relative perturbations can be merged in sums [12, Chap. 3]: if $\tilde{s} = \sum_{j=1}^m s_j(1+\eta_j)$ with $s_j \ge 0$ and $|\eta_j| \le \gamma_k$, then $\tilde{s} = s(1+\eta)$ with $s = \sum_{j=1}^m s_j$ and $|\eta| \le \gamma_k$.
**Initialization.** The computed $\tilde{\alpha}_i = \text{RN}(p_1 - p_i)$ satisfy $\tilde{\alpha}_i = \alpha_i(1 + \delta_1)$, with $|\delta_1| \le u$. Also $\tilde{c}_0$ computed from the expression of $c_0$ by left-to-right evaluation is given by:

$$ \tilde{c}_0 = c_0 \psi (1 + \beta), \quad \psi := \exp\left(\eta \sum_{i=1}^{d} m_i\right), \qquad (13) $$

with $|\eta| \le \gamma_{d-1}$ and $|\beta| \le \gamma_{d+7}$, assuming $(d+7)u < 1$. This follows from setting $P := \prod_{i=1}^d p_i$, $M := \sum_{i=1}^d m_i$ and noting that with positive operands and left-to-right evaluation,

$$ \tilde{P} = P(1 + \theta_{d-1}), \quad |\tilde{M} - M| \le \gamma_{d-1}M. $$

Thus $\exp(-\tilde{M}) = \exp(-M)\psi$ with $|\eta| \le \gamma_{d-1}$. Assuming faithful rounding for $\sqrt{}$, $\exp$ and $\Gamma(d/2 + 1)$ and evaluating $\tilde{c}_0$ left-to-right yields $\tilde{c}_0 = c_0 \psi \sqrt{1 + \theta_{d-1}} (1 + \theta_2)(1 + \delta_\sqrt{})(1 + \delta_e)(1 + \delta_\Gamma)$, with $|\delta_\sqrt{}|, |\delta_e|, |\delta_\Gamma| \le 2u$ and $\theta_2$ accounting for the final multiplication and division.

**Propagation** Let $\tilde{c}_n, \tilde{a}_{i,n}, \tilde{b}_{i,n}$ be the numerical values produced by CDFPOSREC. Using the roundoff error algebra, we obtain the following relative rounding error representation.

**Lemma 3.1** (Relative errors on $c_n$). For all $n \ge 1$ and under the assumption that $u((2d+4)n + d + 9) < 1$, we have

$$ \tilde{c}_n = \psi(1 + \varepsilon_n^{(c)})c_n, \quad |\varepsilon_n^{(c)}| \le \gamma_{(2d+4)n+d+7}, $$

$$ \tilde{a}_{i,n} = \psi(1 + \varepsilon_{i,n}^{(a)})a_{i,n}, \quad |\varepsilon_{i,n}^{(a)}| \le \gamma_{(2d+4)n+d+8}, $$

$$ \tilde{b}_{i,n} = \psi(1 + \varepsilon_{i,n}^{(b)})b_{i,n}, \quad |\varepsilon_{i,n}^{(b)}| \le \gamma_{(2d+4)n+d+9}. $$

*Proof.* We proceed by induction and consider first the left-to-right evaluation of $c_1$. The numerator is a sum of $2d$ nonnegative terms. Forming any term uses at most two multiplications, the accumulation uses $2d-1$ additions, and the scaling by $1/(n(n+d/2))$ uses three basic operations. Hence the update introduces at most $k \le 2+(2d-1)+3 = 2d+4$ new rounding factors. The first summand is thus $p_1 m_1 c_0 \psi (1+\beta)(1+\theta_k)$, and $\psi > 0$ and $|\beta| \le \gamma_{d+7}$ given in equation (13). The other terms in the summation all contain $c_0 \psi (1+\beta)(1+\theta_{k'})$ with $k' \le 2d+4$. Since all the terms are positive, relative errors can be merged in sums:

$$ |\tilde{c}_1 - \psi c_1| \le \psi \gamma_{2d+4+d+7} c_1, $$

which is equivalent to

$$ \tilde{c}_1 = \psi(1 + \varepsilon_1^{(c)})c_1, \quad |\varepsilon_1^{(c)}| \le \gamma_{2d+4+d+7}. $$

Similarly, for $\tilde{a}_{i,1}$ the most perturbed summand is $\tilde{c}_1$ affected by an extra addition and for $\tilde{b}_{i,1}$ the most perturbed summand is $\tilde{a}_{i,1}$ affected by an extra addition, which gives
the required inequalities. Now suppose that step $n-1$ holds. At step $n$, the first summand of $c_n$ is affected (as explained before) by $k \le 2d+4$ operations, so we get the term $p_1m_1c_{n-1}\psi(1+\varepsilon_{n-1}^{(c)})(1+\theta_k)$, while the second one is $p_1a_{1,n-1}\psi(1+\varepsilon_{1,n-1}^{(a)})(1+\theta_{k-1})$ and the remaining terms accounting for smaller perturbations. This implies that $|\varepsilon_n^{(c)}| \le \gamma_{n(2d+4)+d+7}$. A similar analysis provides the remaining bounds. $\square$

Finally, we obtain the following relative error.

**Lemma 3.2.** Denote $C_N := \sum_{n=0}^{N} c_n$, $g_N := e^{-p_1} C_N$ and $\tilde{C}_N$, $\tilde{g}_N$ the FP output (with $\tilde{C}_N$ obtained from left-to-right evaluation, and exp faithfully rounded). Under the assumption in Lemma 3.1 and ((2d+4)N+d+11)u < 1, we have

$$ \tilde{g}_N = \psi(1 + \varepsilon_N^{(g)}) g_N, \quad |\varepsilon_N^{(g)}| \le \gamma_{(2d+4)N+d+11}, $$

$$ \text{where } \psi = \exp(\eta \sum_{i=1}^{d} m_i) \text{ with } |\eta| \le \gamma_{d-1}. $$

*Proof.* In the running sum, the summands

$$ \tilde{c}_n = \psi c_n (1 + \theta_{(2d+4)n+d+7}), $$

are further affected by $N, N, N-1, \dots, 1$ additions, and since all terms are nonnegative, the sum-of-positives rule yields

$$ \tilde{C}_N = \psi C_N (1 + \varepsilon_N^C), \quad |\varepsilon_N^C| \le \gamma_{(2d+4)N+d+8}. $$

Finally, $\exp(-p_1)$ is faithfully rounded (relative error $\le 2u$, absorbed as a $\theta_2$-term) and the final multiplication contributes one more rounding: $\tilde{g}_N = e^{-p_1} \psi C_N (1+\varepsilon_N^C)(1+\theta_3)$, which implies $|\varepsilon_N^{(g)}| \le \gamma_{(2d+4)N+d+11}$. $\square$

*Proof of Theorem 1.1.* According to Lemma 3.2, a working precision of $t = \tau + \ln(\psi - 1) + \ln(\varepsilon_N^{(g)}) + O(1) = \tau + \ln(N) + \ln(d) + \ln(M) + O(1)$ is sufficient to guarantee a relative error bounded by $2^{-\tau}$, with $N = O(p_1+M+\sum_i \ln(p_i/p_i)+\tau)$ by Prop. 2.2(ii), hence the claimed overall asymptotic bit complexity. $\square$

# 4 A faster algorithm for high dimension

We now present Algorithm CDFFAST, which runs in arithmetic complexity $\tilde{O}(N+d)$ instead of $O(Nd)$ for CDFREC and CDFPosREC. Our analysis is currently limited to an absolute error setting, making this algorithm impractical for accurately estimating very low probabilities with a reasonable<fnref n="1" /> working precision $t$. By “reasonable”, we mean $t/N = o(1)$ but still $(M+d)^\alpha N^\beta u = o(1)$ where $u = 2^{-t}$ and for relatively small constants $\alpha, \beta$. The main steps are:

<footnote n="1">¹The problem becomes dramatically simpler if we ask for an accuracy $\tau$ (and hence a working precision $t$) at least of the order of $N+d$.</footnote>
1. Approximate the power series of $\ln(\hat{f}^*)$ up to order $N-1$, where $\hat{f}^*(\lambda) = (1 - p_1\lambda)\hat{f}(\lambda)$;

2. Approximate the power series of $\hat{f}^*(\lambda)$ by using a numerical Fast Fourier Transform approach for the exponentiation on the circle of radius $\frac{\nu}{p}$ in $\mathbb{C}$ for a well-chosen $\nu \in (0, 1)$;

3. Recover the coefficients of $\hat{f}$ from those of $\hat{f}^*$ and apply the inverse modified Laplace transform.

Section 4.1 focuses on a new algorithm VANDERT related to transposed multipoint evaluation to perform step (1) numerically with arithmetic complexity $\tilde{O}(N+d)$. Steps (2) and (3) of CDFFAST are presented in Section 4.2 together with the final error bound and complexity estimates.

**Fast Fourier Transform** Complexity results in this section heavily rely on the Fast Fourier Transform (FFT). Let $DFT(P, N)$ denote the Discrete Fourier Transform that maps a polynomial $P(z)$ of degree less than $N$ to the vector $(P(\omega^j))_{j=0}^{N-1}$, where $\omega = \exp(\frac{2i\pi}{N})$, and iDFT($v, N$) the inverse operation computing the unique polynomial of degree less than $N$ that interpolates the $(v_j)_{j=0}^{N-1}$ at $\omega^j$. These maps can be computed using Fast Fourier Transform (FFT) algorithms (hereafter denoted by FFT and iFFT) using $O(N \ln(N))$ FP operations. This notably makes it possible to multiply two degree $N$ polynomials in $O(N \ln(N))$ FP operations.

Let $\|P\|_1 = \sum_n |a_n|$ for $P(z) = \sum_n a_n z^n$ a polynomial or convergent power series on the unit disk, and $\|v\|_\infty = \max_j |v_j|$ for $v = (v_j)_j$.

**Lemma 4.1** (see, e.g., [12, §24.1]).

(i) $\| \text{DFT}(P, N) \|_\infty \le \|P\|_1$ and $\| \text{iDFT}(v, N) \|_1 \le N \|v\|_\infty$ for all $P(z)$ and $(v_j)_j$.

(ii) $\| \text{FFT}(P, N) - \text{DFT}(P, N) \|_\infty = O(\|P\|_1 \ln(N)u)$ for all $P(z)$.

(iii) $\| \text{iFFT}(v, N) - \text{iDFT}(v, N) \|_1 = O(\|v\|_\infty N \ln(N)u)$ for all $v$.

(iv) For $P, Q$ of degree $N$ and $R \approx PQ$ computed by FFT multiplication, $\|R-PQ\|_1 = O(\|P\|_1 \|Q\|_1 N \ln(N)u)$.

## 4.1 Transposed multipoint evaluation problem

Given $\alpha = (\alpha_i)_{i=1}^d$ $d$ points in $[0, 1]$ with weights $w = (w_i)_{i=1}^d$, the transposed multipoint evaluation problem consists in computing the first $N$ coefficients of the power series:

$$ \sum_{i=1}^{d} \frac{w_i}{1 - \alpha_i z} \% z^N = \sum_{n=0}^{N-1} \left( \sum_{i=1}^{d} w_i \alpha_i^n \right) z^n. \quad (14) $$

This is equivalent to applying the $N \times d$ transposed Vandermonde matrix of the $\alpha_i$ to $w$. Fast algorithms in arithmetic complexity $\tilde{O}(\max(N, d))$ for this problem are known
in computer algebra (see, e.g., [4]). When executed in FP arithmetic, however, a factor of $O(2^{\min(N,d)})$ in the rounding error bound is responsible for the loss of roughly $\min(N, d)$ bits of accuracy, making this approach impractical for large $N, d$ but with a reasonable FP precision.

This problem is the transpose of the well-known multipoint numerical evaluation problem, for which Moroz [20] recently proposed a quasi-optimal algorithm in the input/output size. Given $P \in \mathbb{R}_{<N}[z]$ a polynomial of degree less than $N$ and $d$ points $\alpha_i$ in the complex unit disk (which we here simplify to $[0, 1]$), this algorithm evaluates $P(\alpha_i)$ with absolute error $O(2^{-\tau})$ using only $\tilde{O}(\max(N, d))$ FP operations in precision $t = \tilde{O}(\tau)$. The main idea is to “zoom in” on a number $K = \lceil \log_2 (\frac{3eN}{\tau}) \rceil$ of overlapping intervals $[\gamma_k - \rho_k, \gamma_k + \rho_k]$ covering $[0, 1]$, where:

$$ 
\begin{align} 
r_k &= 1 - 2^{-k} \quad (0 \le k < K), \qquad r_K = 1, \nonumber \\ 
\gamma_k &= \frac{1}{2}(r_{k-1} + r_k) = \begin{cases} 1 - 3 \cdot 2^{-k-1} & (1 \le k < K), \\ 1 - 2^{-K} & (k = K), \end{cases} \tag{15} \\ 
\rho_k &= \frac{3}{4}(r_{k-1} - r_k) = \begin{cases} 3 \cdot 2^{-k-2} & (1 \le k < K), \\ 3 \cdot 2^{-K-1} & (k = K). \end{cases} \nonumber 
\end{align} 
$$

(1) Shift and truncate $P$ as $Q_k(z) = P(\gamma_k + \rho_k z) \% z^\tau$ ($1 \le k \le K$) using a shift algorithm in $\tilde{O}(N)$ FP operations in precision $t$;

(2) Evaluate $Q_k$ on $\beta_i = \frac{\alpha_i - \gamma_k}{\rho_k}$ for $i \in \mathcal{I}_k = \{i \mid r_{k-1} < \alpha_i \le r_k\}$, using a multipoint evaluation algorithm as in [4] in $\tilde{O}(\max(\#\mathcal{I}_k, \tau))$ FP operations in precision $t$.

Following Tellegen's principle [4], we transpose this scheme:

(2)' For $1 \le k \le K$, compute $V_k(z) = \sum_{i \in \mathcal{I}_k} \frac{w_i}{1-\beta_i z} \% z^\tau$ using an algorithm VANDERT$_0((\beta_i)_{i \in \mathcal{I}_k}, (w_i)_{i \in \mathcal{I}_k}, m)$ like [4, §6.2] in $\tilde{O}(\max(\#\mathcal{I}_k, m))$ FP operations in precision $t$.

(1)' Compute the transposed shift $U_k(z) = \text{shift}^T(V_k, \gamma_k, \rho_k, N)$, and return the sum of the $U_k$.

The transposition of $\text{shift}(-, \gamma, \rho, N) : \mathbb{R}[z]_{<N} \to \mathbb{R}[z]_{<N}, P(z) \mapsto P(\gamma + \rho z)$, maps $\frac{1}{1-\beta z} \% z^N = \sum_{n<N} \beta^n z^n$ to $\frac{1}{1-\alpha z} \% z^N = \sum_{n<N} \alpha^n z^n$ where $\alpha = \gamma + \rho\beta$. It is fully characterized by:

$$ \text{shift}^T(U(z), \gamma, \rho, N) = \frac{1}{1-\gamma z} U\left(\frac{\rho z}{1-\gamma z}\right) \% z^N, \quad (16) $$

$$ = \mathcal{B}^{-1}(U(\rho z) \cdot \exp(\gamma z)) \% z^N, \quad (17) $$

with $\mathcal{B} : P(z) = \sum_n a_n z^n \mapsto \mathcal{B}P(z) = \sum_n \frac{a_n}{n!} z^n$ the Borel transform. This leads to Algorithm VANDERT, with bounds on the complexity and requested FP precision given by the following Lemma.
**Algorithm 3: VANDERT(α, w, N, τ)**

**input:** α = (αᵢ)ᵢ=₁ᵢ ∈ [0, 1]ᵀ, w = (wᵢ)ᵢ=₁ᵢ ∈ ℝ, N, τ ∈ N>0.

**output:** U(z) ∈ ℝ[z]≤N approximating (14).

1 U(z) ← 0; (r_k)_{k=0}^{K}, (\gamma_k)_{k=1}^{K}, (\rho_k)_{k=1}^{K} \leftarrow \text{Eq. (15)}

2 for k = 1...K do

3 $\beta_i \leftarrow \frac{\alpha_i - \gamma_k}{\rho_k}$ **for** $i \in \mathcal{I}_k = \{i \mid r_{k-1} < \alpha_i \le r_k\}$

4 $V(z) \leftarrow \text{VANDERT}_0((\beta_i)_{i \in \mathcal{I}_k}, (w_i)_{i \in \mathcal{I}_k}, \tau)$

5 $U(z) \leftarrow U(z) + \mathcal{B}^{-1}[\mathcal{B}V(\rho_k z) \cdot \exp(\gamma_k z)] \% z^N$

6 end

7 return U(z)

**Proposition 4.2.** Let α₁, ..., α₉ ∈ [0, 1] and U = ∑ᵢ wᵢ (1 - αᵢz) % zⁿ.

Then VANDERT(α, w, N, τ) computes $\tilde{U}$ s.t. $\|\tilde{U} - U\|_1 \le \|w\|_1 N 2^{-\tau}$ using $\tilde{O}(N + d)$ FP operations in precision $t = \tilde{O}(\tau)$.

*Proof.* The claimed arithmetic complexity follows from that of VANDERT₀ [4, §6.2] and
FFT polynomial multiplication. It remains to show the error bound using a working
FP precision of $t = \tilde{O}(\tau)$.

We first prove that $\|U - U^*\|_1 \le \|w\|_1 N 3^{-\tau}$ where $U^*$ would be the result of VANDERT(α, w, N, τ) without rounding errors, i.e., $U^* = \sum_k \text{shift}^T(V_k, \gamma_k, \rho_k, N)$ with $V_k = \sum_{i \in \mathcal{I}_k} \frac{w_i}{1-\beta_i z} \% z^\tau$. It is sufficient for this to observe how shift$^T$ acts on a truncated tail $E_i = \sum_{n=\tau}^{N-1} \beta_i^n z^n = \frac{\beta_i^\tau z^\tau}{1-\beta_i z} \% z^\tau$. Using Eq. (16),

$$
E_i^* := \text{shift}^T(E_i, \gamma_k, \rho_k, N) = \frac{\left(\beta_i \frac{\rho_k z}{1-\gamma_k z}\right)^\tau}{1 - \alpha_i z} \% z^N.
$$

First, use the bound $\|(1 - \alpha_i z)^{-1} \% z^N\|_1 \le N$. If $k < K$, then $\frac{\rho_k}{1 - \gamma_k} = \frac{1}{2}$ and so
$$ \| (\beta_i \rho_k z)^{\tau} (1 - \gamma_k z)^{-\tau} \% z^N \|_1 \le \frac{|(\beta_i \rho_k)^{\tau}|}{(1 - \gamma_k)^{\tau}} \le 3^{-\tau} $$
Hence, $\|E_i^*\|_1 \le N 3^{-\tau}$. If $k = K$, $\rho_k = \frac{3}{2} 2^{-K}$, $|\beta_i| \rho_k \le 2^{-K}$ and $1 - \gamma_k = 2^{-K}$. Using the identity $(1 - \gamma_k z)^{-\tau} = \sum_{n \ge 0} \binom{n+\tau-1}{\tau-1} \gamma_k^n z^n$,
we bound $\|z^\tau (1 - \gamma_k z)^{-\tau} \% z^N\|_1 \le \sum_{n=\tau-1}^{N-2} \binom{n}{\tau-1} = \binom{N-1}{\tau}$. This gives $\|E_i^*\|_1 \le N \binom{N-1}{\tau} 2^{-K\tau} \le N \binom{Ne2^{-K}}{\tau} \le N 3^{-\tau}$ by definition of K.

We finally prove that $t = \tilde{O}(\tau)$ is a sufficient working FP precision to guarantee
the claimed error bound. We assumed that VANDERT₀ computes $V_k$ with error $2^{-\tau}$
using $\tilde{O}(\tau)$ bits of precision. For shift$^T$ to also require no more than $\tilde{O}(\tau)$ bits of
precision, the product and $\mathcal{B}/\mathcal{B}^{-1}$ transforms in line 5 must be implemented carefully
to avoid scalings with $n!$ ($n < N$) that would cause the loss of $\tilde{O}(N)$ bits. Decompose
$\exp(\gamma_k z)$ as $\sum \frac{A_j(z)}{((\tau-1)j)!} z^{(\tau-1)j}$ with $\deg(A_j) \le \tau - 2$, and compute separately each
product $\mathcal{B}V_k(\rho_k z)^j \cdot A_j(z) z^{(h-1)j}$ (with error bound given by Lemma 4.1), followed by
$\mathcal{B}^{-1}$. By partially cancelling the factorials, remaining scaling factors are brought down to a size of $N^{O(\tau)}$, thus requiring only $\tilde{O}(\tau)$ extra bits of precision. $\square$

## 4.2 FFT-based exponentiation and final rounding error bound

We explain step by step how to combine the transposed multipoint evaluation routine VANDERT with FFT/iFFT-based exponentiation of a power series into Algorithm CDFFAST. We carefully derive asymptotic error bounds for each step, eventually establishing the complexity stated in Theorem 1.2 for a given absolute target error.

Let again $N$ denote the number of terms $c_n$ to be computed in the power series of $f(\xi)$, $M = \sum_{i=1}^d m_i$ and $\alpha_i = p_1 - p_i < p_1$. Define

$$ \hat{f}^*(\lambda) = (1 - p_1\lambda)\hat{f}(\lambda) = \left(\prod_{i=1}^d \sqrt{p_i}\right) \exp(-M) \exp(\psi(\lambda)), $$

$$ \text{where } \psi(\lambda) = p_1 m_1 \lambda + \sum_{i=2}^{d} \left( \frac{p_i m_i \lambda}{1 - \alpha_i \lambda} - \frac{1}{2} \ln(1 - \alpha_i \lambda) \right), $$

so that $\hat{f}^*(\lambda)$ and $\psi(\lambda)$ are analytic over the closed disk $|\lambda| \le \frac{1}{p_1}$. Note that $\psi(0) = 0$, $\psi(\lambda) \gg 0$, and $\hat{f}^*(\lambda) \gg 0$. Since $\hat{f}(\lambda)$ has a pole at $\frac{1}{p_1}$, its coefficients $\hat{c}_n$ grow like $p_1^n$. We therefore need to rescale them appropriately to prevent the large degree coefficients to “pollute” the low degree ones in FFT/iFFT operations as predicted by the error bounds in Lemma 4.1. For this, we use the change of variable $\lambda = \frac{\nu}{p_1}z$ with $\nu = 2^{-t/N}$ ($\nu$ close to but smaller than 1), and we define $\hat{f}_\nu(z) = \hat{f}(\lambda) = \sum_n \hat{c}_{\nu,n} z^n$, $\hat{f}_\nu^*(z) = \hat{f}^*(\lambda) = \sum_n \hat{c}_{\nu,n}^* z^n$, and $\psi_\nu(z) = \psi(\lambda)$:

$$ \psi_\nu(z) = \nu m_1 z + z \sum_{i=2}^{d} \frac{w_i}{1 - \alpha_i' z} + \int_0 \sum_{i=2}^{d} \frac{w_i'}{1 - \alpha_i' z} dz, \quad (18) $$

$$ \alpha_i' = \frac{\nu}{p_1} \alpha_i, \quad w_i = \frac{\nu}{p_1} p_i m_i, \quad w_i' = \frac{\alpha_i'}{2}, \quad (2 \le i \le d), $$

where $\int_0 dz$ denotes the formal antiderivative w.r.t. $z$ with a zero constant coefficient. Call $\psi_{\nu,N}(z) = \psi_\nu(z) \% z^N$ its truncation to degree $N-1$. Since $\frac{p_i}{1-\alpha_i \lambda} = p_1$ when $\lambda = \frac{1}{p_1}$, we have

$$ \|\psi_{\nu,N}\|_1 \le \psi_\nu(1) \le \psi\left(\frac{1}{p_1}\right) = M + \frac{1}{2} \sum_{i=1}^{d} \ln \frac{p_1}{p_i}. \quad (19) $$

Using that $\psi_{\nu,N}$ is truncated, that $\|w'\|_1 \le \frac{d-1}{2}$ and that $\sum_{n=1}^{N-1} \frac{1}{n} \le \ln(N) + 1$, yields an alternative bound independent of the ratio $\frac{p_i}{p_1}$,

$$ \|\psi_{\nu,N}\|_1 \le M + \frac{d-1}{2}(\ln(N) + 1). \quad (20) $$
In lines 2–4, $\psi_{\nu,N}$ is approximated as $\tilde{\psi}_{\nu,N}$ using VANDERT twice using a higher precision $t' = \tilde{O}(t)$.

**Lemma 4.3.** $\|\tilde{\psi}_{\nu,N} - \psi_{\nu,N}\|_1 = O((M+d)Nu)$.

*Proof.* This follows from Prop. 4.2 with $\|w\|_1 = M - m$ and $\|w'\|_1 \le \frac{d-1}{2}$, and taking into account the additional rounding errors when taking the antiderivative and summing the three polynomials in line 4. $\square$

We are now ready for the FFT/iFFT-based exponentiation of $\psi_\nu$, since $\hat{f}_{\nu,N}^*(z) := \hat{f}_\nu^*(z)\%z^N = \hat{c}_0 \exp(\psi_{\nu,N}(z))\%z^N$ and $\tilde{\psi}_{\nu,N} \approx \psi_{\nu,N}$. First, $\tilde{v} = (\tilde{v}_j)_{j=0}^{N-1} = \text{FFT}(\tilde{\psi}_{\nu,N}, N)$ is computed in line 5 to approximate $v = (v_j)_{j=0}^{N-1} = \text{DFT}(\psi_{\nu,N}, N) = (\psi_{\nu,N}(\omega^j))_{j=0}^{N-1}$ where $\omega = \exp(\frac{2i\pi}{N})$.

**Lemma 4.4.** $\|\tilde{v} - v\|_\infty = O((M+d)Nu)$.

*Proof.* The bound is the sum of two terms:

* the DFT of the input error, which is bounded by $\|\tilde{\psi}_{\nu,N} - \psi_{\nu,N}\|_1 = O((M+d)Nu)$ by Lemma 4.3;

* the FFT error, which by Lemma 4.1 and Eq. (20) is

$$O(\|\tilde{\psi}_{\nu,N}\|_1 \ln(N)u) = O(\|\psi_{\nu,N}\|_1 \ln(N)u) = O((M+d) \ln(N)^2 u).$$

<figure><img src="image_2.png" /></figure>

Now, line 6 computes approximations $\tilde{e}_j = \hat{c}_0 \exp(\tilde{v}_j)$ of the values $e_j = \hat{c}_0 \exp(v_j) = \hat{c}_0 \exp(\psi_{\nu,N}(\omega^j))$.

**Lemma 4.5.** $\|e\|_\infty \le p_1^{d/2}$ and $\|\tilde{e} - e\|_\infty = O(p_1^{d/2}(M+d)Nu)$.

*Proof.* We have $|e_j| \le \hat{c}_0 \exp(\|\psi_{\nu,N}\|_1) = p_1^{d/2}$ by Eq. (19), and

$$|\tilde{e}_j - e_j| \le |e_j|(\exp(|\tilde{v}_j - v_j|) - 1) \le \|e\|_\infty(\exp(|\tilde{v} - v|_\infty) - 1) \\ = O(\|e\|_\infty \|\tilde{v} - v\|_\infty) = O(p_1^{d/2}(M+d)Nu),$$

plus the rounding errors for $\hat{c}_0$, exp and multiplication. $\square$

<figure><img src="image_3.png" /></figure>

We eventually apply the inverse FFT onto the $\tilde{e}_j$ to approximate $\hat{f}_{\nu,N}^*(z)$ by $\tilde{f}_{\nu,N}^* = \text{iFFT}(\tilde{e}, N)$ in line 7.

**Lemma 4.6.** $\|\tilde{f}_{\nu,N}^* - \hat{f}_{\nu,N}^*\|_1 = O(p_1^{d/2}(M+d)N^2u)$.

*Proof.* The bound is the sum of:

* the iDFT of the input error, i.e., $N\|\tilde{e} - e\|_\infty = O(p_1^{d/2}(M+d)N^2u)$ by Lemma 4.5;
• the iFFT error, which using Lemmas 4.1 and 4.5 is

$$O(\|e\|_{\infty} N \ln(N) u) = O(p_1^{d/2} N \ln(N) u);$$

• the aliasing error: Calling $\tilde{c}_{\nu,n}$ the (nonnegative) coefficients of $\hat{c}_0 \exp(\psi_{\nu,N}(z))$, the degree $n$ coefficient of iDFT$(e, N)$ is not $\tilde{c}_{\nu,n}$ but $\sum_{j \ge 0} \tilde{c}_{\nu,n+jN}$. The resulting error is therefore bounded by $\sum_{n > N} \tilde{c}_{\nu,n}$, and hence by $\sum_{n > N} \hat{c}_{\nu,n}^*$ since $\hat{c}_0 \exp(\psi_{\nu,N}(z)) \ll \hat{c}_0 \exp(\psi_{\nu}(z)) = \hat{f}_{\nu}^*(z)$. Because $\hat{f}_{\nu}^*(z)$ is analytic and bounded for $|z| \le \frac{1}{\nu}$ by $\hat{f}_{\nu}^*(\frac{1}{\nu}) = \hat{f}^*(\frac{1}{p_1}) = p_1^{d/2}$, the Cauchy bound yields $\hat{c}_{\nu,n}^* \le p_1^{d/2} \nu^n$. Therefore, the aliasing error is bounded by $p_1^{d/2} \frac{\nu^N}{1-\nu} = O(p_1^{d/2} Nu)$, since $\nu^N = u$ and $1 - \nu \sim \frac{t}{N} \ln(2)$. $\square$

The above proof shows that $\nu$ should not be too close to 1 to guarantee a low aliasing error. Nor should it be too far in what follows. Since $\hat{f}_{\nu}(\lambda) = \frac{1}{1-\nu z} \hat{f}_{\nu}^*(\lambda)$, we have the recurrence relation

$$\hat{c}_{\nu,n} = \sum_{j \le n} \nu^{n-j} \hat{c}_{\nu,j}^* = \nu \hat{c}_{\nu,n-1} + \hat{c}_{\nu,n}^*.$$

Let $\tilde{c}_{\nu,n}$ be the approximations of $\hat{c}_{\nu,n}$ computed numerically from the $\tilde{c}_{\nu,n}^*$ using this recurrence. The sum of the error due to the difference between the $\tilde{c}_{\nu,n}^*$ and $\hat{c}_{\nu,n}^*$ (bounded by $\|\tilde{f}_{\nu,N}^* - \hat{f}_{\nu,N}^*\|_1$) plus the one resulting from the $2N$ arithmetic operations to unroll the recurrence (bounded by $O(\|\tilde{f}_{\nu,N}^*\|_1 Nu)$) yields

$$|\tilde{c}_{\nu,n} - \hat{c}_{\nu,n}| = O(p_1^{d/2} (M+d) N^2 u).$$

The last instructions in lines 9–11 eventually compute the expected approximation $\tilde{s}_N$ of $s_N$,

$$s_N = e^{-p_1} \sum_{n=0}^{N-1} \frac{\hat{c}_{\nu,n} p_1^n}{\Gamma(n + \frac{d}{2} + 1) \nu^n}, \quad \tilde{s}_N = e^{-p_1} \sum_{n=0}^{N-1} \frac{\tilde{c}_{\nu,n} p_1^n}{\Gamma(n + \frac{d}{2} + 1) \nu^n}.$$

**Proposition 4.7.** $|\tilde{s}_N - s_N| = O((M+d)N^2 2^{-\frac{t}{2}})$.

*Proof.* Neglecting rounding errors to compute $\tilde{s}_N$ from the $\tilde{c}_{\nu,n}$ ($O(N)$ arithmetic operations on positive FP numbers), we have

$$\begin{aligned} |\tilde{s}_N - s_N| & \le \|\tilde{c}_\nu - \hat{c}_\nu\|_{\infty} e^{-p_1} \sum_{n=0}^{N-1} \frac{(p_1/\nu)^n}{\Gamma(n + \frac{d}{2} + 1) \nu^n}, \\ & \le \|\tilde{c}_\nu - \hat{c}_\nu\|_{\infty} p_1^{-\frac{d}{2}} e^{-p_1} \sum_{n=0}^{N-1} \frac{(p_1/\nu)^{n+\frac{d}{2}}}{\Gamma(n + \frac{d}{2} + 1) \nu^n} \le \|\tilde{c}_\nu - \hat{c}_\nu\|_{\infty} p_1^{-\frac{d}{2}} \exp\left(\frac{p_1}{\nu} - p_1\right), \end{aligned}$$
**Algorithm 4:** CDFFAST(d, (pᵢ)ᵢ=₁, (mᵢ)ᵢ=₁, N, t)

**input:** d ∈ N>0; (pᵢ)ᵢ=₁, (mᵢ)ᵢ=₁ ∈ (ℝ>0)ᵈ; N, t ∈ N>0

**output:** a numerical approximation s of g(1).

\# use FP precision t unless otherwise specified

1. $\nu \leftarrow 2^{-t/N}$; $\alpha', w, w'$ as in Eq. (18)

2. $S_1(z) \leftarrow \text{VANDERT}(\alpha', w, N - 1, t)$ \quad \# use FP prec. $t' = \tilde{O}(t)$

3. $S_2(z) \leftarrow \text{VANDERT}(\alpha', w', N - 1, t)$ \quad \# idem

4. $\psi_{\nu,N}(z) \leftarrow z(\nu m_1 + S_1(z)) + \int_0 S_2(z) dz$

5. $v = (v_j)_{j=0}^{N-1} \leftarrow \text{FFT}(\psi_{\nu,N}(z), N)$

6. $e = (e_j)_{j=0}^{N-1} \leftarrow (\hat{c}_0 \exp(v_j))_{j=0}^{N-1}$ \quad \# $\hat{c}_0 = \left( \prod_{i=1}^d \sqrt{p_i} \right) \exp(-M)$

7. $\hat{f}_{\nu,N}^*(z) = \sum_{n=0}^{N-1} \hat{c}_{\nu,n}^* z^n \leftarrow \text{iFFT}(e, N)$

8. $\hat{c}_{\nu,0} \leftarrow \hat{c}_{\nu,0}^*$; **for** n = 1...N-1 **do** $\hat{c}_{\nu,n} \leftarrow \nu \hat{c}_{\nu,n-1} + \hat{c}_{\nu,n}^*$

9. $r \leftarrow 1/\Gamma(\frac{d}{2} + 1)$; $s \leftarrow 0$

10. **for** n = 0...N-1 **do** s ← s + r$\hat{c}_{\nu,n}$; $r \leftarrow \frac{p_1}{(n+d/2+1)\nu}r$

11. **return** s ← $\exp(-p_1)s$

where we used the identities $\sum_{n \ge 0} \frac{x^n}{\Gamma(n+1)} = e^x$ and $\sum_{n \ge 0} \frac{x^{n+\frac{1}{2}}}{\Gamma(n+\frac{3}{2})} = e^x \operatorname{erf}(\sqrt{x}) \le e^x$. The
choice of $\nu = 2^{-t/N}$ implies

$$
\frac{p_1}{\nu} - p_1 = \frac{1 - \nu}{\nu} p_1 = \frac{tp_1}{N} \ln(2) + o(t) \le \frac{t}{2} \ln(2),
$$

assuming $N$ was chosen greater than $2p_1$, as per Prop. 2.2(i). $\square$

*Proof of Theorem 1.2.* According to Prop. 4.7, choose a working precision $t = 2\tau + O(\ln(N) + \ln(M+d))$ for CDFFAST and $t' = \tilde{O}(t)$ for the two calls to VANDERT, to have a final absolute error $r_N \le 2^{-\tau}$. The arithmetic complexity is dominated by that of VANDERT and degree $N$ FFT/iFFT, hence in $\tilde{O}(N+d)$, with $N = O(p_1 + \tau)$ by Prop. 2.2(i). $\square$

5 Implementation and numerical experiments

The Julia and Matlab® implementation of Algorithms CDFREC and CDFPosREC is
available in [1]. We present an accuracy and performance comparison with the existing
state of the art routine for computing the generalized chi-square upper-tail cdf of the
Matlab toolbox 'Generalized chi-square distribution' of A. Das [8].

The function `gx2cdf` of this toolbox incorporates methods based on the numerical inversion of the characteristic function (`imhof`, `ifft`), on a convergent series expansion (`Ruben`) and on ray tracing (`ray`). Due to its very poor performance, this last one has not been included in the comparisons. A reference value is also computed in Julia
by CDFPOSREC with $N = 10000$ terms and extra precision (1024 bits). In [8], 16 test cases are presented, borrowed from previously published works [13] (cases #1 to #12) and [16] (cases #13 to #16). The cases #8, #10 and #12 are excluded from our study, since their quadratic form of Gaussian variables is not positive semidefinite and algorithms CDFREC and CDFPOSREC do not apply in these cases. For each remaining case, the cdf of the quadratic form is evaluated at three different values. The specific parameter values for each case are given in Appendix A, and we only report the dimension and the ratio $p_1/\min_i p_i$ below (note that even if these quantities are the same for several cases, other omitted parameters are different).

For each evaluation, the algorithms based on series expansions (CDFREC, CDF-POSREC and Ruben) are run with a sufficient number of terms $N^*$ (estimated by our bounds) in order to obtain a target accuracy of $1e-07$. The parameters of the two other tested methods (imhof, ifft) are the ones given in Table 2 of the reference [8].

For each tested method, the relative error wrt. the reference value is given; timing ratios are given wrt. the reference timing $t_r$ of the CDFPOSREC algorithm and were obtained by repeated averaged simulations. These were performed on Matlab R2025b® with a 12th Gen Intel® Core(™) i7-1265U (1.80 GHz) and 32 GB of RAM. All the results are gathered in Table 1.

They show that CDFPOSREC is overall the most efficient, keeping in mind the tradeoff accuracy/timing and that Ruben has similar results in terms of accuracy but is always slower with a factor of order 10 for some examples. Algorithm CDFREC is in general the fastest but has some numerical difficulties on particular examples (cases #7 and #11). This is most probably due to some cancellations when evaluating the coefficients of the series in double precision, since a Julia evaluation with 1024 bits does not have this issue. Algorithms Imhof and ifft are less accurate and slower (in particular for the second one with a factor of order $10^5$ for some examples). For each method, the bold values highlight the cases where it is slower or less accurate.

# 6 Conclusion and perspectives

We revisited the finite-precision evaluation of the cdf of positive semidefinite quadratic forms in Gaussian random variables through a symbolic-numeric lens based on modified Laplace transforms and D-finite series. Starting from a closed-form expression of the modified Laplace transform, we reduced cdf evaluation to summing a truncated power series with provably nonnegative coefficients and explicit truncation bounds. While this already yields a reliable scheme, we highlighted that the classical scalar recurrence for generating the coefficients may still suffer from finite-precision cancellation as the dimension grows.

To address this, we introduced a coupled recurrence (CDFPOSREC) whose update coefficients are all nonnegative. This “positivity throughout” property substantially improves robustness and enables a sharper forward-error analysis, leading to the relative-accuracy complexity bound stated in Theorem 1.1. We also proposed a second approach (CDFFAST) aimed at reducing the arithmetic cost of coefficient generation
<table><caption>Table 1: Numerical results comparisons of different methods on a set of generalized chi-square distributions from [8].</caption><thead><tr><th rowspan="2">Case</th><th rowspan="2">d</th><th rowspan="2">p<sub>1</sub>/min, p<sub>1</sub></th><th rowspan="2">Ref. Val.</th><th rowspan="2">N*</th><th colspan="2">CDFPOSREC</th><th colspan="2">CDFREC</th><th colspan="2">Ruben</th><th colspan="2">Imhof</th><th colspan="2">IFFT</th></tr><tr><th>Rel. Err.</th><th>t<sub>r</sub> (s)</th><th>Rel. Err.</th><th>t/t<sub>r</sub></th><th>Rel. Err.</th><th>t/t<sub>r</sub></th><th>Rel. Err.</th><th>t/t<sub>r</sub></th><th>Rel. Err.</th><th>t/t<sub>r</sub></th></tr></thead><tbody><tr><td rowspan="3">1</td><td rowspan="3">3</td><td rowspan="3">2.45</td><td>0.94578615</td><td>17</td><td>3.1e-12</td><td>5.87e-04</td><td>3.1e-12</td><td>6.03e-01</td><td>3.1e-12</td><td>7.59e+00</td><td>1.7e-09</td><td>2.33e+01</td><td>2.5e-05</td><td>7.50e+02</td></tr><tr><td>0.50643823</td><td>33</td><td>3.7e-13</td><td>2.57e-04</td><td>3.7e-13</td><td>9.05e-01</td><td>3.7e-13</td><td><b>2.56e+01</b></td><td>1.0e-09</td><td><b>1.20e+02</b></td><td>2.4e-05</td><td>7.90e+02</td></tr><tr><td>0.12395907</td><td>95</td><td>1.9e-11</td><td>1.16e-03</td><td>1.9e-11</td><td>2.16e-01</td><td>1.9e-11</td><td>4.88e+00</td><td>5.0e-09</td><td>9.15e+01</td><td>4.0e-04</td><td>2.38e+02</td></tr><tr><td rowspan="3">2</td><td rowspan="3">6</td><td rowspan="3">2.45</td><td>0.99354712</td><td>14</td><td>4.0e-12</td><td>4.84e-04</td><td>4.0e-12</td><td>1.26e+00</td><td>4.0e-12</td><td>9.28e+00</td><td>3.5e-09</td><td>1.11e+01</td><td>9.0e-07</td><td><b>1.13e+05</b></td></tr><tr><td>0.399795</td><td>108</td><td>5.6e-12</td><td>2.38e-03</td><td>7.5e-12</td><td>2.77e-01</td><td>5.6e-12</td><td>5.58e+00</td><td>7.0e-09</td><td>3.53e+00</td><td>1.5e-10</td><td>2.22e+04</td></tr><tr><td>0.016102973</td><td>326</td><td>1.1e-11</td><td>7.08e-03</td><td>1.2e-08</td><td>1.17e-01</td><td>1.1e-11</td><td>3.59e+00</td><td>6.2e-07</td><td>9.94e-01</td><td>1.2e-04</td><td>7.69e+03</td></tr><tr><td rowspan="3">3</td><td rowspan="3">12</td><td rowspan="3">2.45</td><td>0.99731927</td><td>54</td><td>1.8e-12</td><td>3.41e-04</td><td>1.8e-12</td><td>1.03e+00</td><td>1.8e-12</td><td>7.27e+00</td><td>3.5e-11</td><td>7.47e+00</td><td>2.0e-05</td><td>4.66e+02</td></tr><tr><td>0.43525063</td><td>271</td><td>3.0e-12</td><td>6.31e-03</td><td>2.0e-05</td><td>1.58e-01</td><td>3.0e-12</td><td>3.61e+00</td><td>8.8e-10</td><td>6.60e-01</td><td>3.9e-09</td><td>5.22e+01</td></tr><tr><td>0.0087690053</td><td>652</td><td>3.4e-12</td><td>1.80e-02</td><td><b>5.2e+01</b></td><td>7.53e-02</td><td>3.1e-12</td><td>3.25e+00</td><td>1.4e-07</td><td>2.51e-01</td><td>4.0e-03</td><td>1.85e+01</td></tr><tr><td rowspan="3">4</td><td rowspan="3">12</td><td rowspan="3">2.45</td><td>0.96664038</td><td>54</td><td>4.6e-12</td><td>3.56e-04</td><td>4.6e-12</td><td>7.93e-01</td><td>4.6e-12</td><td>7.56e+00</td><td>1.8e-09</td><td>8.37e+00</td><td>1.0e-05</td><td>3.90e+02</td></tr><tr><td>0.41955462</td><td>163</td><td>5.5e-12</td><td>8.61e-04</td><td>3.0e-11</td><td>3.32e-01</td><td>5.5e-12</td><td>4.79e+00</td><td>1.2e-09</td><td>3.25e+00</td><td>1.3e-07</td><td>1.57e+02</td></tr><tr><td>0.0087153638</td><td>434</td><td>5.2e-12</td><td>2.41e-03</td><td>1.2e-05</td><td>2.05e-01</td><td>5.3e-12</td><td>3.76e+00</td><td>2.8e-07</td><td>1.21e+00</td><td>2.9e-03</td><td>5.70e+01</td></tr><tr><td rowspan="3">5</td><td rowspan="3">8</td><td rowspan="3">1.53</td><td>0.99388203</td><td>36</td><td>3.6e-12</td><td>7.22e-04</td><td>3.6e-12</td><td>8.03e-01</td><td>3.6e-12</td><td>5.91e+00</td><td>1.1e-08</td><td>4.73e+00</td><td>3.8e-06</td><td>4.91e+03</td></tr><tr><td>0.40865788</td><td>181</td><td>7.7e-12</td><td>2.78e-03</td><td>4.7e-13</td><td>1.98e-01</td><td>7.7e-12</td><td>3.77e+00</td><td>7.8e-09</td><td>1.24e+00</td><td>4.9e-07</td><td>1.14e+03</td></tr><tr><td>0.022081647</td><td>362</td><td>1.6e-11</td><td>6.54e-03</td><td>1.1e-08</td><td>1.07e-01</td><td>1.6e-11</td><td>3.41e+00</td><td>1.3e-07</td><td>6.56e-01</td><td>2.4e-04</td><td>5.26e+02</td></tr><tr><td rowspan="3">6</td><td rowspan="3">2</td><td rowspan="3">1.53</td><td>0.95487281</td><td>18</td><td>2.5e-12</td><td>4.15e-04</td><td>2.5e-12</td><td>6.46e-01</td><td>2.5e-12</td><td>9.62e+00</td><td>4.2e-08</td><td>6.05e+01</td><td>5.1e-06</td><td>9.51e+03</td></tr><tr><td>0.40756543</td><td>108</td><td>2.5e-12</td><td>5.21e-04</td><td>2.5e-12</td><td>3.18e-01</td><td>2.5e-12</td><td>7.13e+00</td><td>4.6e-09</td><td><b>2.23e+02</b></td><td>1.7e-06</td><td>4.44e+03</td></tr><tr><td>0.022343129</td><td>271</td><td>1.3e-11</td><td>1.20e-03</td><td>1.3e-11</td><td>1.59e-01</td><td>1.3e-11</td><td>5.28e+00</td><td>4.9e-08</td><td><b>1.30e+02</b></td><td>2.2e-04</td><td>1.93e+03</td></tr><tr><td rowspan="3">7</td><td rowspan="3">24</td><td rowspan="3">3.46</td><td>0.98905831</td><td>244</td><td>3.4e-13</td><td>2.22e-03</td><td><b>4.6e-03</b></td><td>3.83e-01</td><td>3.4e-13</td><td>3.79e+00</td><td>1.4e-11</td><td>1.37e+00</td><td>1.1e-05</td><td>1.46e+03</td></tr><tr><td>0.34526541</td><td>652</td><td>9.5e-12</td><td>3.72e-03</td><td><b>6.0e+11</b></td><td>2.08e-01</td><td>9.5e-12</td><td>4.07e+00</td><td>3.2e-10</td><td>1.46e+00</td><td>4.8e-06</td><td>7.38e+02</td></tr><tr><td>0.015399638</td><td>1141</td><td>2.1e-11</td><td>6.63e-03</td><td><b>1.3e+28</b></td><td>1.39e-01</td><td>2.1e-11</td><td>3.32e+00</td><td>1.1e-08</td><td>4.99e-01</td><td>1.1e-03</td><td>4.20e+02</td></tr><tr><td rowspan="3">9</td><td rowspan="3">10</td><td rowspan="3">1.53</td><td>0.95631841</td><td>126</td><td>8.7e-13</td><td>2.88e-03</td><td>8.6e-13</td><td>2.89e-01</td><td>8.7e-13</td><td>3.97e+00</td><td>1.3e-09</td><td>1.57e+00</td><td>2.2e-06</td><td>2.01e+04</td></tr><tr><td>0.41523898</td><td>289</td><td>4.8e-12</td><td>6.97e-03</td><td>1.1e-09</td><td>1.75e-01</td><td>4.8e-12</td><td>3.69e+00</td><td>1.5e-08</td><td>7.64e-01</td><td>3.6e-07</td><td>8.33e+03</td></tr><tr><td>0.046230859</td><td>471</td><td>3.7e-13</td><td>5.65e-03</td><td>3.3e-07</td><td>1.25e-01</td><td>2.8e-13</td><td>3.61e+00</td><td>2.3e-08</td><td>7.25e-01</td><td>5.7e-05</td><td>7.57e+03</td></tr><tr><td rowspan="3">11</td><td rowspan="3">34</td><td rowspan="3">2.65</td><td>0.98415909</td><td>652</td><td>4.9e-13</td><td>1.46e-02</td><td><b>4.8e+13</b></td><td>1.72e-01</td><td>4.9e-13</td><td>4.00e+00</td><td>1.4e-12</td><td>5.08e-01</td><td>1.4e-06</td><td>3.70e+02</td></tr><tr><td>0.42637747</td><td>1304</td><td>7.8e-12</td><td>4.08e-02</td><td><b>7.5e+44</b></td><td>9.55e-02</td><td>7.8e-12</td><td>3.23e+00</td><td>1.3e-10</td><td>1.50e-01</td><td>1.8e-07</td><td>1.69e+02</td></tr><tr><td>0.011662614</td><td>2174</td><td>4.2e-11</td><td>2.95e-02</td><td><b>2.6e+87</b></td><td>7.95e-02</td><td>4.1e-11</td><td>3.35e+00</td><td>2.3e-09</td><td>1.52e-01</td><td>1.8e-04</td><td>1.60e+02</td></tr><tr><td rowspan="3">13</td><td rowspan="3">4</td><td rowspan="3">2.24</td><td>0.45746062</td><td>108</td><td>3.2e-12</td><td>4.36e-03</td><td>3.2e-12</td><td>8.68e-01</td><td>3.2e-12</td><td>8.03e+00</td><td>1.8e-09</td><td>1.62e+01</td><td>6.6e-07</td><td>3.28e+04</td></tr><tr><td>0.031108931</td><td>326</td><td>1.4e-11</td><td>8.91e-03</td><td>2.0e-10</td><td>4.87e-01</td><td>1.4e-11</td><td>5.78e+00</td><td>2.9e-08</td><td>9.24e+00</td><td>1.2e-05</td><td>1.93e+04</td></tr><tr><td>0.0068853922</td><td>434</td><td>4.3e-12</td><td>1.09e-02</td><td>2.7e-10</td><td>9.89e-01</td><td>5.3e-12</td><td>7.09e+00</td><td>5.1e-07</td><td>1.35e+01</td><td>9.5e-05</td><td>1.23e+04</td></tr><tr><td rowspan="3">14</td><td rowspan="3">2</td><td rowspan="3">1.52</td><td>0.95487281</td><td>18</td><td>2.5e-12</td><td>1.05e-03</td><td>2.5e-12</td><td>2.18e+00</td><td>2.5e-12</td><td><b>2.52e+01</b></td><td>4.2e-08</td><td><b>1.11e+02</b></td><td>5.6e-07</td><td><b>1.31e+05</b></td></tr><tr><td>0.40756543</td><td>108</td><td>2.5e-12</td><td>2.80e-03</td><td>2.5e-12</td><td>8.02e-01</td><td>2.5e-12</td><td><b>1.05e+01</b></td><td>4.6e-09</td><td>8.18e+01</td><td>4.5e-08</td><td>5.08e+04</td></tr><tr><td>0.022343129</td><td>271</td><td>1.3e-11</td><td>5.22e-03</td><td>1.3e-11</td><td>8.70e-01</td><td>1.3e-11</td><td><b>1.03e+01</b></td><td>4.9e-08</td><td>6.21e+01</td><td>5.2e-05</td><td>3.12e+04</td></tr><tr><td rowspan="3">15</td><td rowspan="3">3</td><td rowspan="3">14.1</td><td>0.34793927</td><td>2174</td><td>9.7e-12</td><td>4.86e-02</td><td>2.3e-11</td><td>1.90e-01</td><td>9.7e-12</td><td>3.40e+00</td><td>1.4e-09</td><td>4.02e+00</td><td>2.5e-05</td><td>3.60e+03</td></tr><tr><td>0.033475151</td><td>8698</td><td>1.0e-11</td><td>1.24e-01</td><td>4.5e-08</td><td>9.53e-02</td><td>1.2e-11</td><td>5.75e+00</td><td>3.1e-07</td><td>1.60e+00</td><td>3.8e-05</td><td>1.28e+03</td></tr><tr><td>0.0067478997</td><td>13047</td><td>7.0e-12</td><td>2.14e-01</td><td>1.1e-06</td><td>9.04e-02</td><td>1.9e-12</td><td>4.84e+00</td><td>2.4e-06</td><td>1.38e+00</td><td>1.2e-05</td><td>6.34e+02</td></tr><tr><td rowspan="3">16</td><td rowspan="3">10</td><td rowspan="3">1.53</td><td>0.95631841</td><td>126</td><td>8.7e-13</td><td>5.57e-03</td><td>7.3e-13</td><td>2.14e+00</td><td>8.7e-13</td><td><b>1.99e+01</b></td><td>1.3e-09</td><td>3.43e+01</td><td>2.2e-07</td><td>3.96e+04</td></tr><tr><td>0.41523898</td><td>289</td><td>4.8e-12</td><td>1.40e-02</td><td>1.6e-09</td><td>1.31e+00</td><td>4.8e-12</td><td>7.40e+00</td><td>1.5e-08</td><td>1.62e+01</td><td>3.6e-08</td><td>1.81e+04</td></tr><tr><td>0.046230859</td><td>471</td><td>3.5e-13</td><td>6.10e-03</td><td>1.3e-07</td><td>6.14e-01</td><td>2.8e-13</td><td>4.48e+00</td><td>2.3e-08</td><td>4.42e+00</td><td>5.7e-06</td><td>7.06e+03</td></tr></tbody></table>

to quasi-linear in $\max(N, d)$ by combining transposed multipoint evaluation and FFT-based series exponentiation. Although currently supported only by an absolute-error analysis (Theorem 1.2), it suggests a promising path toward scalable high-dimensional evaluations.

Our numerical experiments confirm that CDFPOSREC achieves an excellent accuracy/time trade-off on standard benchmarks, and that it mitigates the numerical instabilities observed with the alternating-sign recurrence in difficult cases.

Several directions remain open. First, we will implement the fast scheme CDFFAST for absolute-error control and then extend it to relative-error guarantees. Second, we aim to integrate these methods into statistical computing environments and to move beyond the positive semidefinite case, where oscillatory cdfs require different analytic preconditioning. Third, it would be interesting to study the related parameter-dependent optimization problem, for which holonomic gradient descent has been proposed. Finally, we will investigate applying CDFFAST within fast multipole methods.
## References

[1] ARZELIER, D., BRÉHARD, F., AND JOLDES, M. Quadratiques: Fast and reliable evaluation of the distribution of quadratic forms of gaussian random variables, 2026. GitLab repository, accessed 2026-02-06.

[2] ARZELIER, D., BRÉHARD, F., JOLDES, M., AND M., M. Rounding error analysis of an orbital collision probability evaluation algorithm. In *Proceedings - IEEE 31st Symposium on Computer Arithmetic, ARITH 2024* (2024), IEEE, pp. 96-103.

[3] BARNDORFF-NIELSEN, O. *Parametric statistical models and likelihood*, vol. 50 of *Lecture Notes in Statistics*. Springer Verlag, Berlin, Germany, 1986.

[4] BOSTAN, A., LECERF, G., AND SCHOST, E. Tellegen's principle into practice. In *Proceedings of the 2003 International Symposium on Symbolic and Algebraic Computation* (New York, NY, USA, 2003), ISSAC '03, Association for Computing Machinery, pp. 37-44.

[5] BUTLER, R. *Saddlepoint approximations with applications*. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, UK, 2007.

[6] DANIELS, H. Saddlepoint approximations in statistics. *The annals of mathematical statistics* 25, 4 (December 1954), 631-650.

[7] DANIELS, H. Tail probability approximations. *International Statistical Review* 55, 1 (April 1987), 37-48.

[8] DAS, A. New methods to compute the generalized chi-square distribution. *Journal of Statistical Computations* 95, 3 (May 2025), 1-36.

[9] DAVIES, R. Numerical inversion of a characteristic function. *Biometrika* 60, 2 (August 1973), 415-417. https://doi.org/10.1093/biomet/60.2.415.

[10] DAVIES, R. Algorithm as 155: The distribution of a linear combination of $\xi^2$ random variables. *Journal of the Royal Statistical Society. Series C (Applied Statistics)* 29, 3 (1980), 323-333. https://doi.org/10.2307/2346911.

[11] GAWRONSKI, W., MÜLLER, J., AND REINHARD, M. Reduced cancellation in the evaluation of entire functions and applications to the error function. *SIAM J. Numerical Analysis* 45, 6 (2007), 2564-2576. https://www.jstor.org/stable/40233192.

[12] HIGHAM, N. *Accuracy and stability of numerical algorithms*. SIAM, Philadelphia, PA, USA, 2002.

[13] IMHOF, J. Computing the distribution of quadratic forms in normal variables. *Biometrika* 48, 3/4 (December 1961), 419-426.
[14] JOHNSON, N., AND KOTZ, S. *Continuous univariate distributions, volume 2.*
Wiley Series in Probability and Mathematical Statistics. John Wiley & Sons, New
York, NY, USA, 1970. https://doi.org/10.1080/00224065.1996.11979675.

[15] KUONEN, D. Saddlepoint approximations for distributions of quadratic forms in normal variables. *Biometrika* **86**, 4 (12 1999), 929–935.

[16] LIU, H., TANG, Y., AND ZHANG, H. H. A new chi-square approximation to the distribution of non-negative definie quadratic forms in non-central normal variables. *Computational Statistics and Data Analysis* **53**, 4 (2024), 853–856.

[17] MASSON, M., ARZELIER, D., BRÉHARD, F., JOLDES, M., AND SALVY, B. Fast and reliable computation of the instantaneous orbital collision probability. *Journal of Guidance, Control, and Dynamics* **47**, 8 (August 2024), 1660–1673.

[18] MATHAI, A. M., AND PROVOST, S. B. *Quadratic forms in random variable: Theory and Applications*. Statistics: Texbooks and Monographs. Marcel Dekker, New York, USA, 1992. https://doi.org/10.2307/2290674.

[19] MEZZAROBBA, M. Rounding error analysis of linear recurrences using generating series. *Electronic Transactions on Numerical Analysis* **58** (2023), 196–227.

[20] MOROZ, G. New data structure for univariate polynomial approximation and applications to root isolation, numerical multipoint evaluation, and other problems. In *2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS)* (2022), pp. 1090–1099.

[21] OAKU, T. Algorithms for *d*-modules, integration, and generalized functions with applications to statistics. In *The 50th Anniversary of Gröbner Bases*, vol. 77. Mathematical Society of Japan, 2018, pp. 253–353.

[22] RUBEN, H. Probability content of regions under spherical normal distributions iv: the distribution of homogeneous and non-homogeneous quadratic functions of normal variables. *The Annals of Mathematical Statistics* **33** (1962), 542–570.

[23] RUBEN, H. A new result on the distribution of quadratic forms. *The Annals of Mathematical Statistics* **34**, 3 (1963), 1582–1584.

[24] SERRA, R., ARZELIER, D., JOLDES, M., LASSERRE, J., RONDEPIERRE, A., AND SALVY, B. Fast and accurate computation of orbital collision probability for short-term encounters. *Journal of Guidance Control and Dynamics* **39**, 9 (2016), 1009–1021. https://doi.org/10.2514/1.G001353.

[25] TEMME, N. The uniform asymptotic expansion of a class of integrals related to cumulative distribution functions. *SIAM J. Math. Anal.* **13**, 2 (March 1982), 239–253.

[26] WIDDER, D. *An Introduction to Transform Theory*. Acad. Press, 1971.
A Input parameters for the numerical experiments

<table><caption>Table 2: Input parameters *d*, *p*<sub>i</sub> and *m*<sub>i</sub> for the set of cases from Table 2 in [8]. The *p*<sub>i</sub>s and *m*<sub>i</sub>s are presented in a vector with the multiplicity of each different value for the *p*<sub>i</sub>s. The multiplicity of the different values of *m*<sub>i</sub>s may be deduced from the previous rule.</caption><thead><tr><th rowspan="2">Case #</th><th colspan="2">Input parameters</th></tr><tr><th>p<sub>i</sub></th><th>m<sub>i</sub></th></tr></thead><tbody><tr><td rowspan="3">1</td><td>[0.05(1), 0.0288675(1), 0.0204124(1)]</td><td>[0, 0, 0]</td></tr><tr><td>[0.926013(1), 0.534634(1), 0.378043(1)]</td><td>[0, 0, 0]</td></tr><tr><td>[4.47214(1), 2.58199(1), 1.82574(1)]</td><td>[0, 0, 0]</td></tr><tr><td rowspan="3">2</td><td>[0.141421(2), 0.0816497(2), 0.057735(2)]</td><td>[0, 0, 0]</td></tr><tr><td>[4.47214(2), 2.58199(2), 1.82574(2)]</td><td>[0, 0, 0]</td></tr><tr><td>[23.2379(2), 13.4164(2), 9.48683(2)]</td><td>[0, 0, 0]</td></tr><tr><td rowspan="3">3</td><td>[1.58114(2), 0.912871(4), 0.645497(6)]</td><td>[0, 0, 0]</td></tr><tr><td>[17.6777(2), 10.2062(4), 7.21688(6)]</td><td>[0, 0, 0]</td></tr><tr><td>[65.7267(2), 37.9473(4), 26.8328(6)]</td><td>[0, 0, 0]</td></tr><tr><td rowspan="3">4</td><td>[1.58114(6), 0.912871(4), 0.645497(2)]</td><td>[0, 0, 0]</td></tr><tr><td>[8.21584(6), 4.74342(4), 3.3541(2)]</td><td>[0, 0, 0]</td></tr><tr><td>[35.7771(6), 20.6559(4), 14.6059(2)]</td><td>[0, 0, 0]</td></tr><tr><td rowspan="3">5</td><td>[2.58199(2), 1.69031(6)]</td><td>[0, 0]</td></tr><tr><td>[28.8675(2), 18.8982(6)]</td><td>[0, 0]</td></tr><tr><td>[81.6497(2), 53.4522(6)]</td><td>[0, 0]</td></tr><tr><td rowspan="3">6</td><td>[0.912871(1), 0.597614(1)]</td><td>[0.547723, 2.50998]</td></tr><tr><td>[13.4164(1), 8.7831(1)]</td><td>[0.223607, 1.0247]</td></tr><tr><td>[53.033(1), 34.7183(1)]</td><td>[0.141421, 0.648074]</td></tr><tr><td rowspan="3">7</td><td>[5.03115(2), 3.55756(6), 2.90474(4), 2.05396(10), 1.45237(2)]</td><td>[0, 0, 0, 0]</td></tr><tr><td>[21.9089(2), 15.4919(6), 12.6491(4), 8.94427(10), 6.32456(2)]</td><td>[0, 0, 0, 0]</td></tr><tr><td>[50.7198(2), 35.8643(6), 29.2831(4), 20.7063(10), 14.6416(2)]</td><td>[0, 0, 0, 0]</td></tr><tr><td rowspan="3">9</td><td>[8.4533(3), 5.53399(7)]</td><td>[0, 0]</td></tr><tr><td>[29.2119(3), 19.1237(7)]</td><td>[0, 0]</td></tr><tr><td>[60.5117(3), 39.6142(7)]</td><td>[0, 0]</td></tr><tr><td rowspan="3">11</td><td>[16.4317(8), 9.48683(11), 6.7082(8), 6.21059(7)]</td><td>[0, 0, 0]</td></tr><tr><td>[46.4758(8), 26.8328(11), 18.9737(8), 17.5662(7)]</td><td>[0, 0, 0, 0]</td></tr><tr><td>[100(8), 57.735(11), 40.8248(8), 37.7964(7)]</td><td>[0, 0, 0, 0]</td></tr><tr><td rowspan="3">13</td><td>[4.47214(1), 2.23607(2), 2(1)]</td><td>[0.0894427, 0, 0.25]</td></tr><tr><td>[23.2379(1), 11.619(2), 10.3923(1)]</td><td>[0.0516398, 0, 0.144338]</td></tr><tr><td>[35.7771(1), 17.8885(2), 16(1)]</td><td>[0.0447214, 0, 0.125]</td></tr><tr><td rowspan="3">14</td><td>[0.912871(1), 0.597614(1)]</td><td>[0.547723, 2.50998]</td></tr><tr><td>[13.4164(1), 8.7831(1)]</td><td>[0.223607, 1.0247]</td></tr><tr><td>[53.033(1), 34.7183(1)]</td><td>[0.141421, 0.648074]</td></tr><tr><td rowspan="3">15</td><td>[20(2), 1.41776(1)]</td><td>[0, 0.352668]</td></tr><tr><td>[160(2), 11.3421(1)]</td><td>[0, 0.176334]</td></tr><tr><td>[293.939(2), 20.8368(1)]</td><td>[0, 0.143976]</td></tr><tr><td rowspan="3">16</td><td>[8.4533(3), 5.53399(7)]</td><td>[0, 0]</td></tr><tr><td>[29.2119(3), 19.1237(7)]</td><td>[0, 0]</td></tr><tr><td>[60.5117(3), 39.6142(7)]</td><td>[0, 0]</td></tr></tbody></table>
B Tail Bounds and Explicit Truncation Indices

**Lemma B.1** (Tail bound based on Laplace transform). Fixing $r \in (0, 1/p)$, the remainder of the sum $\sum_{k \ge 0} c_k$ satisfies

$$ \sum_{k \ge N} c_k \le \frac{\hat{f}(r)}{\Gamma(N + \frac{d}{2} + 1) r^N}, $$

for every

$$ N \ge N_0(r) := \left\lfloor \frac{1}{r} - \left(\frac{d}{2} + 1\right) \right\rfloor. $$

*Proof.* Setting $\alpha := \frac{d}{2} + 1$, then by definition $r(N + \alpha) \ge 1$ if $N \ge N_0(r)$. Since $N + \alpha + k \ge N + \alpha$ for all $k \ge 0$, we have

$$ r(N + \alpha + k) \ge r(N + \alpha) \ge 1, $$

and for $n = N + k, \forall k \ge 0$, we have therefore,

$$ \frac{\Gamma(n + \alpha)r^n}{\Gamma(N + \alpha)r^N} = (N + \alpha + k)(N + \alpha + k - 1)\cdots(N + \alpha)r^k \ge 1, \forall n \ge N. $$

Since $c_n \ge 0$, multiplying by $c_n$ and summing over $n \ge N$ gives

$$ \sum_{n \ge N} c_n \le \frac{1}{\Gamma(N + \alpha) r^N} \sum_{n \ge N} c_n \Gamma(n + \alpha) r^n. \qquad (21) $$

Since $|r| < 1/p$, $\sum_{n \ge 0} c_n \Gamma(n + \alpha) r^n$ is convergent and is equal to $\hat{f}(r)$, majorizing the last sum in (21). This completes the proof. $\square$

**Lemma B.2** (Geometric majorant for the coefficients of $\varphi$). By definition $p \ge \alpha_i > 0$ for $i = 2, \dots, d$ and let

$$ \varphi(\lambda) = \sum_{k \ge 0} \varphi_k \lambda^k, \quad |\lambda| < \frac{1}{p}. $$

Define the majorant coefficients

$$ \bar{\varphi}_k := \varphi_0^{k+1}, \quad (k \ge 0), $$

then, for all $k \ge 0$, $\varphi_k \le \bar{\varphi}_k$.

*Proof.* We prove by induction that $\varphi_k \le \bar{\varphi}_k$ for all $k \ge 0$. By the geometric series expansions we have

$$ \varphi_0 = m_1 p + p + \sum_{i=2}^{d} \left( \frac{\alpha_i}{2} + p_i m_i \right), $$
and for every $k \ge 1$,

$$ \varphi_k = p^{k+1} + \sum_{i=2}^{d} \left( \frac{\alpha_i}{2} \alpha_i^k + p_i m_i (k+1) \alpha_i^k \right). $$

The base case is trivial and for the induction step assume $\varphi_k \le \bar{\varphi}_k = \varphi_0^{k+1}$. Multiplying by $\varphi_0 > 0$ gives

$$ \bar{\varphi}_{k+1} = \varphi_0^{k+2} \ge \varphi_0 \varphi_k. $$

Using the explicit expression of $\varphi_k$ and the positivity of all terms, we lower bound $\varphi_0 \varphi_k$ by keeping only the contributions of $\varphi_0$ coming from $p$ and $\sum_{i=2}^d p_i m_i$:

$$ \varphi_0 \varphi_k \ge \left( p + \sum_{i=2}^{d} p_i m_i \right) \left( p^{k+1} + \sum_{i=2}^{d} \left( \frac{\alpha_i}{2} \alpha_i^k + p_i m_i (k+1) \alpha_i^k \right) \right). $$

Expanding the right-hand side and keeping the terms needed below yields

$$ \varphi_0 \varphi_k \ge p^{k+2} + \sum_{i=2}^{d} \left( \frac{\alpha_i}{2} \alpha_i^k p + p_i m_i (k+1) \alpha_i^k p + p^{k+1} p_i m_i \right). $$

Since $p > \alpha_i$, $\frac{\alpha_i}{2} \alpha_i^k p \ge \frac{\alpha_i}{2} \alpha_i^{k+1}$ and $p^{k+1} p_i m_i \ge p_i m_i \alpha_i^k p$, so we get

$$ \bar{\varphi}_{k+1} \ge \varphi_0 \varphi_k \ge \varphi_{k+1}. $$

□

It follows that the series

$$ \bar{\varphi}(\lambda) := \sum_{k \ge 0} \bar{\varphi}_k \lambda^k = \sum_{k \ge 0} \varphi_0^{k+1} \lambda^k = \frac{\varphi_0}{1 - \varphi_0 \lambda}, \quad |\lambda| < \frac{1}{\varphi_0}. $$

is a convergent majorizing series for $\varphi$ i.e., $0 \le \varphi_k \le \bar{\varphi}_k$ for all $k \ge 0$. Equivalently, in the coefficientwise order on power series, $0 \ll \varphi \ll \bar{\varphi}$. Using (9), we define the majorant equation

$$ \bar{f}'(\lambda) = \bar{\varphi}(\lambda) \bar{f}(\lambda), \quad \bar{f}(0) = \hat{f}_0. $$

whose solution is denoted by $\bar{f}$. By the majorant-series comparison results of [19, Lemmas 6.3 and 6.5], the coefficientwise inequality $0 \ll \varphi \ll \bar{\varphi}$ implies

$$ 0 \ll \hat{f} \ll \bar{f}. $$

Finally, the majorant ODE integrates explicitly and yields

$$ \bar{f}(\lambda) = \frac{\hat{f}_0}{1 - \varphi_0 \lambda}. $$

Hence,

$$ 0 \ll \hat{f}(\lambda) \ll \frac{\hat{f}_0}{1 - \varphi_0 \lambda}, $$

on the common domain where both series are defined (in particular for $|\lambda| < \min\{1/p, 1/\varphi_0\}$).
**Lemma B.3** (Tail bound 2). For every $N \ge 0$, the remainder of the sum $\sum_{n \ge 0} c_n$ satisfies

$$ \sum_{n \ge N} c_n \le \frac{\hat{f}(0)e^{\varphi_0} \varphi_0^N}{\Gamma(N + d/2 + 1)}. $$

*Proof.* If follows from the definition of $\hat{f}$ and an exponential tail bound:

$$ \sum_{n \ge N} c_n \le \sum_{n \ge N} \frac{\hat{f}(0) \varphi_0^n}{\Gamma(n + d/2 + 1)} \le \frac{\hat{f}(0) \varphi_0^N}{\Gamma(N + d/2 + 1)} \sum_{n \ge 0} \frac{\varphi_0^n}{n!}. $$

□

Since $c_n \ge 0$, it is sufficient to impose

$$ e^{-p} \sum_{n \ge N+1} c_n \le \delta, \qquad (22) $$

to get a bound $\delta$ on the truncation error. Using the two tail bounds as well as fixing $r := \frac{1}{2p} \in (0, 1/p)$ in the first one, a sufficient condition for (22) is that $N$ satisfies at least one of the following:

(i) $N \ge \lceil 2p - (\frac{d}{2} + 1) \rceil - 1$ and

$$ e^{-p} \frac{\hat{f}\left(\frac{1}{2p}\right)}{\Gamma\left(N + \frac{d}{2} + 2\right) \left(\frac{1}{2p}\right)^{N+1}} \le \delta; $$

(ii)

$$ e^{-p} \frac{\hat{f}(0) e^{\varphi_0} \varphi_0^{N+1}}{\Gamma\left(N + \frac{d}{2} + 2\right)} \le \delta. $$

**Lemma B.4.** Let $A > 0$, $C > 0$, and $\delta \in (0, 1)$. For an integer $N \ge 1$, set

$$ E_N := \frac{C A^N}{N!}. $$

*Define*

$$ N_1 := \lceil 2eA \rceil. $$

*If*

$$ N \ge \max\left\{N_1, \left\lceil \log_2\left(\frac{C}{\delta\sqrt{2\pi N_1}}\right) \right\rceil\right\}, $$

then $E_N \le \delta$.
*Proof.* Using the Stirling inequality

$$N! \geq \sqrt{2\pi} N^{N+1/2} e^{-N},$$

we obtain

$$E_N \leq \frac{C}{\sqrt{2\pi N}} \left(\frac{eA}{N}\right)^N.$$

If $N \geq N_1$, then $\frac{eA}{N} \leq \frac{1}{2}$ and hence

$$E_N \leq \frac{C}{\sqrt{2\pi N}} 2^{-N} \leq \frac{C}{\sqrt{2\pi N_1}} 2^{-N}.$$

If moreover

$$N \geq \log_2\left(\frac{C}{\delta\sqrt{2\pi N_1}}\right),$$

then $\frac{C}{\sqrt{2\pi N_1}} 2^{-N} \leq \delta$, proving $E_N \leq \delta$. $\square$

**Lemma B.5.** Let $A > 0$, $C > 0$, $\delta \in (0, 1)$, and let $\alpha \geq 1$. For an integer $N \geq 0$, set

$$E_N := \frac{C A^{N+1}}{\Gamma(N + 1 + \alpha)}.$$

Define

$$N_1 := \lceil 2eA \rceil.$$

If

$$N \geq \max\{N_1 - 1, \left\lfloor \log_2\left(\frac{C}{\delta\Gamma(\alpha)\sqrt{2\pi N_1}}\right) \right\rfloor - 1\},$$

then $E_N \leq \delta$.

*Proof.* Since $\alpha \geq 1$, we have $\Gamma(N + 1 + \alpha) \geq \Gamma(\alpha)(N + 1)!$ and therefore

$$E_N \leq \frac{C}{\Gamma(\alpha)} \frac{A^{N+1}}{(N+1)!}.$$

Apply Lemma B.4 with

$$A \mapsto A, \quad C \mapsto \frac{C}{\Gamma(\alpha)}, \quad \delta \mapsto \delta, \quad \text{and} \quad N \mapsto N + 1.$$

With $N_1 = \lceil 2eA \rceil$ as in Lemma B.4, the conditions become

$$N + 1 \geq \max\{N_1, \left\lfloor \log_2\left(\frac{C}{\delta\Gamma(\alpha)\sqrt{2\pi N_1}}\right) \right\rfloor\}.$$

Under this assumption, Lemma B.4 yields

$$\frac{C}{\Gamma(\alpha)} \frac{A^{N+1}}{(N+1)!} \leq \delta,$$

hence $E_N \leq \delta$. $\square$
**Lemma B.6 (Explicit truncation index).** Fix $\delta \in (0, 1)$ and set $\alpha := \frac{d}{2} + 1$. Define

$$N_{(ii)}(\delta) := \max\left\{\lceil 2e \varphi_0 \rceil - 1, \left\lceil \log_2 \left( \frac{e^{-p} \hat{f}(0)e^{\varphi_0}}{\delta \Gamma(\alpha) \sqrt{2\pi \lceil 2e \varphi_0 \rceil}} \right) \right\rceil - 1 \right\},$$

and

$$N_{(i)}(\delta) := \max\left\{\lceil 2p - \alpha \rceil - 1, \lceil 4e p \rceil - 1, \left\lceil \log_2 \left( \frac{e^{-p} \hat{f}\left(\frac{1}{2p}\right)}{\delta \Gamma(\alpha) \sqrt{2\pi \lceil 4e p \rceil}} \right) \right\rceil - 1 \right\}.$$

Finally, set

$$N_{\star}(\delta) := \min\{N_{(i)}(\delta), N_{(ii)}(\delta)\}.$$

Then, for any integer $N \ge N_{\star}(\delta)$, the truncation error satisfies

$$\left| e^{-p} \sum_{n \ge 0} c_n - e^{-p} \sum_{n=0}^{N} c_n \right| \le \delta.$$

*Proof.* Since $c_n \ge 0$, the truncation error equals $e^{-p} \sum_{n \ge N+1} c_n$. Fixing $r = \frac{1}{2p}$, the two tail bounds give the two sufficient conditions (i) and (ii) stated above for $e^{-p} \sum_{n \ge N+1} c_n \le \delta$.

For (ii), apply Lemma B.5 with $\alpha = \frac{d}{2} + 1$, $A = \varphi_0$, and $C = e^{-p} \hat{f}(0)e^{\varphi_0}$, which yields that (ii) holds for all $N \ge N_{(ii)}(\delta)$.

For (i), apply Lemma B.5 with $\alpha = \frac{d}{2} + 1$, $A = 2p$, and $C = e^{-p} \hat{f}\left(\frac{1}{2p}\right)$, and add the validity constraint $N \ge \lceil 2p - \alpha \rceil - 1$. This yields that (i) holds for all $N \ge N_{(i)}(\delta)$.

Therefore, for any $N \ge \min\{N_{(i)}(\delta), N_{(ii)}(\delta)\} = N_{\star}(\delta)$, at least one of (i) or (ii) holds, hence $e^{-p} \sum_{n \ge N+1} c_n \le \delta$, proving the claim. $\square$

## C Evaluation error for CDFREC recurrence

We use the same (Higham) floating-point model, with precision $t$, unit roundoff $u = 2^{-t}$, rounding to nearest, and no overflow/underflow, assuming multiplications by powers of 2 are exact and appearing real constants and elementary functions are faithfully rounded (relative error $< 2u$).

**Exact and computed quantities.** Let $(c_n)_{n \ge 0}$ be the exact coefficients defined by (11) with $c_0 = \hat{f}(0)/\Gamma(\frac{d}{2} + 1)$. Let $(\tilde{c}_n)_{n \ge 0}$ be the floating-point values produced by Algorithm 1, including rounding in the parameters $\alpha_i$ and in the derived coefficients $P_j, Q_j$, as well as in the loop.
**Local and global rounding errors.** Define, for $n \ge 1$, the exact functional corresponding to the recurrence step:

$$F_n(x_{n-1}, \dots, x_{n-2d}) := \frac{1}{n} \sum_{i=1}^{\min(2d,n)} (-1)^{i-1} \frac{P_{i-1} + Q_i(n-i)}{\mathrm{Fac}(n,i)} x_{n-i},$$

where

$$\mathrm{Fac}(n, i) := \prod_{k=1}^{i} \left( \frac{d}{2} + n + 1 - k \right) = \prod_{j=0}^{i-1} \left( \frac{d}{2} + n - j \right), \quad \mathrm{Fac}(n, 1) = \frac{d}{2} + n.$$

The global (absolute) error is

$$\delta_n := c_n - \tilde{c}_n, \quad n \ge 0,$$

and the local (absolute) error injected at iteration $n$ is

$$\varepsilon_n := \tilde{c}_n - F_n(\tilde{c}_{n-1}, \dots, \tilde{c}_{n-2d}), \quad n \ge 1.$$

Since $c_n = F_n(c_{n-1}, \dots, c_{n-2d})$ and $\tilde{c}_n = F_n(\tilde{c}_{n-1}, \dots, \tilde{c}_{n-2d}) + \varepsilon_n$, the global errors satisfy the inhomogeneous recurrence

$$\delta_n = F_n(\delta_{n-1}, \dots, \delta_{n-2d}) - \varepsilon_n, \quad n \ge 1, \qquad (23)$$

with $\delta_0 = c_0 - \tilde{c}_0$ capturing the rounding in $c_0$.

**Lifting to generating series.** Recall

$$\hat{f}(\lambda) = \sum_{n \ge 0} c_n \Gamma\left(n + \frac{d}{2} + 1\right) \lambda^n.$$

Set $\Gamma_n := \Gamma(n + \frac{d}{2} + 1)$ and define similarly the generating series of global and local errors:

$$\hat{\delta}(\lambda) := \sum_{n \ge 0} \delta_n \Gamma_n \lambda^n, \quad \hat{\varepsilon}(\lambda) := \sum_{n \ge 1} \varepsilon_n \Gamma_n \lambda^n.$$

**Lemma C.1.** The series $\hat{\delta}$ satisfies the first-order inhomogeneous differential equation

$$Q(\lambda) \hat{\delta}'(\lambda) - P(\lambda) \hat{\delta}(\lambda) = -\hat{\varepsilon}'(\lambda), \qquad (24)$$

with initial value $\hat{\delta}(0) = \Gamma(\frac{d}{2} + 1)\delta_0$.

*Proof.* Note that $\Gamma_n = \mathrm{Fac}(n, i)\Gamma_{n-i}$ and we have the formal power series identities in the ring $\mathbb{R}[[\lambda]]$:

$$\lambda \hat{\delta}'(\lambda) = \sum_{n \ge 1} n \delta_n \Gamma_n \lambda^n, \quad \lambda \hat{\varepsilon}'(\lambda) = \sum_{n \ge 1} n \varepsilon_n \Gamma_n \lambda^n, \qquad (25)$$

$$(\sum_{i \ge 1} a_i \lambda^i) \hat{\delta}(\lambda) = \sum_{n \ge 1} \left( \sum_{i=1}^{\min(n,m)} a_i \delta_{n-i} \Gamma_{n-i} \right) \lambda^n, \qquad (26)$$

$$(\sum_{i \ge 1} a_i \lambda^i) \lambda \hat{\delta}'(\lambda) = \sum_{n \ge 1} \left( \sum_{i=1}^{\min(n,m)} a_i (n-i) \delta_{n-i} \Gamma_{n-i} \right) \lambda^n, \qquad (27)$$
for any polynomial $\sum_{i=1}^{m} a_i \lambda^i$ (Cauchy product). Moreover, with

$$Q(\lambda) = \sum_{j=0}^{2d-1} (-1)^j Q_j \lambda^j, \quad P(\lambda) = \sum_{j=0}^{2d-1} (-1)^j P_j \lambda^j,$$

we have the index-shift expansions

$$\sum_{i \ge 1} (-1)^i P_{i-1} \lambda^i = -\lambda P(\lambda), \quad \sum_{i \ge 1} (-1)^i Q_i \lambda^i = Q(\lambda) - 1, \qquad (28)$$

Multiply (23) by $n\Gamma_n\lambda^n$, and sum over $n \ge 1$. Using the identities above one obtains

$$Q(\lambda) \lambda \hat{\delta}'(\lambda) - P(\lambda) \lambda \hat{\delta}(\lambda) = -\lambda \hat{\varepsilon}'(\lambda),$$

and conclude by formally dividing by $\lambda$. $\square$

**Coefficientwise positivity and majorants.** Since $p_1 > 0$ and $\alpha_i > 0$, the expansions of $(1 - p_1\lambda)^{-1}$, $(1 - \alpha_i\lambda)^{-1}$ and $(1 - \alpha_i\lambda)^{-2}$ have nonnegative coefficients. In particular

$$\varphi(\lambda) \gg 0, \quad Q(\lambda)^{-1} = \frac{1}{(1 - p_1\lambda) \prod_{i=2}^{d} (1 - \alpha_i\lambda)^2} \gg 0.$$

For the alternating polynomials $P(\lambda) = \sum_{j=0}^{2d-1} (-1)^j P_j \lambda^j$ and $Q(\lambda) = \sum_{j=0}^{2d-1} (-1)^j Q_j \lambda^j$ with $P_j, Q_j \ge 0$, we introduce the positive-coefficient majorants

$$P^+(\lambda) := \sum_{j=0}^{2d-1} P_j \lambda^j, \quad Q^+(\lambda) := \sum_{j=0}^{2d-2} Q_{j+1} \lambda^j.$$

**Relative perturbation model for $P_j$ and $Q_j$.** Assume the precomputed coefficients satisfy a componentwise relative error model: there exist $\varepsilon_j^P, \varepsilon_j^Q$ such that

$$\tilde{P}_j = P_j(1 + \varepsilon_j^P), \quad \tilde{Q}_j = Q_j(1 + \varepsilon_j^Q), \quad |\varepsilon_j^P| \le \gamma_{\text{coef}}, \quad |\varepsilon_j^Q| \le \gamma_{\text{coef}},$$

for all relevant indices $j$.

**Local error bound.** Since the inner loop has length $m = \min(2d, n) \le 2d$, the number of floating-point operations performed in one iteration is uniformly bounded by a constant $K$ independent of $n$. By Higham's model, the rounding contribution of these operations can be written as a single multiplicative factor $(1 + \theta)$ with $|\theta| \le \gamma_K$. Combining this with the coefficient perturbations and using $(1 + \gamma_a)(1 + \gamma_b) - 1 \le \gamma_{a+b}$, we absorb both effects into the single parameter

$$\gamma := (1 + \gamma_K)(1 + \gamma_{\text{coef}}) - 1 \le \gamma_{K+K_{\text{coef}}}, \quad \text{assuming } u(K + K_{\text{coef}}) < 1,$$

and $\gamma_{\text{coef}} = \gamma_{K_{\text{coef}}}$ for some $K_{\text{coef}}$ independent of $n$.
We now bound K directly from Algorithm 1. For each inner-loop index i, the update

$$ \text{Fac} \leftarrow (d/2 + n + 1 - i) \text{Fac} $$

uses one floating-point multiplication. The term added to $c_n$ is evaluated as

$$ \frac{\tilde{P}_{i-1} + \tilde{Q}_i(n-i)}{\text{Fac}} \tilde{c}_{n-i}, $$

which can be computed using one subtraction (to form $n-i$), one multiplication (for $\tilde{Q}_i(n-i)$), one addition (add $\tilde{P}_{i-1}$), one division (by Fac), and one multiplication (by $\tilde{c}_{n-i}$), followed by one addition to accumulate into $c_n$. Hence there are 7 floating-point operations per inner-loop iteration. Since $m \le 2d$ and the final statement $c_n \leftarrow c_n/n$ costs one floating-point division, we obtain the uniform bound

$$ K \le 7m + 1 \le 14d + 1, \quad \text{and thus} \quad \gamma_K \le \gamma_{14d+1}. $$

Hence, letting

$$ A_{n,i} := \frac{P_{i-1} + Q_i(n-i)}{n \text{Fac}(n,i)}, \quad (1 \le i \le \min(2d,n)), $$

and $\gamma := \gamma_{14d+1+K_{\text{coef}}}$ we have

$$ |\varepsilon_n| \le \gamma \sum_{i=1}^{\min(2d,n)} A_{n,i} |\tilde{c}_{n-i}|, \quad n \ge 1. $$

Using $|\tilde{c}_k| \le c_k + |\delta_k|$ (note that we use the fact that the exact $c_k$ are non-negative) gives

$$ |\varepsilon_n| \le \gamma \sum_{i=1}^{\min(2d,n)} A_{n,i} (c_{n-i} + |\delta_{n-i}|). $$

Multiplying the previous inequality by $n\Gamma_n\lambda^n$ and summing over $n \ge 1$ yields

$$ \hat{\varepsilon}'(\lambda) \ll \gamma \left( (Q^+(\lambda)\varphi(\lambda) + P^+(\lambda)) \hat{f}(\lambda) + Q^+(\lambda) |\hat{\delta}|'(\lambda) + P^+(\lambda) |\hat{\delta}|(\lambda) \right). \quad (29) $$

**Majorant differential inequality (linearized form).** From equation (24) and since $Q(\lambda)^{-1} \gg 0$ as well $\varphi(\lambda) = P(\lambda)/Q(\lambda) \gg 0$, we obtain

$$ |\hat{\delta}|'(\lambda) \ll \varphi(\lambda) |\hat{\delta}|(\lambda) + Q(\lambda)^{-1} |\hat{\varepsilon}'(\lambda)|. $$

Combining with (29) yields

$$ |\hat{\delta}|'(\lambda) \ll \varphi(\lambda) |\hat{\delta}|(\lambda) + \gamma \left( \frac{Q^+(\lambda)}{Q(\lambda)} |\hat{\delta}|'(\lambda) + \frac{P^+(\lambda)}{Q(\lambda)} |\hat{\delta}|(\lambda) + \frac{Q^+(\lambda)\varphi(\lambda) + P^+(\lambda)}{Q(\lambda)} \hat{f}(\lambda) \right). $$
In particular, neglecting the $O(u^2)$ terms (i.e. products of $\gamma = O(u)$ with $|\hat{\delta}|$ and $|\hat{\delta}|'$ under the bootstrap assumption $|\hat{\delta}|, |\hat{\delta}|' = O(u)$), we obtain the linearized majorant inequality

$$|\hat{\delta}|'(\lambda) \ll \varphi(\lambda) |\hat{\delta}|(\lambda) + \gamma \hat{\psi}(\lambda) \hat{f}(\lambda),$$

where

$$\hat{\psi}(\lambda) := \frac{Q^{+}(\lambda)\varphi(\lambda) + P^{+}(\lambda)}{Q(\lambda)} \gg 0.$$

**A simpler majorant for $Q^{-1}$.** Since $p_1 = \max\{p_1, \alpha_2, \dots, \alpha_d\}$ and $\alpha_i > 0$, then for each $i \ge 2$ we have the coefficientwise majorants

$$(1 - \alpha_i \lambda)^{-1} \ll (1 - p_1 \lambda)^{-1}, \quad (1 - \alpha_i \lambda)^{-2} \ll (1 - p_1 \lambda)^{-2},$$

since $\alpha_i \le p_1$ implies $\alpha_i^n \le p_1^n$ for all $n \ge 0$. Multiplying these inequalities yields

$$0 \ll Q(\lambda)^{-1} = \frac{1}{(1 - p_1 \lambda) \prod_{i=2}^{d} (1 - \alpha_i \lambda)^2} \ll (1 - p_1 \lambda)^{-(2d-1)}.$$

Similarly, we obtain the simpler coefficientwise majorant

$$\varphi(\lambda) \ll \varphi^{\#}(\lambda) := m_1 p_1 + \frac{p_1 + \frac{1}{2} \sum_{i=2}^{d} \alpha_i}{1 - p_1 \lambda} + \frac{\sum_{i=2}^{d} p_i m_i}{(1 - p_1 \lambda)^2}.$$

Therefore

$$\hat{\psi}(\lambda) \ll \hat{\psi}^{\#}(\lambda) := \frac{Q^{+}(\lambda)\varphi^{\#}(\lambda) + P^{+}(\lambda)}{(1 - p_1 \lambda)^{2d-1}} \quad (30)$$

Replacing $\hat{\psi}$ by $\hat{\psi}^{\#}$ yields

$$|\hat{\delta}|'(\lambda) \ll \varphi(\lambda) |\hat{\delta}|(\lambda) + \gamma \hat{\psi}^{\#}(\lambda) \hat{f}(\lambda).$$

Since all the series on the right-hand side of the linearized majorant inequality have nonnegative coefficients, Lemma 6.5 in [19] implies that

$$|\hat{\delta}|(\lambda) \ll \hat{\Delta}(\lambda),$$

where $\hat{\Delta}$ is the solution of the inhomogeneous first-order LDE

$$\hat{\Delta}'(\lambda) = \varphi(\lambda)\hat{\Delta}(\lambda) + \gamma \hat{\psi}^{\#}(\lambda) \hat{f}(\lambda), \quad \hat{\Delta}(0) = e_0 \hat{f}(0) \ge |\hat{\delta}(0)|. \quad (31)$$

A technical analysis of the solution of this LDE allows for the following global error bound.

**Proposition C.2.** Under the simplified model above, the total rounding error accumulated while computing $f(1)$ satisfies

$$\sum_{n=0}^{+\infty} |\tilde{c}_n - c_n| = |\delta|(1) \le (e_0 + \gamma C^{\#}) f(1),$$
where $e_0$ is defined in Prop. 2.3, $\gamma$ is the global rounding parameter introduced in the local error bound, and

$$C^{\#} := \int_{0}^{1} W^{\#}(\tau) d\tau, \qquad (32)$$

with $W^{\#}$ an explicit degree $2d - 1$ polynomial with nonnegative coefficients given in Lemma C.3 below.

*Proof.* Using $\hat{f}$ as a solution of the homogeneous part of (31) gives the explicit form

$$\hat{\Delta}(\lambda) = (e_0 + \gamma \hat{\Psi}(\lambda)) \hat{f}(\lambda), \quad \hat{\Psi}(\lambda) := \int_{0}^{\lambda} \hat{\psi}^{\#}(\sigma) d\sigma. \qquad (33)$$

To return to the $\xi$-plane, we use the modified Laplace convention

$$\hat{f}(\lambda) = \frac{\mathcal{L}_f(\lambda^{-1})}{\lambda^{d/2+1}}, \quad \hat{\Delta}(\lambda) = \frac{\mathcal{L}_{\Delta}(\lambda^{-1})}{\lambda^{d/2+1}},$$

Equivalently, setting $\lambda = z^{-1}$ and using (33), we have

$$\mathcal{L}_f(z) = z^{-(d/2+1)} \hat{f}(z^{-1}), \quad \mathcal{L}_{\Delta}(z) = z^{-(d/2+1)} \hat{\Delta}(z^{-1}), \\ \mathcal{L}_{\Delta}(z) = (e_0 + \gamma \hat{\Psi}(z^{-1})) \mathcal{L}_f(z).$$

Associating to $\hat{\Psi}$ the standard Laplace transform $\mathcal{L}_{\Psi}$ defined by $\mathcal{L}_{\Psi}(z) := \hat{\Psi}(z^{-1})$, the standard product-convolution theorem (see [26, Chap. 5, §8]) yields

$$\Delta(\xi) = e_0 f(\xi) + \gamma (\Psi * f)(\xi), \quad (\Psi * f)(\xi) := \int_{0}^{\xi} \Psi(\tau) f(\xi - \tau) d\tau.$$

It remains to bound $(\Psi * f)(1)$.

One has from Widder [26, Chap. 5, Lem. 5.2]

$$\Psi(\xi) = \frac{1}{\xi} \mathcal{L}^{-1} \left\{ \frac{\hat{\psi}^{\#}(\lambda^{-1})}{\lambda^2} \right\} (\xi). \qquad (34)$$

We prove in Lemma C.3 below, that $\Psi(\xi)$ admits a coefficientwise majorant of the form

$$\Psi(\xi) \ll W^{\#}(\xi) e^{p_1 \xi}, \quad (0 \le \xi \le 1),$$

where $W^{\#}$ is an explicit polynomial with nonnegative coefficients obtained by termwise inversion of the elementary fractions. Therefore,

$$(\Psi * f)(1) \le \int_{0}^{1} W^{\#}(\tau) e^{p_1 \tau} g(1 - \tau) d\tau \le e^{p_1} g(1) \int_{0}^{1} W^{\#}(\tau) d\tau,$$

(where we used that $g$ is nondecreasing on $[0, 1]$). Hence

$$|\delta|(1) \le \Delta(1) = e_0 f(1) + \gamma(\Psi * f)(1) \le (e_0 + \gamma C^{\#}) f(1),$$

which concludes the proof. $\square$
**Lemma C.3.** For all $\xi \ge 0$, the series $\Psi(\xi)$ such that $\hat{\Psi}(\lambda) = \mathcal{L}_{\Psi}(\lambda^{-1})$ satisfies

$$ \Psi(\xi) \ll W^{\sharp}(\xi) e^{p_1 \xi}, $$

where $W^{\sharp}$ is a degree-2d-1 polynomial with positive coefficients

$$ W^{\sharp}(\xi) = |W_0| + |W_1| p_1 + \sum_{r=2}^{2d+1} \frac{|W_r|}{(r-1)!} \xi^{r-2}, $$

with $W_0 = -\frac{P_{2d-1}}{p_1^{2d-1}}$, $W_1 = C_{2d-2}$, $W_2 = B_{2d-2} + C_{2d-3}$ and more generally, for $3 \le r \le 2d-1$,

$$ W_r = A_{2d+1-r} + B_{2d-r} + C_{2d-1-r}, $$

and

$$ W_{2d} = A_1 + B_0, \quad W_{2d+1} = A_0. $$

where

$$ A(\lambda) = \sum_{k=0}^{2d-1} A_k(\lambda - p_1)^k, \quad B(\lambda) = \sum_{k=0}^{2d-2} B_k(\lambda - p_1)^k, \quad C(\lambda) = \sum_{k \ge 0} C_k(\lambda - p_1)^k, $$

and

$$ \left\{ \begin{aligned} A(\lambda) &= \left(\sum_{i=2}^{d} p_i m_i\right) \lambda \operatorname{rev}_{2d-2}(Q^{+})(\lambda), \\ B(\lambda) &= \left(p_1 + \frac{1}{2} \sum_{i=2}^{d} \alpha_i\right) \operatorname{rev}_{2d-2}(Q^{+})(\lambda), \\ C(\lambda) &= \frac{\operatorname{rev}_{2d-1}(P^{+})(\lambda)}{\lambda^2} + \frac{m_1 p_1 \operatorname{rev}_{2d-2}(Q^{+})(\lambda)}{\lambda}, \end{aligned} \right. \qquad (35) $$

with the convention that for a degree-m polynomial $R(\lambda) = \sum_{j=0}^{m} r_j \lambda^j$ we define its reversal (or reciprocal) by $\operatorname{rev}_m R(\lambda) := \lambda^m R(\lambda^{-1})$.

*Proof.* Expanding the compact form (30), we have

$$ \hat{\psi}^{\sharp}(\lambda) = \frac{P^{+}(\lambda) + m_1 p_1 Q^{+}(\lambda)}{(1 - p_1 \lambda)^{2d-1}} + \frac{\left(p_1 + \frac{1}{2} \sum_{i=2}^{d} \alpha_i\right) Q^{+}(\lambda)}{(1 - p_1 \lambda)^{2d}} + \frac{\left(\sum_{i=2}^{d} p_i m_i\right) Q^{+}(\lambda)}{(1 - p_1 \lambda)^{2d+1}}, $$

and setting $F^{\sharp}(\lambda) := \frac{\hat{\psi}^{\sharp}(\lambda^{-1})}{\lambda^2}$, gives

$$ F^{\sharp}(\lambda) = \frac{\operatorname{rev}_{2d-1}(P^{+})(\lambda)}{\lambda^2(\lambda - p_1)^{2d-1}} + \frac{m_1 p_1 \operatorname{rev}_{2d-2}(Q^{+})(\lambda)}{\lambda(\lambda - p_1)^{2d-1}} \\ + \left(p_1 + \frac{1}{2} \sum_{i=2}^{d} \alpha_i\right) \frac{\operatorname{rev}_{2d-2}(Q^{+})(\lambda)}{(\lambda - p_1)^{2d}} + \left(\sum_{i=2}^{d} p_i m_i\right) \frac{\lambda \operatorname{rev}_{2d-2}(Q^{+})(\lambda)}{(\lambda - p_1)^{2d+1}}, $$
In particular, $F^{\sharp}(\lambda)$ is a rational function whose poles are contained in $\{0, p_1\}$, with pole order at most 2 at $\lambda = 0$ and at most $2d+1$ at $\lambda = p_1$. Consequently, there exist coefficients $W_{-1}, W_0, (W_r)_{1\le r\le 2d+1}$ such that

$$F^{\sharp}(\lambda) = \frac{W_{-1}}{\lambda} + \frac{W_0}{\lambda^2} + \sum_{r=1}^{2d+1} \frac{W_r}{(\lambda - p_1)^r} \quad (36)$$

Moreover, since $F^{\sharp}(\lambda) = O(\lambda^{-2})$ as $\lambda \to \infty$, the coefficient of $1/\lambda$ in the Laurent expansion at infinity vanishes. Expanding (36) at infinity shows that this coefficient is $W_{-1} + W_1$, hence

$$W_{-1} + W_1 = 0.$$

In particular, (36) can be inverted termwise by the usual Laplace rules:

$$\mathcal{L}^{-1}\{F^{\sharp}\}(\xi) = W_{-1} + W_0\xi + \sum_{r=1}^{2d+1} W_r \frac{\xi^{r-1}}{(r-1)!} e^{p_1\xi}.$$

Using $W_{-1} + W_1 = 0$, the constant term cancels in the inverse transform, and we may rewrite it as

$$\mathcal{L}^{-1}\{F^{\sharp}\}(\xi) = W_0\xi + W_1(e^{p_1\xi} - 1) + \sum_{r=2}^{2d+1} W_r \frac{\xi^{r-1}}{(r-1)!} e^{p_1\xi}.$$

Therefore, for $\xi > 0$,

$$|\Psi(\xi)| \le \left| \frac{\mathcal{L}^{-1}\{F^{\sharp}\}(\xi)}{\xi} \right| \le |W_0| + |W_1| \frac{e^{p_1\xi} - 1}{\xi} + \sum_{r=2}^{2d+1} |W_r| \frac{\xi^{r-2}}{(r-1)!} e^{p_1\xi}.$$

Using $\frac{e^{p_1\xi} - 1}{\xi} \le p_1 e^{p_1\xi}$ and $e^{p_1\xi} \ge 1$ yields

$$\Psi(\xi) \ll W^{\sharp}(\xi) e^{p_1\xi}, \quad W^{\sharp}(\xi) := |W_0| + |W_1| p_1 + \sum_{r=2}^{2d+1} \frac{|W_r|}{(r-1)!} \xi^{r-2}.$$

Finally, the coefficients can be identified explicitly. First,

$$W_0 = \lim_{\lambda \to 0} \lambda^2 F^{\sharp}(\lambda) = \frac{\mathrm{rev}_{2d-1}(P^+)(0)}{(-p_1)^{2d-1}} = -\frac{P_{2d-1}}{p_1^{2d-1}}.$$

Next, grouping the terms at the pole $\lambda = p_1$ by pole order,

$$F^{\sharp}(\lambda) = \frac{A(\lambda)}{(\lambda - p_1)^{2d+1}} + \frac{B(\lambda)}{(\lambda - p_1)^{2d}} + \frac{C(\lambda)}{(\lambda - p_1)^{2d-1}},$$

with $A, B, C$ defined in (35) and analytic at $\lambda = p_1$, a Taylor expansion at $\lambda = p_1$ gives

$$F^{\sharp}(\lambda) = \sum_{r=1}^{2d+1} \frac{W_r}{(\lambda - p_1)^r} + \text{(terms analytic at } \lambda = p_1),$$
and hence, writing $A(\lambda) = \sum_{k=0}^{2d-1} A_k(\lambda - p_1)^k$, $B(\lambda) = \sum_{k=0}^{2d-2} B_k(\lambda - p_1)^k$, $C(\lambda) = \sum_{k \ge 0} C_k(\lambda - p_1)^k$, yields, $W_1 = C_{2d-2}$, $W_2 = B_{2d-2} + C_{2d-3}$ and more generally, for $3 \le r \le 2d-1$,

$$W_r = A_{2d+1-r} + B_{2d-r} + C_{2d-1-r},$$

while for the remaining boundary cases one has

$$W_{2d} = A_1 + B_0, \quad W_{2d+1} = A_0.$$

□

**Remark C.4.** *From Lemma C.3, and equation (32), one has*

$$C^{\sharp} = |W_0| + |W_1| p_1 + \sum_{k=1}^{2d} \frac{|W_{k+1}|}{k k!}.$$

Assume that the coefficients $(P_j)_{0 \le j \le 2d-1}$, $(Q_j)_{0 \le j \le 2d-1}$ and the scalars $m_1, \sum_{i=2}^d p_i m_i$, $\sum_{i=2}^d \alpha_i$ are $O(1)$ with respect to $p_1$. Since

$$\mathrm{rev}_{2d-2}(Q^{+})(p_1) = \sum_{j=0}^{2d-2} Q_{j+1} p_1^{2d-2-j} = O(p_1^{2d-2}),$$

we also have

$$(\mathrm{rev}_{2d-2}(Q^{+}))^{(k)}(p_1) = O(p_1^{2d-2-k}), \quad (k \ge 0),$$

and similarly

$$\mathrm{rev}_{2d-1}(P^{+})(p_1) = O(p_1^{2d-1}).$$

So we obtain the orders $W_{2d+1} = A_0 = O(p_1^{2d-1})$, $W_{2d} = A_1 + B_0 = O(p_1^{2d-1})$, $W_{2d-1} = A_2 + B_1 + C_0 = O(p_1^{2d-2})$. More generally, for $r \ge 2$ one has

$$W_r = A_{2d+1-r} + B_{2d-r} + C_{2d-1-r},$$

and the contributions satisfy, as $p_1 \to \infty$,

$$A_{2d+1-r} = O(p_1^{r-2}), \quad B_{2d-r} = O(p_1^{r-1}), \quad C_{2d-1-r} = O(p_1^{r-1}),$$

so that

$$W_r = O(p_1^{r-1}), \quad (2 \le r \le 2d),$$

while the leading pole coefficient is $W_1 = C_{2d-2} = O(1)$, with an additional $O(p_1^{-1})$ contribution from the $P^+$ term. So

$$C^{\sharp} \sim \sum_{k=1}^{2d-1} \frac{p_1^k}{k!}.$$