Outward rounded interval and ball arithmetic

A Dedekind-complete ordered field is a commutative field $\mathbb{R}$ with identity $1>0$, a total order compatible with addition and multiplication, and the least-upper-bound axiom: every nonempty subset of $\mathbb{R}$ that is bounded above has a least upper bound in $\mathbb{R}$. Throughout this batch, all real variables, intervals, limits, derivatives and integrals are taken in such a field. The integers and rationals are interpreted in their usual embedded copies. For $x\in \mathbb{R}$, define $\mathrm{\mid }x\mathrm{\mid }=x$ when $x\ge 0$ and $\mathrm{\mid }x\mathrm{\mid }=-x$ when $x<0$.

For $f:[a,b]\to \mathbb{R}$, continuity at $x\in [a,b]$ means that for every $\epsilon >0$ there is $\delta >0$ such that $\mathrm{\mid }f(t)-f(x)\mathrm{\mid }<\epsilon$ whenever $t\in [a,b]$ and $\mathrm{\mid }t-x\mathrm{\mid }<\delta$. Differentiability at an interior point $c$ means that the two-sided limit

exists in the real field; differentiability on $(a,b)$ means this holds at every interior point. For higher derivatives, define recursively $f^{(j+1)}=(f^{(j)}{)}^{\mathrm{\prime }}$ on the points where the derivative exists. The class $C^{n+1}$ on an open interval means that these derivatives through order $n+1$ exist and are continuous there.

An exact dyadic rational is a number $m{2}^q$ with integers $m,q$. Fix integers $p\ge 2$ and $e_{\min }\le e_{\max }$. Put $\eta =2^{e_{\min }-p+1}$ and $N_{\max }=(2^p-1)2^{e_{\max }-p+1}$. The finite representable set $\mathbb{F}$ consists of zero, both signs of every normal number $k{2}^{e-p+1}$ with integers $e,k$ satisfying $e_{\min }\le e\le e_{\max }$ and $2^{p-1}\le k<2^p$, and both signs of every subnormal number $k\eta$ with integer $k$ satisfying $1\le k<2^{p-1}$. The smallest positive subnormal is $\eta$, the smallest positive normal is $2^{e_{\min }}$, and the largest finite magnitude is $N_{\max }$. The extended endpoint set also contains $-\mathrm{\infty }$ and $+\mathrm{\infty }$ as enclosure sentinels; they are not finite machine numbers.

For a real $x$ with $\mathrm{\mid }x\mathrm{\mid }\le N_{\max }$, define $\mathrm{rd}(x)$ as the greatest member of $\mathbb{F}$ not exceeding $x$, and $\mathrm{ru}(x)$ as the least member of $\mathbb{F}$ not less than $x$. Define $\rho (x)=\mathrm{ru}(x)-\mathrm{rd}(x)$. Set ${ϵ}_{\mathrm{a}\mathrm{b}\mathrm{s}}(0)=0$. For nonzero $x$ in the subnormal range, define the absolute-error bound ${ϵ}_{\mathrm{a}\mathrm{b}\mathrm{s}}(x)=\eta$; for $2^e\le \mathrm{\mid }x\mathrm{\mid }<2^{e+1}$ in the normal range, define ${ϵ}_{\mathrm{a}\mathrm{b}\mathrm{s}}(x)=2^{e-p+2}$. For $x>N_{\max }$, set $(\mathrm{rd}(x),\mathrm{ru}(x))=(N_{\max },+\mathrm{\infty })$; for $x<-N_{\max }$, set $(\mathrm{rd}(x),\mathrm{ru}(x))=(-\mathrm{\infty },-N_{\max })$. In either overflow case set ${ϵ}_{\mathrm{a}\mathrm{b}\mathrm{s}}(x)=+\mathrm{\infty }$ and report overflow, even when the selected directed endpoint is finite. The rounding procedure reports finite, subnormal, or overflow, so no relative-error assumption is implicit.

An exact rational intermediate may be converted by comparing integer numerators and denominators after clearing powers of two. A directed interval endpoint is the resulting greatest lower or least upper member of the extended endpoint set $\mathbb{F}\cup \{-\mathrm{\infty },+\mathrm{\infty }\}$.

In a Dedekind-complete ordered field, the natural numbers are unbounded above.

If $a<b$ in a Dedekind-complete ordered field, every relative open cover of $[a,b]$ has a finite subcover. Relative openness means openness in the subspace topology inherited from the real line.

