Error propagation

A ball function evaluates its result at the midpoint of the input and adds to the radius a bound for the propagated error: how much the exact function can differ between the midpoint and any other point of the input ball. For a function $f$ and a ball $[x \pm r]$ , the task is to bound

|f(x + a) - f(x)| \quad \text{for all } |a| \le r,

in terms of $|x|$ and $r$ . Most bounds below increase with $|a|$ , so they are stated in terms of $|a|$ and applied with $|a| = r$ . The rounding error of the midpoint computation is added separately (see Representation of radii and midpoints). All bounds are evaluated with mag_t arithmetic, rounding upward, from upper bounds for the quantities that occur in numerators and lower bounds for those in denominators.

Arithmetic#

Addition is immediate: $|(x + a) + (y + b) - (x + y)| \le |a| + |b|$ , so arb_add adds the radii.

Theorem 1 (Product). For $x, y, a, b \in \mathbb{C}$ ,

|(x + a)(y + b) - xy| \le |x||b| + |y||a| + |a||b|.

Proof. Expand $(x + a)(y + b) - xy = xb + ya + ab$ and apply the triangle inequality.

The bound is attained for real numbers when all of $x, y, a, b$ have the same sign, so for real balls it is optimal up to rounding. arb_mul computes the radius $|m_x| r_y + |m_y| r_x + r_x r_y$ with mag_fast_mul and mag_fast_addmul, using upper bounds of $|m_x|$ and $|m_y|$ obtained from the top 30 bits of the midpoints.

Theorem 2 (Quotient). For $x, y, a, b \in \mathbb{C}$ with $|b| < |y|$ ,

\left| \frac{x}{y} - \frac{x + a}{y + b} \right| \le \frac{|x||b| + |y||a|}{|y| (|y| - |b|)}.

Proof. We have $\dfrac{x}{y} - \dfrac{x + a}{y + b} = \dfrac{x(y + b) - y(x + a)}{y(y + b)} = \dfrac{xb - ya}{y(y + b)}$ . The triangle inequality bounds the numerator by $|x||b| + |y||a|$ , and the reverse triangle inequality gives $|y + b| \ge |y| - |b| > 0$ .

When $r_y$ is much smaller than $|m_y|$ , arb_div uses the simplified form

\frac{|x| r_y + |y| r_x}{|y|(|y| - r_y)} \le (1 + \varepsilon) \frac{|x/y| r_y + r_x}{|y|}, \quad \frac{|y|}{|y| - r_y} \le 1 + \varepsilon,

with the computed quotient in place of $x/y$ and a fixed factor $1 + 2^{-16}$ for $1 + \varepsilon$ , which is valid when the relative accuracy of $y$ is more than 20 bits. For less accurate input it evaluates the quotient at the endpoints instead, which gives a tighter enclosure.

Roots#

Theorem 3 (Square root). For real $x$ , $a$ with $0 \le |a| \le x$ ,

|\sqrt{x + a} - \sqrt{x}| \le \sqrt{x} \left(1 - \sqrt{1 - \frac{|a|}{x}}\right) \le \frac{\sqrt{x}}{2} \left( \frac{|a|}{x} + \frac{|a|^2}{x^2} \right),

where the first inequality is an equality for $a \le 0$ (and the limit value is 0 when $x = a = 0$ ).

Proof. Since the square root is concave, a perturbation of size $|a|$ changes it most when it is negative, so the error is at most $\sqrt{x} - \sqrt{x - |a|} = \sqrt{x}(1 - \sqrt{1 - t})$ with $t = |a|/x \in [0, 1]$ . For the second inequality we need $\sqrt{1 - t} \ge 1 - (t + t^2)/2$ . The right side is nonnegative on $[0, 1]$ , and squaring turns the inequality into $1 - t \ge 1 - t - t^2 + (t + t^2)^2/4$ , that is, $t^2 (t + 3)(t - 1) / 4 \le 0$ , which holds for $0 \le t \le 1$ .

Theorem 4 (Reciprocal square root). For real $x, a$ with $0 \le |a| < x$ ,

\left| \frac{1}{\sqrt{x + a}} - \frac{1}{\sqrt{x}} \right| \le \frac{|a|}{2 (x - |a|)^{3/2}}.

