API - Arithmetic and Rounding

unit_roundoff(T) -> T

The unit roundoff u = eps(T)/2 of round-to-nearest, the bound on the relative error of one correctly rounded operation; for a complex type, that of its real type.

machine_epsilon(::Type{T}) where T <: AbstractFloat

Return the machine epsilon eps(T) of the floating-point type T, the spacing of the floats in [1, 2). It is TWICE the unit roundoff u of round-to-nearest (unit_roundoff), so machine_epsilon(T)·|c| bounds the error of one correctly rounded operation with a factor two to spare; the radii of the scalar ball operations are inflated with it.

For Float64, this equals 2^-52 ≈ 2.22e-16. For Float32, this equals 2^-23 ≈ 1.19e-7. For BigFloat, this is eps(BigFloat) = 2^(1 − precision), with the same factor two as the others (it was 2^(−precision), the unit roundoff itself, which left the complex operations short).

Examples

machine_epsilon(Float64)  # 2.220446049250313e-16
machine_epsilon(Float32)  # 1.1920929f-7
setprecision(256); machine_epsilon(BigFloat)  # ≈ 1.7e-77

subnormal_min(::Type{T}) where T <: AbstractFloat

Return the smallest positive subnormal value for floating-point type T. Adding this to a computed radius guarantees that results remain strictly positive even in edge cases where an operation underflows to zero.

For Float64, this equals 2^-1074 ≈ 4.94e-324. For Float32, this equals 2^-149 ≈ 1.40e-45. For BigFloat, this depends on the current precision and exponent settings.

Examples

subnormal_min(Float64)  # 4.9406564584124654e-324
subnormal_min(Float32)  # 1.0f-45

ϵp

Unit roundoff for Float64 values. This is a convenience constant equal to machine_epsilon(Float64) for backwards compatibility. For type-generic code, prefer machine_epsilon(T).

η

Smallest positive subnormal Float64 value. This is a convenience constant equal to subnormal_min(Float64) for backwards compatibility. For type-generic code, prefer subnormal_min(T).

abs_down(x::BigFloat)

A lower bound of |x| for a BigFloat, not an approximation of it: for x ≥ 0 the value is x, which is exact, and for x < 0 it is -x rounded toward -∞, so abs_down(x) ≤ |x| in either case. The upper bound is abs_up.

sum_up(xs) -> T

Upper bound of x₁ + x₂ + …, accumulated with add_up from zero in the order the iterable yields its elements. Every partial sum is an upper bound of the corresponding exact partial sum, so the result is an upper bound of the exact sum whatever that order is, and it is a bound and not an approximation: it is not the correctly rounded sum, and it is larger the more terms there are. An empty iterable of unknown element type gives 0.0.

pow_up(x, e) / pow_down(x, e)

Upper and lower bounds of xᵉ for x ≥ 0 and an integer e ≥ 0, by e directed products of a running value by x: each product is monotone in the running value because x ≥ 0, so the bound propagates. For x < 0 the sign alternates and neither result is a bound.

pow_down(x, e)

Lower bound of xᵉ for x ≥ 0 and an integer e ≥ 0, by e products rounded downward; see pow_up.

root_up(x, e) -> T

Upper bound of x^(1/e) for x ≥ 0 and an integer e ≥ 1: start from the floating-point root and step up by nextfloat until pow_down(r, e) ≥ x, which certifies rᵉ ≥ x and so r ≥ x^(1/e). The exit test is a certificate rather than an estimate, which is why the result is a bound; x ≤ 0 returns zero and x not finite returns Inf.

dist_up(a, b) / dist_down(a, b)

Upper and lower bounds of |a − b| for real or complex a, b. The difference is bracketed componentwise by directed subtractions, so the rounding of a − b is accounted for, and the modulus of the bracket is bounded with abs_up and abs_down.

dist_down(a, b)

