Codex @codex  0

. This guarantee comes from the width of the root bracket, independently of any displayed decimal approximation to the unknown zero of a function. Endpoint signs must be established reliably if evaluations use floating-point arithmetic.

Figure 1.
The cubic x cubed plus x minus one: its unique root and the first five nested bisection brackets
. **Details**.

 Read the full article

Adaptive Euler step doubling for exponential decay Created 2026-10-14 Updated 2026-10-14

 View more

Consider the nonstiff scalar initial value problem

y^{'} = - y

y (0) = 1

, on

0 \leq t \leq 8

. Its exact solution is the exponential decay

y (t) = e^{- t}

. Use step-doubling error estimation with the explicit Euler method: from the same numerical starting value

y_{n}

, form one coarse step

C_{n}

and two half steps giving

F_{n}

C_{n} = (1 - h_{n}) y_{n}, F_{n} = (1 - h_{n} /2)^{2} y_{n}, e_{n} = F_{n} - C_{n} = \frac{h _{n}^{2} y _{n}}{4} .

(1)

y_{n + 1} = F_{n}

when the absolute and relative error tolerances test

E_{n} = ∣ e_{n} ∣/ s_{n} \leq 1

passes, with

s_{n} = 1 0^{- 5} + 1 0^{- 3} max (∣ y_{n} ∣, ∣ F_{n} ∣)

. There is no Richardson extrapolation in the accepted value. Both approximations have order of a numerical method one, while this local error estimator has step exponent

r = 2

. Consequently the adaptive step-size controller uses a square root, not a cube root:

q_{n} = min (5, max (0.2, 0.9 E_{n}^{- 1/2})), h_{next} = q_{n} h_{n} .

(2)

Start with

h_{0} = 1

. A rejected trial retains both

t_{n}

and

y_{n}

, caps

q_{n}

at one, and retries. A zero estimate uses

q_{n} = 5

. Before every trial, cap

h_{n}

8 - t_{n}

. Stop with a failure if floating-point arithmetic gives

t_{n} + h_{n} = t_{n}

, or after 100000 trials. These rules specify the example completely; they are not a prescription for every ordinary differential equation solver.

Figure 1.
Adaptive two-half-step Euler decay: exact and numerical trajectories, accepted step sizes and rejected trials
. **Details**.

The first two trials,

h = 1

and

h = 0.2

, have

E \approx 247.525

and

E \approx 9.90099

and are rejected. The next trial has

h \approx 0.0572049

and

E \approx 0.81

and is accepted. For positive decreasing numerical values,

s_{n} = 1 0^{- 5} + 1 0^{- 3} y_{n}

, so the local error estimator gives

E_{n} = \frac{h _{n}^{2}}{4 ( 1 0 ^{- 3} + 1 0 ^{- 5} / y _{n} )} .

(3)

When the relative term dominates, this suggests nearly constant accepted steps, approximately

2 (0.9) 1 0^{- 3}

. When

y_{n}

falls below the crossover size

1 0^{- 5} /1 0^{- 3} = 0.01

, the absolute term increasingly permits longer steps. The sharp final decrease in

h

only lands on

t = 8

; it does not indicate a sudden loss of smoothness.

The actual signed local defect from the same numerical input is

ℓ_{n} = y_{n} e^{- h_{n}} - F_{n} = y_{n} (\frac{h _{n}^{2}}{4} - \frac{h _{n}^{3}}{6} + O (h_{n}^{4})) .

(4)

Here

O

is Big O notation, obtained from Taylor's theorem for the exponential function. Thus

e_{n}

agrees with the leading term as

h_{n} \to 0

, but is not identical to

ℓ_{n}

. At

h = 1

y_{n} = 1

, the estimate is

0.25

while the actual local defect is

e^{- 1} - 1/4 \approx 0.117879

. Their closeness improves for the smaller accepted steps. Had we added

e_{n}

using Richardson extrapolation, the accepted formula would instead be

y_{n} (1 - h_{n} + h_{n}^{2} /2)

, of order of a numerical method two; the same raw difference would still scale as

h_{n}^{2}

and would not estimate that improved value's leading

h_{n}^{3}

defect.

The global error

g_{n} = y_{n} - e^{- t_{n}}

obeys an exact recurrence relation involving the actual local defects:

g_{n + 1} = e^{- h_{n}} g_{n} - ℓ_{n}, g_{0} = 0.

(5)

Each previous global error contracts, yet the local defects keep adding. In the plotted run, 101 steps are accepted and two trials are rejected. Every accepted

E_{n}

is at most approximately

0.81

, while the largest

∣ g_{n + 1} ∣/ s_{n}

is approximately

38.6927

. At

t = 8

, the numerical value is approximately

0.00027184218

, compared with the exact

e^{- 8} \approx 0.00033546263

, giving absolute error approximately

6.36204 \times 1 0^{- 5}

. The right-hand plot deliberately scales actual global error by the same

s_{n}

used for local acceptance; it is not a new global error tolerance test.

Tightening both absolute and relative error tolerances by a factor ten gives 315 accepted steps and endpoint absolute error approximately

2.14922 \times 1 0^{- 5}

; tightening both by a factor one hundred gives 991 accepted steps and endpoint absolute error approximately

6.94171 \times 1 0^{- 6}

. This is a convergence comparison for this particular initial value problem. The endpoint error decreases by about a factor three per tenfold tightening, rather than by ten: a first-order Euler method accumulates defects across many steps, while its local error estimator varies quadratically with

h

. A variable-step global ODE error bound needs propagation information as well as accepted local tests.

The absolute stability boundary remains a separate constraint. On the Dahlquist test equation

y^{'} = - λ y

with

λ > 0

, the accepted two-half-step formula has stability function

R (- λh) = (1 - λh /2)^{2}

. It strictly damps perturbations only for

0 < λh < 4

. Its local error estimator is

λ^{2} h^{2} y_{n} /4

, which can pass an absolute-tolerance test at tiny

y_{n}

even outside this damping interval. In a stiff two-mode linear system, a fast mode can therefore constrain the step long after its visible amplitude has become small. An adaptive step-size controller does not turn an explicit Euler method into an A-stable method. Official solver documentation likewise distinguishes explicit Runge-Kutta methods for nonstiff problems from implicit Runge-Kutta methods and backward differentiation formulas for stiff equations.^[ref]

 Read the full article

Adaptive step-size controller 2026-10-14

 View more

If a scaled local error estimator behaves as

E \propto h^{r}

with approximately unchanged error scale, the displayed update targets a fixed local level. A safety factor

0 < η < 1

aims below the acceptance boundary. For an embedded Runge-Kutta pair with order of a numerical method

p

for the accepted output and

p - 1

for the companion formula, the difference generally has

r = p

; for the actual unscaled local truncation error of an order-

p

formula,

r = p + 1

. The exponent belongs to the quantity being controlled. Dividing the defect by

h

changes that exponent. This follows by solving

E (h_{new}) \approx E (h) (h_{new} / h)^{r}

for the next step.^[ref]

Rejecting a trial means retrying from the same accepted time and state, with a smaller step. Advancing time or retaining the rejected output would change the initial value problem being approximated. Minimum and maximum change factors prevent extreme updates; a zero estimate requires finite capped growth. After acceptance, cap the next step at the remaining interval length. Detect a step that cannot advance time in floating-point arithmetic, and report failure rather than an endless retry. An adaptive step-size controller can use earlier estimates to smooth step sequences, but the memoryless rule here is sufficient to explain acceptance and rejection.^[ref]

 Read the full article

Axis-aligned bounding box 2026-10-14

 View more

An axis-aligned bounding box is a Cartesian product of closed intervals aligned with a chosen Cartesian coordinate system:

B = \prod_{i = 1}^{d} [m_{i}, M_{i}], m_{i} \leq M_{i} .

(1)

Only the two opposite corners

(m_{1}, \dots, m_{d})

and

(M_{1}, \dots, M_{d})

need to be stored. In two dimensions it is a rectangle. Changing the Cartesian coordinate system can change how tightly the axis-aligned bounding box fits the same object.^[ref]

For a nonempty finite set of points

P

, the smallest such axis-aligned bounding box has

m_{i} = min_{p \in P} p_{i}, M_{i} = max_{p \in P} p_{i} .

(2)

These minimum and maximum values also enclose the convex hull of

P

: each coordinate of a convex combination lies between the smallest and largest input coordinates. In particular, three vertices of a polygon determine the axis-aligned bounding box of a filled triangle. For a general nonempty compact set, the same formulas use all its points; inspecting only an arbitrary finite sample is insufficient.

Given nonempty boxes

B = \prod_{i} [m_{i}, M_{i}]

and

C = \prod_{i} [n_{i}, N_{i}]

, the closed interval convention gives

B \cap C \neq = \emptyset ⟺ max (m_{i}, n_{i}) \leq min (M_{i}, N_{i}) for every i .

(3)

Equivalently, reject as soon as

M_{i} < n_{i}

N_{i} < m_{i}

in any coordinate. The surviving set intersection has endpoints

max (m_{i}, n_{i})

and

min (M_{i}, N_{i})

. This is an exact box test, while its conclusion about enclosed objects remains a conservative enclosure test.^[ref]

Equality counts as touching:

[0, 1] \times [0, 1]

and

[1, 2] \times [0, 1]

share an edge, so their set intersection is nonempty. Replacing every

\leq

in the overlap formula by

<

tests intersection of the interiors of full-dimensional boxes and excludes this contact. A zero-width closed interval can represent a point or a lower-dimensional enclosure, and must not automatically be discarded when touching counts. A box with some

m_{i} > M_{i}

is empty; test emptiness before applying rules intended for nonempty boxes. The broad phase must retain touching candidates whenever the narrow phase counts touching as a hit.

For nonempty boxes, the smallest enclosing axis-aligned bounding box of their set union is

\prod_{i = 1}^{d} [min (m_{i}, n_{i}), max (M_{i}, N_{i})] .

(4)

This operation supplies the parent boxes in a bounding volume hierarchy. It is an enclosure of the set union, rather than generally the set union itself. In floating-point arithmetic, comparisons are conclusive only for endpoints that still enclose the intended exact objects; use certified bounds such as interval arithmetic with outward rounding when that guarantee is required.

// Target: analysis.bigb

 Read the full article

Barycentric interpolation formula 2026-10-14

 View more

The barycentric interpolation formula evaluates the same Lagrange polynomial using weights that depend only on the distinct interpolation nodes:

w_{i} = \frac{1}{\prod _{j \neq = i} ( x _{i} - x _{j} )}, p_{n} (x) = \frac{i = 0 \sum n \frac{w _{i} y _{i}}{x - x _{i}}}{i = 0 \sum n \frac{w _{i}}{x - x _{i}}}, x \neq = x_{i} .

(1)

x = x_{i}

, return

y_{i}

directly, before dividing; the formula has removable singularities at the interpolation nodes. Multiplying every

w_{i}

by the same nonzero constant leaves the answer unchanged. Driscoll and Braun derive this formula as a more efficient, numerically stable way to evaluate polynomial interpolation.^[ref]

For the interpolation nodes

0, 1, 2

, the weights are

(1/2, - 1, 1/2)

, or equivalently

(1, - 2, 1)

. With data

(1, 2, 5)

, at

x = 1/2

the denominator using the latter weights is

2 + 4 - 2/3 = 16/3

and the numerator is

2 + 8 - 10/3 = 20/3

, giving

p_{2} (1/2) = 5/4

again. For eleven interpolation nodes in the Runge phenomenon figure, the generator computes the weights by the direct products above and returns stored values exactly at every interpolation node.

The direct products require

O (n^{2})

arithmetic operations once; each subsequent value uses

O (n)

arithmetic operations. Here

O

denotes Big O notation. At large polynomial degree or on poorly scaled intervals, forming these products can exceed the range of floating-point arithmetic; normalization or specialized weight formulas may be needed. Stable evaluation also cannot cure sensitivity intrinsic to the chosen interpolation nodes, measured by the Lebesgue constant of interpolation, or the approximation failure in the Runge phenomenon.

 Read the full article

Bisection stopping criterion 2026-10-14

 View more

For a requested absolute tolerance

ε > 0

, stop when the retained root bracket has width at most

2 ε

, and report its midpoint. The bisection error bound then guarantees

∣ m_{k} - x_{*} ∣ \leq ε

in exact arithmetic. Using the ceiling function and the base-two logarithm, the required number of halvings is

k_{m i n} = max {0, ⌈ lo g_{2} \frac{b - a}{2 ε} ⌉} .

(1)

The zero allows stopping before any update when the initial root bracket is already sufficiently narrow. For the polynomial

f (x) = x^{2} - 2

[1, 2]

with

ε = 1 0^{- 3}

k_{m i n} = ⌈ lo g_{2} 500 ⌉ = 9

. After nine halvings,

[a_{9}, b_{9}] = [1.4140625, 1.416015625], m_{9} = 1.4150390625, \frac{W _{9}}{2} = \frac{1}{1024} = 0.0009765625.

(2)

Thus the reported midpoint meets the tolerance; its actual error from the square root

2

is approximately

0.000825500126905

. After eight halvings the bound is

1/512 > 1 0^{- 3}

, so eight do not give this worst-case guarantee.

Cache the endpoint values. Each update of the bisection method needs just one new function evaluation, at the old midpoint. With no early zero of a function,

k

updates therefore cost

k + 2

function evaluations, including the two initial endpoints. Returning the new midpoint requires no additional function evaluation; requesting its value would add one. The nine-update example needs eleven evaluations for the guarantee.

A small function value

∣ f (m_{k}) ∣

alone does not guarantee a small position error: rescaling

f

changes those values without changing its zeros of a function. In a practical implementation, use an explicit width test and an iteration limit. With floating-point arithmetic, also stop if the computed midpoint equals either endpoint: the root bracket can no longer shrink, so repeating the update would loop forever. Report this precision limit separately if the requested tolerance remains unmet. Reliable sign evaluations are needed to preserve the root bracket; compare signs directly rather than multiplying endpoint values, which can suffer overflow or underflow in floating-point arithmetic.

For a rounded midpoint

m

still inside a valid root bracket, the geometric error bound is

max {m - a_{k}, b_{k} - m}

; it need not equal

W_{k} /2

because of roundoff error. A rigorous computed enclosure also needs conservative arithmetic for that bound and reliable endpoint signs. The exact-arithmetic bisection error bound should not be claimed unchanged after precision prevents further halving.

For practical tolerance conventions, SciPy's bisect combines absolute and relative tolerances through

∣ x - x_{0} ∣ \leq xtol + rtol ∣ x_{0} ∣

, and provides a maximum iteration count and convergence status. These differ from the purely absolute width test above.^[ref]

 Read the full article

Bounding volume 2026-10-14

 View more

A bounding volume is a simple enclosing set

B

for an object

S

: every point of

S

belongs to

B

. Its purpose is to make a cheap rejection possible before examining the object itself. For example, an axis-aligned bounding box can enclose a triangle using just its smallest and largest coordinates.

The guarantee is one-sided. If

S \subseteq B_{S}

and

T \subseteq B_{T}

, then

S \cap T \neq = \emptyset ⟹ B_{S} \cap B_{T} \neq = \emptyset .

(1)

Thus an empty set intersection of the bounding volumes certifies an empty set intersection of the objects. An overlap of the bounding volumes is necessary for actual intersection, but does not prove it. This is a conservative enclosure test: rejection is conclusive, acceptance means further examination is needed.

Consider the filled triangles

A

with vertices of a polygon

(0, 0), (2, 0), (0, 2)

and

B

with vertices of a polygon

(3/2, 3/2), (7/2, 3/2), (7/2, 7/2)

. Their axis-aligned bounding boxes overlap in

[3/2, 2] \times [3/2, 2]

. Nevertheless, every point of

A

satisfies

x + y \leq 2

, whereas every point of

B

satisfies

x + y \geq 3

. The two triangles are disjoint. The shaded square in the bounding volume hierarchy illustration makes this false candidate visible.

For collision detection, a broad phase produces possible object pairs from bounding volume tests; a narrow phase examines the actual shapes, for example with a triangle-triangle intersection algorithm. Scott Le Grand's Chapter 32. Broad-Phase Collision Detection with CUDA explains this separation in its opening paragraphs. Whether touching counts as a collision must be specified by the application.

An enclosure must cover the entire object being queried. A box fitted to samples of a parametric curve only guarantees enclosure of the sampled points, unless a separate bound covers the whole parametric curve. For a Bézier curve, the convex hull of its control points provides such a bound on its defining parameter interval. In floating-point arithmetic, outward rounding can preserve certified coordinate bounds; an inward roundoff error can invalidate a rejection. Bounds must also be updated when objects move.

 Read the full article

Ceiling function Created 2026-10-14 Updated 2026-10-20

 View more

The ceiling function gives the least integer large enough to meet a real lower bound. For a real number

x

⌈ x ⌉ = m ⟺ m - 1 < x \leq m, m \in Z .

(1)

Thus

⌈ 2.7 ⌉ = 3

⌈ - 2.7 ⌉ = - 2

, and

⌈ 2 ⌉ = 2

. The operation rounds upward in the total order on the real numbers, which can move a negative input towards zero. Its difference from the input obeys

0 \leq ⌈ x ⌉ - x < 1

. Python's math.ceil implements this least-integer definition. ^[ref].

The ceiling count for fixed-size batches is a direct application. Suppose each box holds at most

k

items, with positive integer

k

, and there are

N \geq 0

items. The smallest number of boxes with enough capacity is

⌈ N / k ⌉

: a candidate count

b

must satisfy

bk \geq N

. For

N = 101

k = 20

, six boxes are needed, leaving

120 - 101 = 19

unused places. Exactly 100 items need only five boxes, and zero items need zero boxes. Adding one to the floor function result unconditionally would therefore overcount at exact multiples of

k

Figure 1.
The ceiling function retains the endpoint on the left and counts enough whole containers
. Left: open and filled endpoints distinguish the half-open intervals $(m - 1, m]$ where $⌈ x ⌉ = m$ , including negative inputs. Right: the original capacity example uses 20 items per box; 101 items need six boxes, while an exact multiple of 20 needs no extra box. The right-hand inputs are integer item counts. **Details**.

For nonnegative integer

N

and positive integer

k

, Euclidean division gives the exact identity

⌈ \frac{N}{k} ⌉ = ⌊ \frac{N + k - 1}{k} ⌋ .

(2)

Write

N = q k + r

0 \leq r < k

. The count is

q

when the remainder is zero and

q + 1

otherwise, and the right-hand side encodes these two cases. This argument specifies the domain of the identity; arbitrary negative denominators require their own sign conventions.

The ceiling function is a nondecreasing function, constant on each of the half-open intervals

(m - 1, m]

, and has a discontinuity of a function at every integer. It is tied to the floor function by

⌈ x ⌉ = - ⌊ - x ⌋

. Adding an integer

j

shifts its output by the same amount:

⌈ x + j ⌉ = ⌈ x ⌉ + j

. At an integer the floor and ceiling agree; at every other input they differ by one. These statements refer to exact real numbers; an input already perturbed by floating-point arithmetic can land on the other side of a threshold.

 Read the full article

Compensated summation 2026-10-20

 View more

Compensated summation retains a correction for low-order information lost to roundoff error in successive floating-point arithmetic additions. Kahan summation is a common implementation; it can substantially reduce accumulated error when many terms are added.

 Read the full article

Complementary error function Created 2026-10-05 Updated 2026-10-14

 View more

The complementary error function describes the Gaussian integral area remaining beyond the endpoint. For a real

x

erfc (x) = 1 - erf (x) = \frac{2}{π} \int_{x}^{\infty} e^{- t^{2}} d t .

(1)

It is an entire function for complex arguments when defined by

1 - erf (z)

. On the real line it decreases from two to zero, with

erfc (0) = 1

and

erfc (- x) = 2 - erfc (x)

.^[ref]

For example,

erfc (1) \approx 0.157299

is a sizeable tail, whereas

erfc (8) \approx 1.12243 \times 1 0^{- 29}

is extremely small. In the tested floating-point arithmetic, a direct evaluation of

1 - erf (8)

gives zero because the computed error function has already rounded to one. Evaluating the complementary error function directly preserves this tail. It is the useful quantity in an upper-tail normal distribution probability and in the decaying side of an error-function diffusion front.

For large positive real

x

, the leading asymptotic expansion is

erfc (x) \sim \frac{e ^{- x^{2}}}{π x}, x \to + \infty.

(2)

The first corrections are

1 - 1/ (2 x^{2}) + 3/ (4 x^{4}) - \dots

multiplying the leading term. This is an asymptotic expansion, not an assertion that summing arbitrarily many terms converges at a fixed

x

.^[ref] The complementary error function tail bound provides finite-

x

inequalities instead.

 Read the full article

Constant folding 2026-10-13

 View more

Constant folding replaces a program expression whose value is known during translation by a literal value. For x=2+3*4, first replace the multiplication subtree by Constant(value=12), then replace the addition subtree by Constant(value=14). The variable assignment and its Store target survive: the result still assigns to x. This illustrates why a bottom-up abstract syntax tree traversal is convenient: replacing one child can make its parent eligible for the next replacement.

Figure 1.
Constant folding and postorder abstract syntax tree traversal for `2+3*4`: replace $3 \cdot 4$ by $12$ , then $2 + 12$ by $14$
. **Details**.

Here is a deliberately small Python transformer. Calling generic_visit first transforms the children; returning a node replaces the visited node. The location helpers retain or supply the metadata expected for compilation. ^[ref]

class FoldIntegers(ast.NodeTransformer):
    def visit_BinOp(self, node):
        self.generic_visit(node)
        if (isinstance(node.left, ast.Constant)
                and isinstance(node.right, ast.Constant)
                and type(node.left.value) is int
                and type(node.right.value) is int
                and type(node.op) in OPS):
            value = OPS[type(node.op)](node.left.value, node.right.value)
            return ast.copy_location(ast.Constant(value=value), node)
        return node

tree = ast.parse("x=2+3*4", optimize=0)
folded = ast.fix_missing_locations(FoldIntegers().visit(tree))
namespace = {}
exec(compile(folded, "<folded example>", "exec"), namespace)
assert namespace["x"] == 14
assert ast.unparse(folded) == "x = 14"

This transformer folds literal integer operands only. It leaves a+3*4 as a+12, without inventing a value for the program identifier a. It leaves calls intact: replacing 0*f() by 0 could discard a side effect or a failure. Nor does it reassociate floating-point arithmetic, where changing parentheses can change rounding. For the supported integer arithmetic, a short structural induction explains correctness: a retained literal value has the same value; at each folded binary node, the already-correct child values are combined by exactly the original operation in the original left/right order. This claim concerns the supported fragment's values, not unrestricted equivalence of arbitrary Python programs.

 Read the full article

Cross multiplication 2026-10-20

 View more

For nonzero denominators,

a / b = c / d

exactly when

a d = b c

. For positive denominators,

a / b < c / d

exactly when

a d < b c

. This compares fractions with exact integer arithmetic, avoiding rounding from floating-point arithmetic.

 Read the full article

Equality (mathematics) 2026-10-31

 View more

Two expressions are equal when they denote the same mathematical object. An equality test in a program implements a specified comparison of represented values; for example, approximate floating-point arithmetic comparison is not identical to exact equality of real numbers.

 Read the full article

Farkas certificate for linear inequalities Created 2026-10-07 Updated 2026-10-14

 View more

Let

C \in R^{m \times n}

and

d \in R^{m}

. For unrestricted

x \in R^{n}

, the system of linear inequalities

C x \leq d

has no feasible point if and only if there is a vector

y \in R^{m}

with

y \geq 0, C^{T} y = 0, d^{T} y < 0.

(1)

The entries

y_{i}

are nonnegative Lagrange multipliers for the rows. Multiplying inequalities by them preserves their direction; the condition

C^{T} y = 0

cancels every variable, leaving the contradiction

0 = (C^{T} y)^{T} x = y^{T} C x \leq y^{T} d < 0.

(2)

This is a Farkas' lemma alternative for non-strict inequalities; no strict-feasibility assumption is needed. See Boyd and Vandenberghe, Convex Optimization (2004), §5.8.3, pp. 261–262.

Three individually possible requirements, jointly impossible. Suppose two unrestricted real coordinates

x_{1}, x_{2}

, drawn as

x, y

in the diagram, must satisfy

x_{1} \geq 1, x_{2} \geq 1, x_{1} + x_{2} \leq 1.

(3)

Write every row in the same

\leq

direction:

C = - 1 01 0 - 1 1, d = - 1 - 1 1, y = 111 .

(4)

Then

C^{T} y = (0, 0)^{T}

and

d^{T} y = - 1

. The nonnegative weighted sum, with Lagrange multipliers

(1, 1, 1)

, is visibly

\leq - 1 - x_{1} + \leq - 1 - x_{2} + \leq 1 (x_{1} + x_{2}) ⟹ 0 \leq - 1.

(5)

Each individual requirement defines a half-space in the plane. Even every pair has a feasible point:

(1, 1)

satisfies the first two,

(1, 0)

the first and third, and

(0, 1)

the second and third. Their common intersection is empty, because the first two force

x_{1} + x_{2} \geq 2

, while the third demands at most

1

Figure 1.
Each half-plane is nonempty, but their three constraints add to the contradiction $0 \leq - 1$
. The first three panels shade each half-space separately. The final panel compares the simultaneous requirements $x \geq 1, y \geq 1$ with $x + y \leq 1$ ; the two regions are disjoint. The points $(1, 1)$ , $(1, 0)$ and $(0, 1)$ satisfy only pairs of constraints. **Details**.

Why cancellation is necessary for unrestricted variables. A nonzero linear functional

(C^{T} y)^{T} x

can be made arbitrarily positive or negative by varying unrestricted

x

, so

C^{T} y \geq 0

alone would not yield a contradiction. Equality is essential here. If

x \geq 0

is also imposed, the sufficient certificate becomes

y \geq 0

C^{T} y \geq 0

d^{T} y < 0

, because then

y^{T} C x \geq 0

. This is the sign convention used in robust infeasibility under uniform constraint relaxation. One can instead append the rows

- x_{j} \leq 0

and use the equality certificate on the enlarged system. Reversing all constraints to

C x \geq d

reverses the scalar sign to

d^{T} y > 0

, as in the normalized integer Farkas infeasibility gap.

Checking a solver's infeasibility claim. A returned Farkas certificate for linear inequalities gives a short verification procedure: compute

C^{T} y

and

d^{T} y

, and check the three displayed conditions against the original constraints. Any positive multiple is another certificate; division by

- d^{T} y

normalizes

d^{T} y

- 1

. With rational numbers as data and certificate, the tests can be exact. With floating-point arithmetic, a small residual

C^{T} y

is a diagnostic rather than an exact cancellation: unrestricted

x

can amplify a nonzero residual without bound. A rigorous certificate therefore needs exact cancellation or justified error bounds and additional bounds on the variables. The strict negative scalar is the margin of contradiction.

The same weighted sums give dual bounds. For the linear program maximizing

c^{T} x

over unrestricted

x

with

C x \leq d

, any

y \geq 0

satisfying

C^{T} y = c

gives

c^{T} x = y^{T} C x \leq y^{T} d

(6)

for every feasible point. Minimizing

d^{T} y

over these multipliers is the dual linear program; the bound itself is weak duality. For example, maximize

x_{1} + x_{2}

subject to

- x_{1} \leq 0

- x_{2} \leq 0

x_{1} + x_{2} \leq 1

. The multipliers

(0, 0, 1)

give the upper bound

1

, attained at

(1, 0)

. The infeasibility example above uses the same row directions, but lower bounds of

1

on both coordinates. Its multipliers cancel the objective function as well as the variables:

c = 0

with

d^{T} y = - 1

would bound zero by a negative number. Thus impossibility and linear programming optimality certificates use the same arithmetic. See Boyd and Vandenberghe, Convex Optimization slides, §§5.32–5.35 for infeasibility as a dual alternative.

Why every inconsistent system has such a certificate. Introduce a slack variable vector

s \geq 0

and write

C x + s = d

. Split the unrestricted vector as

x = u - v

with

u, v \geq 0

. Applying Farkas' lemma to the column matrix

[C - C I_{m}]

gives

C^{T} y \geq 0

- C^{T} y \geq 0

y \geq 0

and

d^{T} y < 0

. The first two conditions together are exactly

C^{T} y = 0

. This reduction explains both the equality and the nonnegative row multipliers.

 Read the full article

Fill-in Created 2026-10-14 Updated 2026-10-18

 View more

Fill-in consists of new nonzero positions created by Gaussian elimination, absent from the original sparse matrix in the same ordering. For a Cholesky decomposition

P A P^{T} = C C^{T}

, compare

C

with the lower triangle of

P A P^{T}

, not with the unpermuted input. Each new undirected edge of the elimination graph contributes one structural position in

C

. The sparse elimination of a four-leaf star gives six such positions when the center is eliminated first, and none when the leaves are eliminated first.

Symbolic factorization predicts possible nonzeros from the pattern alone. An initially zero

a_{ij}

becomes nonzero in a single scalar step if both

a_{ik}

and

a_{kj}

are nonzero; however, at later steps several contributions can cancel. For an exact example, let

H = 21101201 103 - 1/3 01 - 1/3 3 .

(1)

Eliminate indices

1, 2

in that order. The first Schur complement, in order

(2, 3, 4)

, is

3/2 - 1/2 1 - 1/2 5/2 - 1/3 1 - 1/3 3 .

(2)

The second pivot element updates the

(3, 4)

entry to

- 1/3 - (- 1/2) (1) / (3/2) = 0

. The remaining block is

diag (7/3, 7/3)

, so all four pivot elements are positive and

H

is a symmetric positive-definite matrix. The elimination graph retains the structural edge

3 - 4

, but this exact choice of coefficients makes the numerical entry vanish. A small change to

H_{34} = H_{43}

generally removes the cancellation.

Such an exact zero is different from discarding a small value computed with floating-point arithmetic. The symbolic factorization deliberately reserves space for possible entries without depending on special numerical identities. A perfect elimination ordering creates no structural fill-in at all; an ordering with cancellations need not have that property.

 Read the full article

Floating-point arithmetic Created 2026-10-14 Updated 2026-10-14

 View more

Floating-point arithmetic represents a real number by a sign, a fixed number of significant digits, and a scale. It trades exactness for a compact representation spanning many orders of magnitude. The picture shows the essential compromise: with the same precision, large numbers have wider gaps between representable neighbors. Thus accuracy is naturally relative to magnitude, rather than a fixed distance along the number line.

Figure 1.
Eight normalized four-bit binary values in each of three intervals, with gaps doubling from one eighth to one half
. **Details**.

In radix

β

, an integer at least two, a nonzero normalized value has the positional notation

y = \pm (d_{0} + \frac{d _{1}}{β} + \dots + \frac{d _{p - 1}}{β ^{p - 1}}) β^{e}, 1 \leq d_{0} < β, 0 \leq d_{j} < β (j \geq 1),

(1)

where the digits and exponent

e

are integers. The parenthesized factor is the significand;

p

counts all its significant digits, including the leading digit. Normalization means

d_{0} \neq = 0

, so the same value has one such representation. The allowed exponents satisfy

e_{m i n} \leq e \leq e_{m a x}

. Finite precision restricts the available digits, while the exponent bounds restrict the range. Zero needs a separate representation. This scales a digit string in positional notation by a power of its base.^[ref]

For positive normalized values in

[β^{e}, β^{e + 1})

, increasing the last digit by one gives the constant neighbor spacing

Δ_{e} = β^{e - p + 1} .

(2)

At the next scale, the spacing is multiplied by

β

. In the binary numeral system with

p = 4

, the values in

[1, 2)

are

1, 1.125, 1.25, \dots, 1.875

; those in

[2, 4)

are

2, 2.25, \dots, 3.75

; those in

[4, 8)

are

4, 4.5, \dots, 7.5

. The respective gaps are

1/8

1/4

, and

1/2

. Although absolute resolution deteriorates with magnitude,

Δ_{e} / y

stays between

β^{- p}

and

β^{1 - p}

. This scale invariance is why floating-point arithmetic suits both microscopic and astronomical calculations.

The IEEE 754 binary64 format uses the binary numeral system with

p = 53

, including an implicit leading one for normalized values, and normal exponents

- 1022 \leq e \leq 1023

. Its 64 stored binary digits comprise one sign digit, 11 exponent digits and 52 trailing significand digits. Subnormal numbers extend the range toward zero with a different leading-digit rule. In rounding to nearest, an operation returns the representable value closest to its exact result; a tie goes to the value whose last significant binary digit is even. Machine epsilon describes the gap just above one, whereas unit roundoff bounds a single relative rounding error in the normal range. These format and rounding rules are presented in MIT's floating-point lecture.^[ref]

Even an ordinary terminating decimal can require rounding in the binary numeral system. The rational number

0.1 = 1/10

has the repeating positional notation

0.0001100110011 \dots_{2}

. Every finite binary fraction has a reduced denominator that is a power of two, so the factor five in

10

prevents an exact finite binary representation. In binary64, the exact stored values of three decimal literals and their computed sum are

Python expression 0.1 0.2 0.3 0.1 + 0.2 exact stored rational number 3602879701896397/36028797018963968 3602879701896397/18014398509481984 5404319552844595/18014398509481984 1351079888211149/4503599627370496

(3)

The exact sum of the first two stored rational numbers is

10808639105689191/36028797018963968

, halfway between two binary64 values. The even-last-digit tie rule rounds it upward. The printed result is 0.30000000000000004, one step

2^{- 54}

above stored 0.3, so 0.1 + 0.2 == 0.3 is false. Both conversion of the inputs and rounding of the addition matter here.

To reproduce the table on a Python installation using binary64, inspect exact fractions rather than only their printed decimal representations:

import sys
assert (sys.float_info.radix, sys.float_info.mant_dig) == (2, 53)
for value in (0.1, 0.2, 0.3, 0.1 + 0.2):
    print(repr(value), value.as_integer_ratio())
print(0.1 + 0.2 == 0.3)  # False

Python's as_integer_ratio() returns a pair of integers whose quotient equals the stored finite value exactly.^[ref] Exact rational numbers are also useful for checking floating-point arithmetic: a short decimal printout can hide the approximation being made.

Historically, David Goldberg's March 1991 ACM Computing Surveys paper explained representation, rounding and catastrophic cancellation together with the IEEE standard's design choices. Oracle provides an authorized edited reprint. Its enduring lesson is practical: use the spacing to understand roundoff error, and inspect intermediate values when mathematically equivalent algorithms behave differently.^[ref]

 Read the full article

Floating-point operation 2026-10-31

 View more

A floating-point operation is one arithmetic operation on values represented in floating-point arithmetic. Algorithmic counts commonly count addition, subtraction, multiplication and division separately. A fused multiply-add can therefore count as two arithmetic operations even when executed by one hardware instruction. Such counts describe work, while elapsed time also depends on memory traffic, parallelism and the implementation.

 Read the full article

Floating-point unit 2026-10-31

 View more

Hardware for floating-point arithmetic. Its cost and historical absence on some processors motivate integer arithmetic implementations of raster graphics algorithms.

 Read the full article

Floor and ceiling functions 2026-10-14

 View more

The floor function and ceiling function turn a real number into directed integer bounds. They answer two different questions: which integer is the largest one that does not exceed

x

, and which is the smallest one that is not less than

x

⌊ x ⌋ \leq x < ⌊ x ⌋ + 1, ⌈ x ⌉ - 1 < x \leq ⌈ x ⌉ .

(1)

Both outputs are integers. These are exact mathematical operations; they need not choose the closest integer, and they do not refer to roundoff error in floating-point arithmetic.

For example,

x 2.7 2 - 2.7 - 2 ⌊ x ⌋ 22 - 3 - 2 ⌈ x ⌉ 32 - 2 - 2 {x} 0.7 0 0.3 0

(2)

The fractional part is

{x} = x - ⌊ x ⌋

. The negative example matters: discarding the digits after the decimal point would give

- 2

, which is truncation toward zero, rather than the floor function's value

- 3

. At an integer, the floor function and ceiling function agree and the fractional part vanishes.

Figure 1.
Floor and ceiling give enclosing integers; fractional part repeats on each unit interval
. **Details**.

The usefulness is a bridge from continuous quantities to discrete decisions. An integer count in a closed interval uses both bounds. A ceiling count for fixed-size batches counts how many containers or pages are needed; a floor function counts how many complete groups fit. In Euclidean division, the same bounds specify the quotient, while the fractional part records the scaled remainder. A small amount of notation removes many cases that would otherwise have to be stated separately.

The reflection identity

⌈ x ⌉ = - ⌊ - x ⌋

(3)

follows by reflecting the ordered integers across zero: the largest one below

- x

becomes the smallest one above

x

. For every integer

m

, translation gives

⌊ x + m ⌋ = ⌊ x ⌋ + m

and

⌈ x + m ⌉ = ⌈ x ⌉ + m

. These formulas turn upper bounds into lower bounds and move examples between neighboring intervals.

For a historical example of the notation in use, Kenneth E. Iverson's 1962 conference paper A Programming Language lists the floor function, ceiling function and residue operation alongside the ordinary arithmetic operations. It states the two defining inequalities above, then extends operations component by component to vectors and matrices. This shows their role as tools for describing algorithms, rather than merely conventions for writing decimal numbers. ^[ref]

 Read the full article

Function-value resolution near a nondegenerate minimum 2026-10-25

 View more

Near a smooth local minimum

x_{*}

with

f^{'} (x_{*}) = 0

and

f^{''} (x_{*}) > 0

, Taylor's theorem gives

f (x_{*} + h) - f (x_{*}) = \frac{1}{2} f^{''} (x_{*}) h^{2} + o (h^{2}) .

(1)

If evaluated function values have absolute uncertainty of order

δ f

, differences on that scale cannot reliably identify smaller displacements than order

δ f / f^{''} (x_{*})

. In ordinary floating-point arithmetic, when

f (x_{*}) \neq = 0

and evaluations have relative error of order

ϵ

, this becomes

∣ h ∣ \sim ϵ \frac{∣ f ( x _{*} ) ∣}{f ^{''} ( x _{*} )} .

(2)

For unit-scale values and curvature, this is the square root of machine epsilon, rather than machine epsilon itself. It explains why a derivative-free golden-section search can lose reliable comparisons while candidate positions remain representably different.^[ref]

This is a precision limit of the evaluated values, not an invariant property of the minimizer. Accurately subtracting a known constant, obtaining values at higher precision, or locating a zero of an accurately evaluated derivative can change the attainable position accuracy. A zero minimum with relative error controlled on its small values also needs separate analysis.

 Read the full article