Proof. By the mean value theorem the difference is $|a| \cdot |f'(\xi)|$ for some $\xi$ between $x$ and $x + a$ , where $f(t) = t^{-1/2}$ and $|f'(t)| = t^{-3/2}/2$ is decreasing, so $|f'(\xi)| \le (x - |a|)^{-3/2}/2$ .

Theorem 5 ( $k$ -th root). For $k > 1$ and real $x, a$ with $0 \le |a| \le x$ ,

|(x + a)^{1/k} - x^{1/k}| \le x^{1/k} \min\left(1, \frac{1}{k} \log\left(1 + \frac{|a|}{x - |a|}\right)\right).

Proof. By concavity the error is largest for $a = -r$ , $r = |a|$ , and

x^{1/k} - (x - r)^{1/k} = x^{1/k} \left[1 - \exp\left(\tfrac{1}{k} \log(1 - r/x)\right)\right] \le x^{1/k} \min\left(1, -\tfrac{1}{k} \log(1 - r/x)\right),

using $1 - e^{-u} \le \min(1, u)$ for $u \ge 0$ . Finally $-\log(1 - r/x) = \log(1 + r/(x - r))$ .

Elementary functions#

The following bound is the general principle; the theorems after it are its sharp special cases.

Theorem 6 (Generic first-order bound). Let $f$ be differentiable on a neighbourhood of the segment from $x$ to $x + a$ (real or complex). Then

|f(x + a) - f(x)| \le \min(2 C_0, C_1 |a|), \quad C_0 = \sup_{|t| \le |a|} |f(x + t)|, \quad C_1 = \sup_{|t| \le |a|} |f'(x + t)|.

Proof. The first bound is the triangle inequality, the second the mean value inequality $|f(x + a) - f(x)| = |\int_0^1 f'(x + sa)\, a\, ds| \le C_1 |a|$ .

Theorem 7 (Exponential). For real $x, a$ : $|e^{x + a} - e^x| = e^x |e^a - 1| \le e^x (e^{|a|} - 1)$ , with equality if $a \ge 0$ .

Proof. For $a \ge 0$ this is an identity. For $a < 0$ , $|e^a - 1| = 1 - e^{-|a|} \le e^{|a|} - 1$ , since $e^u + e^{-u} \ge 2$ .

arb_exp therefore adds $e^m (e^r - 1)$ , computed with mag_expm1 and an upper bound for the computed $e^m$ , to the radius.

Theorem 8 (Logarithm). For real $x, a$ with $0 \le |a| < x$ ,

|\log(x + a) - \log(x)| \le \log\left(1 + \frac{|a|}{x - |a|}\right),

with equality if $a \le 0$ .

Proof. The logarithm is concave and increasing, so the largest change is for $a = -|a|$ , where $\log x - \log(x - |a|) = \log(x/(x - |a|)) = \log(1 + |a|/(x - |a|))$ .

arb_log uses this bound with mag_log1p when the input has more than about 20 bits of relative accuracy, and otherwise evaluates $\log$ at both endpoints with the mag_t logarithm bounds.

Theorem 9 (Sine and cosine). For real $x, a$ : $|\sin(x + a) - \sin x| \le \min(2, |a|, |\cos x| |a| + |a|^2/2)$ , and the same holds for $\cos$ with $|\sin x|$ in place of $|\cos x|$ .

Proof. The bounds 2 and $|a|$ follow from Theorem 6. For the third, Taylor's theorem with the Lagrange remainder gives $\sin(x + a) = \sin x + a \cos x - \tfrac{a^2}{2} \sin \xi$ , and $|\sin \xi| \le 1$ .

arb_sin_cos uses the third form, with $|\cos x|$ bounded from an 8-bit table entry, which is much tighter than $|a|$ near the extrema of the sine. Wide balls (radius above about $2^{-24}$ ) are instead handled by evaluating at the endpoints and locating the extrema in between (arb_sin_cos_wide).

Theorem 10 (Arctangent). For real $x, a$ ,

|\operatorname{atan}(x + a) - \operatorname{atan}(x)| \le \min(\pi, C_1 |a|), \quad C_1 = \sup_{|t| \le |a|} \frac{1}{1 + (x + t)^2}.

If $|a| < |x|$ , one may take $C_1 = (1 + (|x| - |a|)^2)^{-1}$ . Exactly, if $|a| < |x|$ or $|x(x + a)| < 1$ , the difference equals $\operatorname{atan}(|a| / |1 + x(x + a)|)$ up to the choice of sign.