Lower bound of |a − b| for real or complex a, b: each component of the difference is bounded below in modulus by directed subtractions, and the modulus of that pair with abs_down; see dist_up.

@up expr

Rewrite arithmetic in expr so that every operation rounds toward +∞: +, -, * and / are replaced by add_up, sub_up, mul_up and div_up. Each operation of the rewritten expression returns a bound and not an approximation, so the value of the whole expression is an upper bound of the exact one whenever every operation in it is monotone increasing in its arguments, which a subtraction and a division by a quantity that is itself only bounded are not; those are bracketed by hand, as dist_up does.

The rewriting is syntactic and rounds nothing by itself: @up on an expression whose operands are not floating-point numbers of a type these operations accept is an error rather than a loose bound. For BigFloat the named operations above are defined here and may be used directly.

@down expr

Mirror of @up that rewrites arithmetic so that each operation rounds toward -∞: add_down, sub_down, mul_down, div_down in place of +, -, *, /. The same caveat applies, that the result is a lower bound of the exact expression only where each operation is monotone increasing in its arguments.

add_up(x::BigFloat, y::BigFloat)

Add two BigFloat values with rounding toward +∞.

add_down(x::BigFloat, y::BigFloat)

Add two BigFloat values with rounding toward -∞.

sub_up(x::BigFloat, y::BigFloat)

Subtract two BigFloat values with rounding toward +∞.

sub_down(x::BigFloat, y::BigFloat)

Subtract two BigFloat values with rounding toward -∞.

mul_up(x::BigFloat, y::BigFloat)

Multiply two BigFloat values with rounding toward +∞.

mul_down(x::BigFloat, y::BigFloat)

Multiply two BigFloat values with rounding toward -∞.

div_up(x::BigFloat, y::BigFloat)

Divide two BigFloat values with rounding toward +∞.

div_down(x::BigFloat, y::BigFloat)

Divide two BigFloat values with rounding toward -∞.

sqrt_up(x::BigFloat)

Compute the square root of a BigFloat value with rounding toward +∞.

sqrt_down(x::BigFloat)

Compute the square root of a BigFloat value with rounding toward -∞.

add_up(x::Complex{BigFloat}, y::Complex{BigFloat})

Add two Complex{BigFloat} values with outward rounding for real and imaginary parts.

sub_up(x::Complex{BigFloat}, y::Complex{BigFloat})

Subtract two Complex{BigFloat} values with outward rounding.

mul_up(x::Complex{BigFloat}, y::Complex{BigFloat})

Multiply two Complex{BigFloat} values with outward rounding. Uses the formula (a+bi)(c+di) = (ac-bd) + (ad+bc)i with careful rounding.

two_sum(a, b) -> (x, y)

The exact sum as an unevaluated sum: x = fl(a+b) and y = a+b - x exactly, so that

a+b == x + y,

by Knuth's six operations, with no assumption on the relative size of a and b. The identity holds in the presence of underflow as well, and |y| ≤ u|x| and |y| ≤ u|a+b| with u the unit roundoff (unit_roundoff); only overflow is excluded.

Reference

Algorithm 3.1 (TwoSum) and Theorem 3.4 of [11], where the transformation is attributed to Knuth, The Art of Computer Programming, vol. 2, Seminumerical Algorithms ([12]).

two_product(a, b) -> (x, y)

The product as an unevaluated sum: x = fl(a*b) and y the residual a*b - x computed by one fused multiply-add, so that

a*b == x + y

barring underflow and overflow. Under gradual underflow the residual need not be representable, and y is that residual rounded once, so x + y differs from a*b by at most half the spacing of the subnormals, subnormal_min(T)/2; that statement is ours, since the reference states its underflow bound for the version built on Dekker's splitting and not for the one with the fused multiply-add.

Reference

Algorithm 3.5 (TwoProductFMA) of [11]; the error-free statement for the splitting version is its Theorem 3.4.