In a Dedekind-complete ordered field, every $x\ge 0$ has a unique $y\ge 0$ with $y^2=x$.

Input: a nonnegative exact dyadic rational $x$ and integer request $n\ge 1$. Output: an exact dyadic lower and upper bound for $\sqrt{x}$; no finite storage format is needed by this procedure.

1. If $x=0$, return $[0,0]$.
2. Otherwise determine the unique integer $j$ with $2^j\le x<2^{j+1}$ from the exact dyadic representation, put $m=\lfloor j\mathrm{/}2\rfloor$ and $\delta =2^{m-n}$.
3. Set integer indices $L=0$ and $U=2^{n+1}$. While $U-L>1$, put $Q=\lfloor (L+U)\mathrm{/}2\rfloor$. If $(Q\delta {)}^2\le x$, replace $L$ by $Q$; otherwise replace $U$ by $Q$.
4. Return the exact dyadic pair $[L\delta ,U\delta ]$. A caller that stores these endpoints in a selected finite format must separately round the lower endpoint down and the upper endpoint up, propagating the resulting subnormal or overflow status.

All comparisons and multiplications in the loop are exact integer/dyadic operations.

A real interval is a closed set $[a,b]\subset \mathbb{R}$ with finite dyadic endpoints $a\le b$. Its width is $b-a$, midpoint is $(a+b)\mathrm{/}2$, and its midpoint-radius form is $m+[-r,r]$ with $m=(a+b)\mathrm{/}2$ and $r=(b-a)\mathrm{/}2$. The empty set is the status empty; the whole real line is entire. overflow, domain-error, split-domain, and inconclusive are computation statuses, not finite intervals.

A complex rectangle is $X+iY=\{x+iy:x\in X,y\in Y\}$ for real intervals $X,Y$. Complex addition and multiplication use $(x+iy)+(u+iv)=(x+u)+i(y+v)$ and $(x+iy)(u+iv)=(xu-yv)+i(xv+yu)$. For $z=x+iy$, its complex modulus is $\mathrm{\mid }z\mathrm{\mid }=\sqrt{x^2+y^2}$, and for $z\mathrm{\ne }0$ its reciprocal is $1\mathrm{/}z=(x-iy)\mathrm{/}(x^2+y^2)$. A complex ball $B(c,r)$ is $\{z\in \mathbb{C}:\mathrm{\mid }z-c\mathrm{\mid }\le r\}$ with dyadic complex center $c$ and dyadic $r\ge 0$. For any two represented sets $S,T$, the inclusion order is $S⊑T$ iff $S\subseteq T$.

The width of a complex rectangle with coordinate widths $w_X,w_Y$ is its Euclidean diameter $\sqrt{w_X^2+w_Y^2}$; the width of a ball $B(c,r)$ is its diameter $2r$. Width comparisons against a dyadic tolerance can be decided by exact rational comparisons after squaring the nonnegative quantities.

The range of a real or complex function $f$ on a set $S$ is $f(S)=\{f(x):x\in S\cap \mathrm{dom}(f)\}$. A finite expression DAG is a finite directed acyclic graph with input nodes, constant nodes, and operation nodes; constants are exact dyadics or are supplied with certified enclosing intervals. Sharing a node means the exact same intermediate value is reused.

For finite $x$ with $\mathrm{\mid }x\mathrm{\mid }\le N_{\max }$, the directed roundings defined in Dyadic floating-point format and directed rounding satisfy $\mathrm{rd}(x)\le x\le \mathrm{ru}(x)$ and each endpoint error is at most the stated ${ϵ}_{\mathrm{a}\mathrm{b}\mathrm{s}}(x)$ (with zero error when $x$ is representable). For $x>N_{\max }$ the pair $(N_{\max },+\mathrm{\infty })$ encloses $x$; for $x<-N_{\max }$ the pair $(-\mathrm{\infty },-N_{\max })$ encloses $x$. Overflow has error bound $+\mathrm{\infty }$ and cannot be reported as a finite-error result.