Proof. The first statement is Theorem 6 with $f' = 1/(1 + t^2)$ and the range of $\operatorname{atan}$ of length $\pi$ . If $|a| < |x|$ , every point between $x$ and $x + a$ has absolute value at least $|x| - |a|$ , which gives the monotone choice of $C_1$ . The exact form is the subtraction formula $\operatorname{atan} u - \operatorname{atan} v = \operatorname{atan}((u - v)/(1 + uv))$ , valid when $uv > -1$ .

arb_atan uses $|a| / (1 + (|x| - |a|)^2)$ , capped by $\pi$ , for accurate input and the endpoints otherwise.

Sums and series#

Many special functions are evaluated by truncating a series. The truncated sum is computed in ball arithmetic, and a bound for the tail is added to the radius. The following bounds are the most common.

Theorem 11 (Geometric bound). If $|c_k| \le C$ for all $k \ge N$ and $|z| \le D < 1$ , then

\left| \sum_{k=N}^{\infty} c_k z^k \right| \le \frac{C D^N}{1 - D}.

Proof. Compare termwise with the geometric series $\sum_{k \ge N} C D^k$ .

Theorem 12 (Integral bound). If $f$ is nonnegative and decreasing on $[N, \infty)$ , then

\int_N^{\infty} f(t)\, dt \le \sum_{k=N}^{\infty} f(k) \le f(N) + \int_N^{\infty} f(t)\, dt.

Proof. For $k \ge N$ , $f(k + 1) \le \int_k^{k+1} f \le f(k)$ ; sum over $k$ .

For example, mag_hurwitz_zeta_uiui bounds $\zeta(s, a) \le a^{-s} + a^{1-s}/(s - 1)$ in this way, and mag_exp_tail uses the geometric bound for the exponential series: if $N \ge 2|x|$ , consecutive terms of $\sum_{k \ge N} x^k/k!$ decrease by a factor at most $1/2$ , so the tail is at most twice its first term.

For analytic functions, the coefficients can be bounded by Cauchy's integral formula. If $f(z) = \sum_{k=0}^{\infty} c_k z^k$ is analytic on a neighbourhood of the disk $|z| \le R$ , then

c_k = \frac{1}{2\pi i} \int_{|z| = R} \frac{f(z)}{z^{k+1}}\, dz,

which gives the following two bounds.

Corollary 13 (Cauchy bounds). With $C = \max_{|z| = R} |f(z)|$ , the coefficients satisfy $|c_k| \le C / R^k$ , and for $0 \le r < R$ and $|z| \le r$ ,

\left| \sum_{k=N}^{\infty} c_k z^k \right| \le \frac{C D^N}{1 - D}, \quad D = \frac{r}{R}.

Proof. Bound the integrand by $C / R^{k+1}$ on a contour of length $2\pi R$ ; the tail bound is then Theorem 11 applied to $c_k R^k \cdot (z/R)^k$ .

The Euler-Maclaurin formula and its remainder bounds, used for the zeta function, are discussed in Zeta functions, polylogarithms and L-functions.

Use in FLINT#

arb_add, arb_sub, arb_mul, arb_addmul and their variants use the sum and Theorem 1; acb_mul applies Theorem 1 to complex balls in rectangular form componentwise.
arb_div and arb_inv use Theorem 2 or the endpoints.
arb_sqrt, arb_rsqrt and arb_root_ui use Theorem 3, Theorem 4 and Theorem 5.
arb_exp, arb_log, arb_sin_cos and arb_atan use the corresponding theorems, with endpoint evaluation for inaccurate input.
The tail bounds are implemented by mag_geom_series, mag_exp_tail, mag_hurwitz_zeta_uiui and mag_polylog_tail, and are used throughout arb_hypgeom, acb_hypgeom and the zeta function code.

void arb_mul(arb_t z, const arb_t x, const arb_t y, slong prec)

void arb_div(arb_t z, const arb_t x, const arb_t y, slong prec)

void arb_sqrt(arb_t z, const arb_t x, slong prec)

void arb_exp(arb_t z, const arb_t x, slong prec)

void arb_log(arb_t z, const arb_t x, slong prec)

void arb_sin_cos(arb_t s, arb_t c, const arb_t x, slong prec)

void arb_atan(arb_t z, const arb_t x, slong prec)