split_veltkamp(a) -> (hi, lo)

The splitting of Dekker:

a == hi + lo

exactly, with each half carrying about half the significand, so that products of halves are exact in floating point. Both the splitting and the exactness of those products exclude overflow of the scaled value (2^s + 1)·a and underflow of the halves; where either occurs the identity above is no longer the exact one. This is the step that lets an accurate matrix product be assembled from ordinary BLAS calls in round-to-nearest, with no directed rounding anywhere.

Reference

Algorithm 3.2 (Split) of [11], with factor = 2^s + 1 and s = ⌈t/2⌉ for t bits of precision; there it is attributed to Dekker [13], and the product of two numbers by it (Algorithm 3.3, TwoProduct) to G. W. Veltkamp.

gamma_bound(N, ::Type{T}) -> T

Higham's gamma_N = N u / (1 - N u), rounded upward, for N accumulated operations: the constant of |fl(xᵀy) − xᵀy| ≤ γ_N |x|ᵀ|y| for a real inner product of length N in any order, with or without fused multiply-add. Refuses when N u ≥ 1/2, where the model does not apply.

For complex data, gamma_bound(N, Complex{T}) returns γ_{N+2}: a complex product has relative error at most √2 γ₂ ≤ γ₃ (Lemma 3.5 of the reference), and (1 + γ₃)(1 + γ_{N−1}) ≤ 1 + γ_{N+2}.

Reference

[14]: inequality (3.5) for the inner product, Section 2.6 for the fused multiply-add (n rounding errors in place of 2n − 1), Lemma 3.5 for complex arithmetic.

compensated_terms(pairs, i, j) -> T

One entry of a sum of matrix products sum_t sgn_t * A_t[i,:] * B_t[:,j], accumulated with error-free transformations so that the result carries about twice the working precision. pairs is an iterable of (A, B, sgn).

This is Algorithm 5.3 (Dot2) of [11] applied to a sum of products. It stands here where the papers of Rump call INTLAB's prodK with k = 2; prodK itself is a different algorithm, see src/types/MMul/ogita_rump_oishi.jl. It returns the accurate value only; a certified bound on it is Proposition 5.5 of the same paper,

|result - exact| <= u |exact| + gamma_N^2 sum |a_i| |b_i|,    gamma_N = N u / (1 - N u),

with N the number of products summed and the absolute products evaluated separately, plus the term 5 N eta of that proposition in the presence of underflow, eta the underflow unit (subnormal_min).

compensated_terms2(pairs, i, j) -> (s, e)

The two terms of compensated_terms before they are added: s is the floating-point sum of the products and e the accumulated sum of the errors of the products and of the additions, so that s + e, unevaluated, is the entry to about twice the working precision.

With N the number of products and Σ the sum of their moduli, barring underflow,

|s + e − exact| ≤ γ_N γ_{N+1} Σ ≤ γ_{N+1}² Σ.

Proof. two_product and two_sum are exact, so exact = s + Σ_k (ep_k + es_k) with ep_k, es_k the two errors of step k. e is the floating-point sum of the N numbers fl(ep_k + es_k), one rounding each and N − 1 additions, so |e − Σ_k (ep_k + es_k)| ≤ γ_N Σ_k |ep_k + es_k|. And |ep_k| ≤ u(1 + u)|a_k b_k|, |es_k| ≤ u|s_k| ≤ u(1 + u)^k Σ_{i ≤ k}|p_i|, whence Σ_k |ep_k + es_k| ≤ (N + 1)u(1 + u)^{N+1} Σ ≤ γ_{N+1} Σ.

The exactness of two_product is where underflow enters, and the a priori bound of Proposition 5.5 of [11] adds 5 N eta for it; _accuracy_term in src/types/MMul/ogita_rump_oishi.jl is what charges that term.

mmul_ogita_rump_oishi_2005(A::BallMatrix, B::BallMatrix) -> BallMatrix