For real intervals $X=[a,b]$ and $Y=[c,d]$, compute the displayed exact endpoint formulas and round the lower endpoint down and the upper endpoint up in the format of Dyadic floating-point format and directed rounding. If either exact endpoint overflows, return overflow with the outward extended endpoints; never substitute a finite endpoint for an overflow sentinel. Addition is $[a+c,b+d]$, subtraction is $[a-d,b-c]$, and multiplication is $[\min (ac,ad,bc,bd),\max (ac,ad,bc,bd)]$. Division is defined only when $0\mathrm{\notin }Y$ and is $X[1\mathrm{/}d,1\mathrm{/}c]$, with reciprocal endpoints computed as exact rationals before outward rounding. Squaring is $[0,\max (a^2,b^2)]$ when $0\in X$ and $[\min (a^2,b^2),\max (a^2,b^2)]$ otherwise. Square root is defined only for $X\subseteq [0,\mathrm{\infty })$ and is enclosed by applying Exact dyadic bisection for square-root endpoints to the exact dyadic endpoints with request $n=p$, where $p$ is the selected format precision, and rounding outward for storage.

If $0\in Y$ for division, return split-domain with the two subdomains $Y\cap (-\mathrm{\infty },0)$ and $Y\cap (0,\mathrm{\infty })$, or entire when the caller does not request subdivision. No finite quotient interval is returned in this case. An invalid square-root domain returns domain-error or an explicitly restricted nonnegative subdomain.

Complex-rectangle addition, subtraction and multiplication use $(x+iy)(u+iv)=(xu-yv)+i(xv+yu)$, applying the real interval operations to each component. For a rectangle $X+iY$ excluding $0$, define its reciprocal by $(X-iY)\mathrm{/}(X^2+Y^2)$, computing the square hulls, their sum, and the component quotient endpoint formulas as exact rationals before rounding the final reciprocal components outward. The exact denominator interval is strictly positive because at least one coordinate interval excludes zero; do not round it before division, since underflow could introduce zero. Propagate any final endpoint overflow. Division multiplies by this reciprocal.

For a ball $B(c,r)$, let $X_B=[\mathrm{\Re }c-r,\mathrm{\Re }c+r]$ and $Y_B=[\mathrm{\Im }c-r,\mathrm{\Im }c+r]$ be its coordinate projection intervals. For a real interval $J$, let $M(J)=\max (\mathrm{\mid }\inf J\mathrm{\mid },\mathrm{\mid }\sup J\mathrm{\mid })$ and let $d(J)=0$ if $0\in J$, otherwise $d(J)=\min (\mathrm{\mid }\inf J\mathrm{\mid },\mathrm{\mid }\sup J\mathrm{\mid })$. Put $M_B=M(X_B)+M(Y_B)$ and $m_B=\max (d(X_B),d(Y_B))$. Then $\mathrm{\mid }z\mathrm{\mid }\le M_B$ for every $z\in B(c,r)$ and $\mathrm{\mid }z\mathrm{\mid }\ge m_B$ for every such $z$. Define ball addition and subtraction as $B(0,\mathrm{ru}(M_{B_1}+M_{B_2}))$, multiplication as $B(0,\mathrm{ru}(M_{B_1}{M}_{B_2}))$, and division as $B(0,\mathrm{ru}(M_{B_1}\mathrm{/}{m}_{B_2}))$ only when $m_{B_2}>0$. These origin-centered output balls are conservative but inclusion-isotone. If $m_{B_2}=0$, return split-domain, entire, or inconclusive; if a radius overlows, return overflow.

Dyadic endpoints are rational, so every real arithmetic endpoint and every complex-division component before square root is rational. Reciprocals are formed only for nonzero rational endpoints. In a complex rectangle excluding zero, at least one coordinate interval omits zero; its exact squared interval has a strictly positive lower endpoint, so the unrounded $X^2+Y^2$ does also. The prescribed exact rational component quotients therefore have nonzero denominators even when storage rounding would underflow that lower bound. For ball division, $m_{B_2}>0$ gives $\mathrm{\mid }y\mathrm{\mid }\ge m_{B_2}>0$ for every $y\in B_2$. Every endpoint overflow is reported explicitly and cannot produce a finite-result claim.

For real $x,y$, $\mathrm{\mid }xy\mathrm{\mid }=\mathrm{\mid }x\mathrm{\mid }\mathrm{\mid }y\mathrm{\mid }$ and $\mathrm{\mid }x+y\mathrm{\mid }\le \mathrm{\mid }x\mathrm{\mid }+\mathrm{\mid }y\mathrm{\mid }$. For complex $z,w$, $\mathrm{\mid }zw\mathrm{\mid }=\mathrm{\mid }z\mathrm{\mid }\mathrm{\mid }w\mathrm{\mid }$, $\mathrm{\mid }z+w\mathrm{\mid }\le \mathrm{\mid }z\mathrm{\mid }+\mathrm{\mid }w\mathrm{\mid }$, and if $w\mathrm{\ne }0$ then $\mathrm{\mid }z\mathrm{/}w\mathrm{\mid }=\mathrm{\mid }z\mathrm{\mid }\mathrm{/}\mathrm{\mid }w\mathrm{\mid }$.