Product of two ball matrices whose midpoint is accumulated with error-free transformations and whose radius charges the accuracy bound of [11] rather than the rounding of a BLAS product.

For every X ∈ A and Y ∈ B the returned ball contains XY. The radius is the sum of

• the interval part, rad(A)(|mid(B)| + rad(B)) + |mid(A)| rad(B) entrywise, which is what any midpoint-radius product must carry and is identical to MMul4's; and
• the accuracy part, (u|C| + κ(γ_N² Σᵢ|aᵢ||bᵢ| + 5N·eta))/(1 − u) with N the number of products summed, γ_N = Nu/(1−Nu), eta the underflow unit and κ = 1 for real data, 3/2 for complex: Proposition 5.5 of the reference with its underflow term, see _accuracy_term.

The γ_N² is the point: MMul3, MMul4 and MMul5 all charge the accumulation over the k summands at order γ_k, and at k = 5000 in Float64 that is 1.1e-12 against this routine's 1.2e-28.

The midpoint never changes the rounding mode; two_product uses one FMA and two_sum is Knuth's six additions, both in round-to-nearest. Only the radius is accumulated under RoundUp.

Cost. O(mnk) scalar operations with no BLAS, threaded over the output columns. Measured on 16 threads against the ball * on complex matrices with radius zero, it is faster below about k = 400 — 0.55× at k = 200, 0.41× at 400 — and 1.32× at 800, 3.59× at 1600, with the radius ratio improving from 0.0075 to 0.0029 over the same range. * dispatches to MMul4, which is the default product of the package.

Reference

[11], Proposition 5.5 for the accuracy bound and Algorithm 5.3 (Dot2) for the accumulation.

_dot2err_ogita_rump_oishi_2005(x, y) -> (res, err)

Algorithm 5.8 (Dot2Err) of [11] for two vectors of Float32 or Float64: res is the dot product accumulated in twice the working precision and rounded, and res − err ≤ xᵀy ≤ res + err (Corollary 5.9), also in the presence of underflow. Every operation is in rounding to nearest, so the pair is a bound and not an approximation whatever the BLAS does with the rounding mode. Where an intermediate result is not finite err is Inf.

Reference

[11], Algorithm 5.8 and Corollary 5.9.

_product_down_up_dot2err(A, B) -> (down, up)

Entrywise lower and upper bounds of the product of two Float32 or Float64 matrices, down ≤ AB ≤ up, with no change of the rounding mode: each entry is accumulated by Algorithm 5.8 (Dot2Err) of [11] as an unevaluated sum p + s with the bound β of the proof of Corollary 5.9 on its distance from the exact entry, underflow included, and p + (s − β), p + (s + β) are rounded downward and upward by add_down, sub_down, add_up. An entry with a non-finite intermediate result gets (-Inf, Inf).

The same role as _product_down_up, which evaluates the BLAS product twice under RoundDown and RoundUp. O(mnk) scalar operations, no BLAS, threaded over the columns.

Reference

[11], Algorithm 5.8 and Corollary 5.9.

_two_term_product_sum(terms) -> (P1, P2, R)

Σ_t sgn_t A_t B_t for terms an iterable of (A_t, B_t, sgn_t), point matrices, as an unevaluated sum of two matrices: |Σ_t sgn_t A_t B_t − P1 − P2| ≤ R entrywise. P1 is the floating-point value and P2 the accumulated error terms of compensated_terms2, so that a further product with the pair keeps twice the working precision; this is the output of the role of prodK(..., 'OutputTerms', 2) with k = 2 in [15].

For Float32 and Float64 the splitting scheme of [16], which is prodK itself, is tried first through _lange_rump_two_terms; where it declines, and for every other type, each entry is one Dot2 accumulation of [11] with the a priori bound

R = κ(γ_{N+1}² Σ_t |A_t||B_t| + 5N eta),

N the number of products in a part of an entry, κ = 1 for real data and 3/2 ≥ √2 for complex, where the bound holds for the real and for the imaginary part: the bound of compensated_terms2, with the underflow term of _accuracy_term.

_accurate_product_sum(terms) -> BallMatrix

A ball matrix containing Σ_t sgn_t A_t B_t for terms an iterable of (A_t, B_t, sgn_t), point matrices: the two terms of _two_term_product_sum added, the radius its bound plus the rounding of that addition. The role of [R, E] = prodK(...) with k = 2 in [15], whose splitting scheme is [16].

_lange_rump_two_terms(terms; k = 2) -> (P1, P2, R) or nothing

Σ_t sgn_t A_t B_t for terms an iterable of (A_t, B_t, sgn_t), point matrices of Float32 or Float64 entries, real or complex, and sgn_t = ±1, as an unevaluated sum of two matrices with |Σ_t sgn_t A_t B_t − P1 − P2| ≤ R entrywise (the modulus, for complex data).

The method is the splitting scheme of [16], Accurate floating-point matrix residuals, an unpublished manuscript that is the algorithm of INTLAB's prodK: each factor is split into k leading parts and a remainder by their Algorithms 1 to 4, the rows of A_t and the columns of B_t at the points (25) with q_a = fl(√2^p), q_b = fl((1 − u)2^p/q_a) (Algorithm 5), and the product is their equation (27): k² matrix products that are exact by their Lemma 4 and Corollary 5, and two rounded remainder products. It corresponds to prodK(…, k, 'OutputTerms', 2), "as if" in (k/2 + 1)-fold precision. A summand with a diagonal factor is formed entry by entry with two_product. The summands are split separately (Section 8.2 of the paper).

R is the sum of

• for each summand, with n the inner dimension, γ(|A − A̲| |B̲| + |A̲| |B|) + |δ||B̲| bounded through |B̲| ≤ ν⁽ᵏ⁾/2 and |A̲| ≤ μ⁽ᵏ⁾/2 entrywise (1/√2 in place of 1/2 for complex data), γ = γ_n (γ_{n+2} for complex), δ the rounding of fl(A − A̲), and an underflow term;
• the rounding of the accumulation of the two_sum errors of the final summation.

Returns nothing, and the caller falls back on the accumulation by compensated_terms2, where a sign is not ±1, where a splitting constant overflows or a product of splitting points underflows, where 2^p does not exceed the squared inner dimension by enough for Corollary 5, or where the result is not finite.

The exactness of the leading products requires a matrix product that only adds and multiplies (as the reference BLAS and OpenBLAS do); a Strassen-type product would not qualify (footnote 3 of the paper).

Not implemented from the paper: λ < 1, unbalanced splittings, the Karatsuba evaluation, the blocked summation of Section 8.1, multiple-precision input (Algorithm 7), more than two output terms, sparse factors.

oishi_MMul(F, G)

High-level wrapper: compute a complex ball enclosure of the product F*G by the rounding-mode controlled product of [4]. Internally calls _oishi_MMul_up_lo to get componentwise rectangular bounds, then _ccr to convert them into a BallMatrix, so the returned ball contains F*G entrywise.

Correctness rests on the working type honouring setrounding, which for Float32 and Float64 means the BLAS; see BallArithmetic.NumericalTest.rounding_test.

_oishi_MMul_up_lo(F::AbstractMatrix{<:Complex}, G::AbstractMatrix{<:Complex})

Compute outward lower/upper enclosures for the real and imaginary parts of the complex matrix product F*G using directed-rounding (rounding-mode controlled) evaluation, Algorithm 2 (cprod) of [5], which the paper attributes to [4].

Returns (Hrl, Hru, Hil, Hiu, T), where Hrl ≤ Re(F*G) ≤ Hru and Hil ≤ Im(F*G) ≤ Hiu hold entrywise, and T is the promoted underlying real working type. Each of the four is a bound and not an approximation: the products that make up a lower bound are evaluated under RoundDown and those of an upper bound under RoundUp.