For every nonnegative exact dyadic input $x$ and $n\ge 1$, the exact dyadic pair returned by Exact dyadic bisection for square-root endpoints contains $\sqrt{x}$ and has width $2^{\lfloor j\mathrm{/}2\rfloor -n}$ when $x>0$, where $j$ is the exponent selected by that algorithm. The bisection loop terminates in at most $n+1$ iterations. Storing the endpoints after outward rounding preserves enclosure, although storage rounding may increase the width. For fixed $n$, both exact endpoint functions are nondecreasing in $x\ge 0$; directed storage rounding preserves that monotonicity.

Let $f$ be a real function that is $n+1$ times continuously differentiable on an open interval containing the finite real interval $X=[a,b]$. Here $n\ge 0$ is an integer, and continuous differentiability means that the derivatives through order $n+1$ exist there and are continuous. Let $c\in X$ be a dyadic point. Put

Supply dyadic intervals $C_j$ containing $f^{(j)}(c)\mathrm{/}j!$ for $0\le j\le n$ and a dyadic interval $D$ containing $f^{(n+1)}(X)$. Put $M=\max (\mathrm{\mid }\inf D\mathrm{\mid },\mathrm{\mid }\sup D\mathrm{\mid })$ and $R=\max (\mathrm{\mid }a-c\mathrm{\mid },\mathrm{\mid }b-c\mathrm{\mid })$. The centered Taylor-model enclosure is ${\mathcal{T}}_n(f,X,c,(C_j{)}_{j=0}^n,D)={\sum }_{j=0}^n{C}_j(X-c{)}^j+[-\frac{M{R}^{n+1}}{(n+1)!},\frac{M{R}^{n+1}}{(n+1)!}],$ where products, powers, sums, and the final sum are evaluated by outward-rounded interval operations. The first term carries input variation through the polynomial; the symmetric second interval is the explicit derivative-based remainder. If a required coefficient or derivative enclosure is unavailable, or any operation or endpoint conversion returns a non-finite status, no finite Taylor-model certificate is returned.

Let ${\sum }_{n\ge 0}{a}_n=A$ and ${\sum }_{n\ge 0}{b}_n=B$ be absolutely convergent real series. Define

Then ${\sum }_{n\ge 0}{c}_n$ converges absolutely and has sum $AB$. Thus the product of the two series can be grouped by total degree without changing its sum.

Every function differentiable at $c$ is continuous at $c$. If $f,g$ are real functions differentiable at $c$, then $f+g$ and $fg$ are differentiable there and

If $g(c)\mathrm{\ne }0$, then $f\mathrm{/}g$ is differentiable at $c$ and

If $f$ is differentiable at $c$, $g$ is differentiable at $f(c)$, and the composition is defined near $c$, then $g\circ f$ is differentiable at $c$ and $(g\circ f{)}^{\mathrm{\prime }}(c)=g^{\mathrm{\prime }}(f(c))f^{\mathrm{\prime }}(c)$. For every integer $m\ge 1$, $(x^m{)}^{\mathrm{\prime }}=m{x}^{m-1}$.

For every real $R\ge 0$, the factorial series ${\sum }_{n=0}^{\mathrm{\infty }}{R}^n\mathrm{/}n!$ converges. If $N\ge 0$ and $q=R\mathrm{/}(N+2)<1$, its tail after index $N$ satisfies ${\sum }_{n=N+1}^{\mathrm{\infty }}\frac{R^n}{n!}\le \frac{R^{N+1}}{(N+1)!(1-q)}\mathrm{.}$ Consequently the series defining $\exp x$, $\sin x$, $\cos x$ and each of their formal derivative series converge absolutely and uniformly on every bounded interval $\mathrm{\mid }x\mathrm{\mid }\le R$. For $0\le r<1$, the geometric series and the odd power series ${\sum }_{k\ge 0}{x}^{2k+1}\mathrm{/}(2k+1)$, with their formal derivative series, converge absolutely and uniformly on $\mathrm{\mid }x\mathrm{\mid }\le r$.