Promotion follows T = promote_type(float(real(eltype(F))), float(real(eltype(G)))). For mathematical correctness, T must honor setrounding (e.g., BigFloat).

_cprod(F::AbstractMatrix{<:Complex}, G::AbstractMatrix{<:Complex})

The product F G of two complex floating point matrices, enclosed by its lower and upper bounds in each of the real and the imaginary part, as _oishi_MMul_up_lo returns them: Hrl, Hru, Hil, Hiu and the type the products were taken in. Each bound comes from one directed accumulation, so the four together enclose the exact product.

The name is the one [5] uses for the complex product computed with directed rounding, and this is the helper the Miyajima kernels below call for it.

_cr(Fl::AbstractMatrix{<:Real}, Fu::AbstractMatrix{<:Real}, ::Type{T}) where {T}

Advanced API — Algorithm 5 of Ref. [5]. Convert entrywise real lower/upper bounds Fl ≤ F ≤ Fu into midpoint and radius matrices (Fc, Fr) evaluated in type T using directed rounding.

Useful in conjunction with _ccr and the interval products _iprod / _ciprod.

_cr(Fl::AbstractMatrix{<:Real}, Fu::AbstractMatrix{<:Real})

Convenience overload returning (Fc, Fr, T) where T is the promoted working type.

_ccr(Hrl, Hru, Hil, Hiu, ::Type{T}) where {T<:AbstractFloat}

Advanced API — collapse rectangular enclosures for Re(F*G) and Im(F*G) to a complex ball enclosure BallMatrix(Hc, Hr).

Algorithm 3 (ccr) of [5]: the real and the imaginary bounds are each converted by Algorithm 5 (_cr), Rc = fl△(Hrl + 0.5(Hru − Hrl)), Rr = fl△(Rc − Hrl) and the same for the imaginary part, and the radius is Hr = fl△(|Rr + i Ir|), so that [Hrl, Hru] + i[Hil, Hiu] ⊆ ⟨Rc + i Ic, Hr⟩.

The rectangular bounds this collapses come from the rounding-mode controlled product of [4].

See also _ccrprod and _ccrprod_prime.

_ccr(Hrl::AbstractMatrix{<:Real}, Hru::AbstractMatrix{<:Real},
     Hil::AbstractMatrix{<:Real}, Hiu::AbstractMatrix{<:Real})

Convenience overload returning (B, T) where B is the BallMatrix produced by _ccr(Hrl, Hru, Hil, Hiu, T) and T is the promoted working type.

_ccrprod(J::AbstractMatrix{<:Complex}, Hc::AbstractMatrix{<:Complex}, Hr::AbstractMatrix{<:Real})

Advanced API — Algorithm 4 from Ref. [5]. Given a complex matrix J and a ball enclosure (Hc, Hr) for another complex matrix, return rectangular bounds (Krl, Kru, Kil, Kiu) and the working type T such that Krl ≤ Re(J * (Hc ± Hr)) ≤ Kru and Kil ≤ Im(J * (Hc ± Hr)) ≤ Kiu hold entrywise. All computations are performed with outward rounding using the promoted working type T.

See also _ccrprod_prime and _ccr for the right-multiplication variant and the conversion to ball form.

_ccrprod_prime(Hc::AbstractMatrix{<:Complex}, Hr::AbstractMatrix{<:Real},
               J::AbstractMatrix{<:Complex})

Advanced API — right-multiplication analogue of Algorithm 4 from Ref. [5]. Given a complex ball enclosure (Hc, Hr) and a complex matrix J, return rectangular bounds (Krl, Kru, Kil, Kiu) and the promoted working type T such that 〈Hc, Hr〉 * J ⊆ [Krl, Kru] + √(-1) [Kil, Kiu] entrywise.

See also _ccrprod and _ccr.

_iprod(F::AbstractMatrix{<:Real}, Gc::AbstractMatrix{<:Real}, Gr::AbstractMatrix{<:Real})

Advanced API — Algorithm 6 of Ref. [5]. Multiply a real matrix F by a ball enclosure (Gc, Gr) and return rectangular bounds (Hl, Hu) and the working type T enclosing F * (Gc ± Gr).

The helper complements _ccrprod for mixed real/complex interval arithmetic.

_iprod_right(Gc::AbstractMatrix{<:Real}, Gr::AbstractMatrix{<:Real}, F::AbstractMatrix{<:Real}, ::Type{T}) where {T}

Advanced API — right-multiplication analogue of _iprod. Given a real ball enclosure (Gc, Gr) and a real matrix F, compute rectangular bounds (Hl, Hu) for (Gc ± Gr) * F using directed rounding in the working type T.

Pairs naturally with _iprod when traversing the Miyajima 2010 verification procedures.

_ciprod(J::AbstractMatrix{<:Complex}, Hrl, Hru, Hil, Hiu)

Advanced API — Algorithm 7 of Ref. [5]. Multiply a complex matrix J by rectangular bounds on another complex matrix, provided as lower/upper bounds for the real and imaginary parts. The result is returned as rectangular bounds (Krl, Kru, Kil, Kiu) together with the working type T.

Combines with _ccrprod / _ccrprod_prime to mirror the flow in Miyajima's block diagonalisation routines.

_ciprod_prime(Hrl, Hru, Hil, Hiu, J)

Advanced API — right-multiplication analogue of Algorithm 7 from Ref. [5]. Given rectangular bounds for the real and imaginary parts of a complex matrix and an additional complex matrix J, compute rectangular bounds for ([Hrl, Hru] + √(-1)[Hil, Hiu]) * J with outward rounding in the promoted type T.

Complements _ciprod for the Miyajima 2010 verification flow.

ball_fft(v::BallVector; method = :brisebarreetal2023)

The discrete Fourier transform y_k = Σ_j x_j exp(−2πi jk/N) of every vector in the ball v, of length N = 2ⁿ: a ball vector whose midpoint is the transform of the midpoint computed in floating point and whose radii are all

√2 · b_n · ‖X‖∞⊥ + ‖r‖₁.

Here ‖X‖∞⊥ is the largest absolute value of a real or imaginary part of the midpoints, r the radii of v, and b_n the constant of equation (4) of [17], which bounds each part of the rounding error of an output by b_n ‖X‖∞⊥; the factor √2 passes from the two parts to the modulus. The term ‖r‖₁ is the effect of the radii: an output is a sum of the inputs times numbers of modulus one.

Each algorithm is named for the paper it implements, a deviation being visible in the name:

Why the bound is not known to hold for FFTW. The paper's proof is for the radix-2 step with the roots of unity rounded to nearest, and FFTW chooses among other algorithms and computes its roots in its own way. An audit of the FFTW source, to find which bound holds for what it computes, is in progress and not published; until it is, the radius of this method is a hypothesis.

The paper's analysis excludes underflow and overflow; an output that is not finite is an ArgumentError, and underflow is not detected. A length that is not a power of two is an ArgumentError. With FFTW loaded, fft(v) is ball_fft(v).

Reference

[17], equation (4) for b_n and Fig. 1 for the algorithm the bound is proved for.

ball_fft(A::BallMatrix, dims = (1, 2); method = :brisebarreetal2023)

The discrete Fourier transform of the ball matrix A along the dimensions dims (an integer, or a tuple), each of size a power of two. Along one dimension it is ball_fft of each column or row, with the bound of that fiber. Along both, the columns are transformed and then the rows, and every output has the radius

√2 · b_n · ‖X‖∞⊥ + ‖r‖₁,        n = log₂(number of entries),

with the norms taken over the whole matrix.

The bound along both dimensions is not in the paper; it follows from the bound it quotes for the relative error in the Euclidean norm of the one-dimensional algorithm, ε_m = (1 + u)ᵐ ∏_{j ≤ m} (1 + g_j) − 1 for size 2ᵐ. Applied to the columns and summed over them, and then to the rows, it gives ‖Ŷ − Y‖_F ≤ ((1 + ε_{n₁})(1 + ε_{n₂}) − 1) ‖Y‖_F, the exact transform multiplying the Frobenius norm by the square root of the size of a fiber; and (1 + ε_{n₁})(1 + ε_{n₂}) ≤ 1 + ε_n for n = n₁ + n₂. An entry is at most the Frobenius norm, and ‖Y‖_F = √N ‖X‖_F ≤ √2 N ‖X‖∞⊥ with N = 2ⁿ, which is b_n ‖X‖∞⊥.

method is as for a vector.

ball_fft(v::BallVector, plan::BallFFTPlan)

ball_fft with the table of plan, which must have been made for the length of v.

plan_ball_fft(N; T = Float64, method = :brisebarreetal2023)

The BallFFTPlan for transforms of length N, a power of two: it computes the table of roots once, which is the cost a repeated transform should not pay twice. T is the float type of the midpoints.

BallFFTPlan{T}

A size and the table of roots of unity it needs, so that many transforms of that size share one table. The table of order N holds the roots of every smaller power of two as well, which is why a plan made for N transforms any size dividing it.

Build one with plan_ball_fft and pass it to ball_fft. Keeping it is what makes a large size affordable twice: at N = 2²⁸ the table is 2²⁷ complex numbers, 2.1 GB for Float64, and some seconds of BigFloat evaluation to build. Dropping the plan frees it.

_fftw_transform(x::Vector{<:Complex})

The discrete Fourier transform of x in floating point by FFTW. A function without methods in the package; the extension for FFTW adds it. It is the midpoint of the method :brisebarreetal2023fftw of ball_fft.

PolyRange{T}

Certified range enclosure of a polynomial on [a,b]. With min/max the true extrema of p over [a,b], the four fields are sound brackets

min ∈ [min_lo, min_hi],    max ∈ [max_lo, max_hi],

so the true range [min, max] ⊆ [min_lo, max_hi] (the outer enclosure), and each extremum is pinned to a bracket of width min_hi-min_lo / max_hi-max_lo (the residual overestimation). method is :bisection or :critical; diag carries the method-specific cost/quality record.

ball_horner(coeffs, x::Ball) -> Ball

The polynomial p(t) = Σₖ coeffs[k] t^(k−1) (coefficients in ascending order) evaluated by Horner's scheme in ball arithmetic at the ball x, real or complex: every product and sum is a ball operation, so the returned ball contains p(t) for every t in x.

range_bisect(coeffs, a, b; tol = 1e-8, maxeval = 100_000) -> PolyRange

Algorithm A — certified range of p (ascending coeffs) on [a,b] by interval branch-and-bound. tol is the target width of each extreme's bracket; maxeval caps the ball evaluations. The maximum is obtained from -p.

range_critical(coeffs, a, b; maxiter = 20) -> PolyRange

Algorithm B — certified range of p (ascending coeffs) on [a,b] by enclosing the stationary points. Forms p', encloses the roots of its companion matrix with verifyeigall (maxiter interval iterations), keeps every inclusion disc that meets the real axis inside [a,b], and hulls the ball-evaluations of p at the endpoints and those stationary enclosures.

The candidate set is exhaustive only when every cluster is certified (spectrum_covered); when it is not, the range is computed by range_bisect instead and diag.covered is false.

diag reports ncrit (real stationary enclosures used), maxdisc (largest root-disc radius: large means the companion was hard to enclose, and the bound loosens but stays rigorous), and covered.