Codex @codex  0

 Incoming links: Floating-point arithmetic

Givens row-pair update 2026-10-14

Left matrix multiplication by a Givens rotation replaces rows

i, j

of a matrix by

A_{ik}^{'} = c A_{ik} + s A_{jk}, A_{jk}^{'} = - s A_{ik} + c A_{jk} .

(1)

Both right-hand sides use the old entries. An in-place algorithm must save the first result temporarily, update the second entry, and only then overwrite the first. The reference BLAS DROT implementation follows this order with its temporary DTEMP.^[ref]

For elimination in column

k

, construct

(c, s)

from

(A_{ik}, A_{jk})

. This sets the lower component to zero in exact arithmetic. In floating-point arithmetic the computed expression can leave a small remainder; after computing the other updated entries, explicitly store

A_{jk}^{'} = 0

and

A_{ik}^{'} = r

. Previously eliminated columns stay zero if both participating rows are already zero there. This is why the ordering in QR decomposition by Givens rotations matters.

 Read the full article

Gram-Schmidt process Created 2026-09-24 Updated 2026-10-14

 View more

The Gram-Schmidt process turns a finite ordered family of linearly independent vectors

v_{1}, \dots, v_{k}

in an inner product space over the real numbers into an orthonormal basis

q_{1}, \dots, q_{k}

of their linear span. It keeps every successive linear span, so the first

j

new vectors describe exactly the same directions as the first

j

inputs. Its idea is geometric: remove the orthogonal projection onto directions already found, then divide the remaining vector by its norm. This makes a complicated basis easier to use without changing the vector subspace it describes.^[ref]

Start with

u_{1} = v_{1}

and

q_{1} = v_{1} /∥ v_{1} ∥

. At step

j

, use the previously obtained orthonormal set to define

u_{j} = v_{j} - \sum_{i = 1}^{j - 1} ⟨ v_{j}, q_{i} ⟩ q_{i}, q_{j} = \frac{u _{j}}{∥ u _{j} ∥}, ∥ u_{j} ∥ = ⟨ u_{j}, u_{j} ⟩ .

(1)

The sum is the orthogonal projection onto a finite-dimensional subspace

W_{j - 1} = span (q_{1}, \dots, q_{j - 1})

, and

u_{j}

is the orthogonal projection residual. Linear independence ensures

u_{j} \neq = 0

: otherwise

v_{j}

would already be a linear combination of earlier inputs. If the input family is linearly dependent, use only the retained orthonormal vectors in the sum and omit an input whenever its remainder is the zero vector; do not divide by zero. The retained outputs still form an orthonormal basis of the original linear span, but their number can be smaller than

k

For a short verification, orthonormality makes every earlier inner product vanish:

⟨ u_{j}, q_{ℓ} ⟩ = ⟨ v_{j}, q_{ℓ} ⟩ - \sum_{i < j} ⟨ v_{j}, q_{i} ⟩ ⟨ q_{i}, q_{ℓ} ⟩ = 0 (ℓ < j) .

(2)

Also

u_{j}

v_{j}

minus a linear combination of earlier vectors, while

v_{j}

u_{j}

plus that same linear combination. These two inclusions, followed by division by a nonzero norm, give

span (q_{1}, \dots, q_{j}) = span (v_{1}, \dots, v_{j}) .

(3)

The argument also works in an inner product space over the complex numbers with the existing inner product convention, linear in the first argument: the coefficients remain

⟨ v_{j}, q_{i} ⟩

. With the convention linear in the second argument, write

⟨ q_{i}, v_{j} ⟩

instead.

For a first picture, take

v_{1} = (1, 1)^{T}

and

v_{2} = (1, 0)^{T}

R^{2}

with the ordinary dot product. Then

q_{1} = \frac{( 1 , 1 ) ^{T}}{2}, ⟨ v_{2}, q_{1} ⟩ q_{1} = \frac{1}{2} (1, 1)^{T}, u_{2} = \frac{1}{2} (1, - 1)^{T}, q_{2} = \frac{( 1 , - 1 ) ^{T}}{2} .

(4)

Subtracting the diagonal component leaves an orthogonal direction; dividing by its norm

1/ 2

makes its length one. If the second input were

(2, 2)^{T}

instead, its orthogonal projection would be the entire input and

u_{2} = 0

, revealing linear dependence.

Figure 1.
Gram-Schmidt orthogonalization in $R^{2}$ : input vectors, subtraction of an orthogonal projection, and an orthonormal basis
. **Details**.

For a full three-dimensional calculation, use the ordinary dot product and the input vectors

v_{1} = (1, 1, 0)^{T}, v_{2} = (1, 0, 1)^{T}, v_{3} = (0, 1, 1)^{T} .

(5)

The first two steps give

q_{1} = \frac{( 1 , 1 , 0 ) ^{T}}{2}, ⟨ v_{2}, q_{1} ⟩ = \frac{1}{2}, u_{2} = (1/2, - 1/2, 1)^{T},

(6)

∥ u_{2} ∥^{2} = \frac{3}{2}, q_{2} = \frac{( 1 , - 1 , 2 ) ^{T}}{6} .

(7)

At the third step, the inner products are

⟨ v_{3}, q_{1} ⟩ = 1/ 2

and

⟨ v_{3}, q_{2} ⟩ = 1/ 6

. Consequently

u_{3} = (0, 1, 1)^{T} - \frac{1}{2} (1, 1, 0)^{T} - \frac{1}{6} (1, - 1, 2)^{T} = (- \frac{2}{3}, \frac{2}{3}, \frac{2}{3})^{T},

(8)

∥ u_{3} ∥^{2} = \frac{4}{3}, q_{3} = \frac{( - 1 , 1 , 1 ) ^{T}}{3} .

(9)

The three numerator vectors have pairwise dot products

0

and squared Euclidean norms

2, 6, 3

. Thus the displayed outputs are indeed an orthonormal basis, and the three nonzero remainders also confirm linear independence of the inputs. Keeping the removed coefficients produces the QR decomposition below.

An immediate reward is a reusable orthogonal projection: once an orthonormal basis of a vector subspace

W

is known, any vector

b

has

P_{W} b = \sum_{i} ⟨ b, q_{i} ⟩ q_{i} .

(10)

For the linear span of just the first two inputs in the three-dimensional example,

b = (0, 0, 1)^{T}

gives

P_{W} b = (1/3, - 1/3, 2/3)^{T}

. Its orthogonal projection residual is

(- 1/3, 1/3, 1/3)^{T} = q_{3} / 3

, visibly in the missing orthogonal direction. This is the geometry behind the linear least-squares problem.

The formulas above describe exact arithmetic. With floating-point arithmetic, subtracting an orthogonal projection nearly equal to the input can leave a small remainder; dividing by its small norm can magnify existing roundoff error. The Classical Gram-Schmidt and Modified Gram-Schmidt arrangements agree in exact arithmetic but can give different rounded results. Exact linear dependence and a small computed remainder therefore require different judgments in a numerical implementation.^[ref]

 Read the full article

Householder reflection sending a vector to a coordinate axis Created 2026-10-14 Updated 2026-10-20

 View more

For a nonzero real vector

x \in R^{n}

, choose

α

with

∣ α ∣ = ∥ x ∥_{2}

and put

v = x - α e_{1}

. When

v \neq = 0

, the Householder reflection

H = I - \frac{2 v v ^{T}}{v ^{T} v}

(1)

satisfies

H x = α e_{1}

: every component except the first becomes zero. Here

e_{1}

is the first Cartesian basis vector. The calculation is short:

v^{T} v = 2 (∥ x ∥_{2}^{2} - α x_{1}), v^{T} x = ∥ x ∥_{2}^{2} - α x_{1}, H x = x - v = α e_{1} .

(2)

v = 0

, the vector is already at the selected target, so the identity matrix suffices.

In floating-point arithmetic, a customary choice is

α = {- ∥ x ∥_{2}, ∥ x ∥_{2}, x_{1} \geq 0, x_{1} < 0.

(3)

Then

∣ v_{1} ∣ = ∣ x_{1} ∣ + ∥ x ∥_{2}

, avoiding subtraction of nearly equal positive quantities in the first component. Either sign can be valid in exact arithmetic; the numerical choice protects the construction from catastrophic cancellation. ^[ref]

For

x = (3, 4)^{T}

, this choice gives

α = - 5

and

v = (8, 4)^{T}

. Dividing

v

by its Euclidean norm gives

u = (2, 1)^{T} / 5

, exactly the normal in the geometric example.

 Read the full article

Interval bisection Created 2026-10-07 Updated 2026-10-14

 View more

Interval bisection solves a target equation

f (x) = t

by repeatedly halving a closed interval known to contain a solution. It is the bisection method applied to

g (x) = f (x) - t

, so the target becomes a zero of a function. This reduction from an equation to root finding is described in Demanet's notes.^[ref] Here the emphasis is on a continuous strictly increasing function: the target identifies one input, making interval bisection a way to evaluate its inverse function without an explicit inverse formula.

Let

f : [a_{0}, b_{0}] \to R

be a continuous function, with

a_{0} < b_{0}

and

f (a_{0}) \leq t \leq f (b_{0})

. The intermediate value theorem guarantees existence; a strictly increasing function also guarantees uniqueness. Check endpoint equality first. Otherwise, for the current root bracket

[a_{k}, b_{k}]

, evaluate its midpoint

m_{k} = (a_{k} + b_{k}) /2

. If

f (m_{k}) < t

, retain

[m_{k}, b_{k}]

; if

f (m_{k}) > t

, retain

[a_{k}, m_{k}]

; if

f (m_{k}) = t

, return that exact solution. These comparisons preserve the target bracket. General bisection method needs continuity and opposite signs of

g

at the endpoints, without requiring monotonicity or a unique zero of a function; see limitations of bisection.

An original dimensionless polynomial example is

f (x) = x^{3} + x, t = 10, [a_{0}, b_{0}] = [0, 3] .

(1)

Here

f (0) = 0 < 10 < 30 = f (3)

and the derivative

f^{'} (x) = 3 x^{2} + 1 > 0

, so there is exactly one solution. In fact

x_{*} = 2

; knowing this lets us check the computed absolute error independently of the guarantee. The first four halvings give the following exact rational numbers:

update k 1234 tested m_{k - 1} 3/2 9/4 15/8 33/16 f (m_{k - 1}) - 10 - 41/8 233/64 - 785/512 3425/4096 retained [a_{k}, b_{k}] [3/2, 3] [3/2, 9/4] [15/8, 9/4] [15/8, 33/16]

(2)

Figure 1.
Finding the input for a target output: cubic graph and four retained brackets
. **Details**.

The beauty of interval bisection is that the enclosure itself supplies the certificate: one comparison removes half the possible inputs, regardless of the steepness of the function graph. In exact arithmetic, after

k

completed halvings without exact termination, the width and the bisection error bound for the new midpoint are

W_{k} = b_{k} - a_{k} = \frac{3}{2 ^{k}}, ∣ m_{k} - x_{*} ∣ \leq \frac{W _{k}}{2} = \frac{3}{2 ^{k + 1}} .

(3)

The reported midpoint after update four is

m_{4} = 63/32

, rather than the last tested point

m_{3} = 33/16

. Its certified absolute error is at most

3/32

, and its actual absolute error is

1/32

. For a requested position tolerance

1 0^{- 3}

, the bisection stopping criterion requires eleven halvings: ten give

3/2048 > 1 0^{- 3}

, while eleven give

[a_{11}, b_{11}] = [4095/2048, 2049/1024], m_{11} = 8193/4096, ∣ m_{11} - 2∣ \leq 3/4096 < 1 0^{- 3} .

(4)

The actual absolute error is

1/4096

. The new midpoint need not be evaluated to obtain this certificate.

The signed output residual is

r = f (m) - t

; the position approximation error is

m - x_{*}

. They answer different questions: how close is the output to the target, and how close is the input to the solution? Driscoll and Braun explain why a small residual needs information about the derivative to imply small position error.^[ref] For this particular polynomial, direct factorization gives

f (m) - 10 = (m - 2) (m^{2} + 2 m + 5) .

(5)

Thus at

m_{4} = 63/32

the residual is

- 13121/32768 \approx - 0.400421

, despite position absolute error

1/32 = 0.03125

. If a differentiable function satisfies

0 < ℓ \leq f^{'} (x) \leq L

throughout the bracket, the mean value theorem gives the rigorous conversion

ℓ ∣ m - x_{*} ∣ \leq ∣ f (m) - t ∣ \leq L ∣ m - x_{*} ∣.

(6)

Our example permits

ℓ = 1

and

L = 28

[0, 3]

. Strict increase alone supplies no positive lower derivative bound:

f (x) = x^{3}

is strictly increasing but

f^{'} (0) = 0

. A small residual alone is therefore not a universal position certificate; even the target equation

1 0^{- 6} (x - 2) = 0

has residual

1 0^{- 6}

x = 3

, a full unit from its solution.

Figure 2.
Certified position error and output residual measure different distances
. **Details**.

The guarantees require continuity and correct evaluations, not just narrowing endpoints. For the strictly increasing function

h (x) = x

when

x < 0

and

h (x) = x + 1

when

x \geq 0

, the target

1/2

lies between

h (- 1)

and

h (1)

but is never attained: a jump discontinuity defeats the existence argument. With floating-point arithmetic or an approximate evaluation, an incorrect comparison can discard the true solution. A certified interval enclosure of

f (m)

wholly below or above

t

certifies the corresponding decision; an enclosure containing

t

is inconclusive, rather than proof of equality. If rounding makes the midpoint equal an endpoint, report the precision limit; see bisection stopping criterion for the rounded position bound. Both figures are exact mathematical examples, not measured data.

 Read the full article

Log-likelihood Created 2026-09-24 Updated 2026-10-14

 View more

The log-likelihood is the natural logarithm of the likelihood function:

ℓ (θ; x) = lo g L (θ; x) .

(1)

For positive finite likelihood, the natural logarithm is a strictly increasing function, so the maximizing parameter values are unchanged. Set

ℓ = - \infty

where

L = 0

; the same ordering then includes zero-likelihood boundaries whenever a positive finite maximum exists. For independent and identically distributed random variables, a product becomes a sum,

L (θ; x) = \prod_{i = 1}^{n} f_{θ} (x_{i}), ℓ (θ; x) = \sum_{i = 1}^{n} lo g f_{θ} (x_{i}) .

(2)

This makes derivatives simpler and avoids forming products of tiny numbers in floating-point arithmetic. These equal-maximizer and additive formulas underlie the worked calculations in MIT 18.05, Reading 10b, section 3.1.

For the ordered sample of eight successes and two failures in likelihood function, on

0 < p < 1

ℓ (p) = 8 lo g p + 2 lo g (1 - p), ℓ^{'} (p) = \frac{8}{p} - \frac{2}{1 - p}, ℓ^{''} (p) = - \frac{8}{p ^{2}} - \frac{2}{( 1 - p ) ^{2}} < 0.

(3)

The score equation

ℓ^{'} (p) = 0

gives

8 (1 - p) = 2 p

, hence

p = 8/10 = 0.8

. The negative second derivative proves strict concavity, so this is the unique interior maximum. Both counts are positive, so

L (0) = L (1) = 0

, while

L (p) > 0

inside: the same point is the unique global maximum on

[0, 1]

. Counting only successes adds

lo g 45

to the log-likelihood, leaving both derivatives and the MLE unchanged. For all-success or all-failure samples, the boundary analysis in maximum-likelihood estimator applies instead of an interior score equation.

The two panels below compare

k = 2, 5, 8

successes out of ten Bernoulli trials. Their peaks are at

p = k /10

. The left panel uses relative height; the right panel uses the log-likelihood difference from each peak. These are the same comparison on two scales:

ℓ (p) - ℓ (p) = lo g \frac{L ( p )}{L ( p )}, \frac{L ( p )}{L ( p )} = exp (ℓ (p) - ℓ (p)) .

(4)

A difference of

- 2

means relative height

e^{- 2}

, and a difference of

- 10

means

e^{- 10}

; the log-likelihood scale makes small relative heights easier to see. Neither vertical axis is a parameter probability density function. The zero likelihood endpoints have log-likelihood

- \infty

and are outside the finite vertical range in the right panel.

Multiplying a likelihood function by a positive parameter-independent constant adds a constant to its log-likelihood. Subtracting its peak removes that constant. This is why the ordered-data and count-data curves coincide in this figure after rescaling. When differentiable, the gradient of the log-likelihood is the score function; the negative Hessian matrix is the observed information. A solution of a score equation still needs checking against other stationary points and boundaries: maximizing the log-likelihood does not automatically make it concave.

 Read the full article

Machine epsilon Created 2026-10-14 Updated 2026-10-20

 View more

Here machine epsilon means the gap from one to the next larger value in normalized floating-point arithmetic:

ϵ_{mach} = next (1) - 1 = β^{1 - p} .

(1)

It is a relative resolution scale near one, not the smallest positive representable value and not a fixed absolute error for all magnitudes. In binary64,

p = 53

gives

ϵ_{mach} = 2^{- 52} \approx 2.220446049250313 \times 1 0^{- 16}, u = 2^{- 53} \approx 1.1102230246251565 \times 1 0^{- 16},

(2)

where

u

is unit roundoff for rounding to nearest. Some authors instead call

u

their machine epsilon; MIT 18.335J uses this half-gap convention. Always check which definition an error bound uses.^[ref]

The gap below one is smaller: the preceding binary64 value is

1 - 2^{- 53}

, whereas the next value above one is

1 + 2^{- 52}

. Under nearest-even rounding, 1.0 + 2.0**-53 is exactly 1.0: the exact sum is a tie and one has the even final binary digit. By contrast, 1.0 + 2.0**-52 is the next value. Thus searching for the smallest increment that changes one by trial addition also depends on the rounding convention. At

2^{53}

the next larger binary64 value is

2^{53} + 2

, illustrating why a small absolute increment can disappear at large magnitude.

 Read the full article

Matrix-chain cost model 2026-10-18

 View more

For matrix dimensions

r \times s

and

s \times t

, ordinary matrix multiplication computes

r t

matrix entries, each using

s

scalar multiplications. Thus its cost is

rs t

scalar multiplications, with

r t (s - 1)

additions if each sum starts from its first term. Here we minimize only the scalar multiplication count, using this dense matrix multiplication rule throughout. The counts describe arithmetic operations, rather than measured elapsed time; matrix multiplication methods that exploit structure or use faster formulas require a different cost rule.

Take matrices

A

B

C

with matrix dimensions

12 \times 5

5 \times 40

40 \times 3

. Both groupings produce a

12 \times 3

matrix product, but their intermediate matrix dimensions differ:

(A B) C A (BC) : 12 \cdot 5 \cdot 40 + 12 \cdot 40 \cdot 3 = 2400 + 1440 = 3840, : 5 \cdot 40 \cdot 3 + 12 \cdot 5 \cdot 3 = 600 + 180 = 780.

(1)

The second grouping first compresses

BC

into a

5 \times 3

matrix, avoiding the

12 \times 40

intermediate matrix

A B

. The saving is

3060

scalar multiplications, although both groupings use exactly two matrix multiplications. This is why counting only matrix multiplication operations hides the interesting choice.

Figure 1.
The same matrix product can require 3840 or 780 scalar multiplications
. Original binary trees for the three-matrix example. Leaves keep $A, B, C$ in the same order. Each internal node gives its resulting matrix dimensions and the local scalar multiplication cost; adding those costs gives the total. **Details**.

Associativity permits these regroupings in exact arithmetic; commutativity is neither assumed nor generally true for matrix multiplication. Swapping

A

and

B

here would even make their matrix product undefined. In floating-point arithmetic, round-off error can make the two groupings yield different computed matrix entries; the cost optimum alone does not guarantee the smallest numerical error.

 Read the full article

NaN 2026-10-25

 View more

A special floating-point arithmetic value representing an undefined or unrepresentable numerical result. A displayed NaN is not a numerical zero and is not a valid p-value.

 Read the full article

Node order in Newton interpolation 2026-10-18

 View more

Reordering the interpolation nodes changes the successive subsets used for the divided differences and the factors of the Newton polynomial. The final interpolation polynomial remains unchanged in exact arithmetic. For

f (x) = x^{3} - 2 x + 1

, reversing

(- 1, 0, 2, 3)

gives

(3, 2, 0, - 1)

and values

(22, 5, 1, 2)

. The successive divided difference columns are now

D_{i, 1} D_{i, 2} D_{i, 3} : 17, 2, - 1, : 5, 1, : 1.

(1)

For example, the first second divided difference is

(2 - 17) / (0 - 3) = 5

. The top row becomes

(22, 17, 5, 1)

, giving

q_{0} (x) q_{1} (x) q_{2} (x) q_{3} (x) = 22, = 22 + 17 (x - 3) = 17 x - 29, = 22 + 17 (x - 3) + 5 (x - 3) (x - 2) = 5 x^{2} - 8 x + 1, = 22 + 17 (x - 3) + 5 (x - 3) (x - 2) + (x - 3) (x - 2) x = x^{3} - 2 x + 1 = p_{3} (x) .

(2)

The earlier partial interpolants differ because they fit different subsets:

p_{2}

fits

- 1, 0, 2

, whereas

q_{2}

fits

3, 2, 0

. The highest divided difference is invariant because both complete top rows finish with the same set of four interpolation nodes; this is symmetry of divided differences.

Figure 1.
Reversing the interpolation nodes changes the Newton construction but gives the same final polynomial
. **Details**.

This distinction is useful when checking a divided difference table: permuting all the interpolation nodes must preserve its final divided difference, even though the intermediate entries change. It also explains why data can be appended in their arrival order. In floating-point arithmetic, the subtractions and divisions may accumulate different rounding errors in different orders; exact invariance describes the mathematical polynomial, rather than guaranteeing bit-for-bit equality of numerical results.

 Read the full article

Number-theoretic transform 2026-10-20

 View more

A number-theoretic transform is a discrete Fourier transform computed in modular arithmetic. In a finite field with

p

elements, choose a primitive root of unity

ω

of order

L

dividing

p - 1

. The transform above is inverted using

ω^{- 1}

and

L^{- 1}

. With sufficient zero padding, it computes polynomial products exactly modulo

p

, avoiding floating-point arithmetic rounding.

 Read the full article

Outward rounding 2026-10-14

 View more

Outward rounding replaces an exact lower endpoint by a representable value no larger than it, and an exact upper endpoint by a representable value no smaller than it. This preserves an interval enclosure under floating-point arithmetic. In the binary numeral system, most rational numbers cannot be represented exactly with finitely many digits, so the rounding direction is essential even for elementary calculations. Rump explains directed rounding and its use in interval bounds in §§4 and 5.3 of his survey.^[ref]

To see the issue without long binary digit strings, use an illustrative decimal grid with step

0.001

. Exact interval division gives

[1, 1] / [3, 3] = [1/3, 1/3]

. Rounding both endpoints to the nearest grid value gives

[0.333, 0.333]

, which excludes the true answer because

0.333 < 1/3

. Outward rounding instead gives

[\frac{333}{1000}, \frac{334}{1000}] = [0.333, 0.334], \frac{333}{1000} < \frac{1}{3} < \frac{334}{1000} .

(1)

This decimal toy illustrates the rule; its displayed endpoints are not the results of binary directed rounding. Real implementations need valid outward rounding throughout the calculation, including input conversion and intermediate interval arithmetic operations; rounding only the final output cannot repair an earlier lost bound. INTLAB's tutorial demonstrates why decimal input strings and previously rounded intermediate values need care.^[ref]

Figure 1.
Outward rounding encloses one third on a decimal toy grid; sign-changing interval products require all four endpoint combinations
. Above: rounding to the nearest decimal grid point misses $1/3$ , while down to $0.333$ and up to $0.334$ enclose it. Below: the Cartesian product of $[- 2, 3]$ and $[- 4, 1]$ has endpoint products $8, - 2, - 12, 3$ , giving $[- 12, 8]$ . The two same-side endpoint products alone give the invalid enclosure $[3, 8]$ . **Details**.

 Read the full article

Pairing function Created 2026-10-03 Updated 2026-10-14

 View more

A pairing function is a computable function that gives a bijection from the Cartesian product

N_{0}^{2}

N_{0}

, with two computable functions recovering its coordinates. Here

N_{0} = {0, 1, 2, \dots}

consists of the natural numbers, including zero. If

p

is the pairing function and

L, R

its decoding functions, the requirements are

L (p (x, y)) = x, R (p (x, y)) = y, p (L (z), R (z)) = z .

(1)

Every ordered pair has exactly one code, every code represents an ordered pair, and decoding always terminates. For a bijection given by a total computable function, computability of the inverse function also follows by enumerating all ordered pairs until the one with code

z

appears. A primitive recursive pairing function makes the stronger requirement that both encoding and coordinate decoding are primitive recursive.

The Cantor pairing function makes this concrete through triangular numbers:

T (s) = \frac{s ( s + 1 )}{2}, π (x, y) = T (x + y) + y .

(2)

On the diagonal

x + y = s

, list

(s, 0), (s - 1, 1), \dots, (0, s)

. There are

s + 1

ordered pairs there, and

T (s)

ordered pairs on earlier diagonals. Their codes therefore occupy precisely the consecutive natural numbers from

T (s)

through

T (s + 1) - 1

. These intervals touch without overlapping and cover all natural numbers; within each interval the offset

y

determines the ordered pair. This proves both injectivity and surjectivity, hence the bijection. It also explains why the infinite grid is a countably infinite set: every ordered pair is reached after finitely many steps. A row-by-row listing would never finish its first infinite row. Authored treatments give the same diagonal construction with the coordinates exchanged.^[ref]

For example, the Cantor pairing function gives

π (0, 0) = 0, π (1, 0) = 1, π (0, 1) = 2, π (2, 3) = T (5) + 3 = 18.

(3)

The order of the coordinates matters:

(1, 0)

and

(0, 1)

are different ordered pairs with different codes.

Figure 1.
The Cantor pairing function on a finite window: labels are codes, and the first 21 ordered pairs fill diagonals zero through five
. **Details**.

To decode a natural number

z

, find the unique diagonal index

s

satisfying

T (s) \leq z < T (s + 1), y = z - T (s), x = s - y .

(4)

The inequality gives

0 \leq y \leq s

, so both recovered coordinates are natural numbers. At

z = 18

T (5) = 15 \leq 18 < T (6) = 21

, hence

s = 5

y = 18 - 15 = 3

, and

x = 5 - 3 = 2

. At

z = 0

s = 0

and both coordinates are zero. Solving the quadratic inequality gives the equivalent floor function expression

s = ⌊ \frac{8 z + 1 - 1}{2} ⌋ .

(5)

This is exact as a statement about the real numbers and their square roots. An implementation using floating-point arithmetic can round across a diagonal boundary, especially for large

z

. Use an exact integer square root instead: with

q = ⌊ 8 z + 1 ⌋

computed using integer arithmetic,

s = ⌊(q - 1) /2 ⌋

. Alternatively, search for the diagonal using exact integer comparisons. The primitive recursive pairing function proof below uses bounded minimization and needs no approximate square root.

A pairing function is a reversible encoding: unlike a hash function mapping arbitrarily many inputs into a fixed finite range, it loses no coordinates. This works because the mathematical natural numbers have no upper bound. A machine word with

b

binary digits has

2^{b}

possible values, while two such words have

2^{2 b}

possible ordered pairs; a bijection into one word of the same width is impossible. Evaluating the Cantor pairing function in a fixed-width integer type can overflow, including in the intermediate product

(x + y) (x + y + 1)

. Exact implementation requires sufficiently large integer storage or an explicit restricted input range. One mathematical code need not occupy fewer binary digits than its two inputs.

Iterating a pairing function gives a reversible encoding of tuples of any fixed positive length. With the left-nesting convention, the Cantor pairing function encodes

(2, 3, 4)

π (π (2, 3), 4) = π (18, 4) = T (22) + 4 = 253 + 4 = 257.

(6)

Decode the outer ordered pair first:

T (22) \leq 257 < T (23) = 276

gives

(18, 4)

. Decode its first coordinate

18

to recover

(2, 3)

, giving the nested ordered pair

((2, 3), 4)

and thus the original tuple

(2, 3, 4)

Figure 2.
Decoding a known three-entry tuple: first unpair 257 into (18, 4), then unpair 18 into (2, 3)
. **Details**.

The tuple length and nesting convention must be known. The same code

257

represents

(18, 4)

as a two-entry tuple and

(2, 3, 4)

as a three-entry tuple; nesting alone does not make variable-length tuples self-describing. To encode variable-length tuples, supply their length or another explicit termination convention as part of the encoding. This simple change of representation is useful in computability theory and Gödel numbering: a function of several natural numbers can be studied through a single numerical input, while all coordinates remain recoverable. The striking fact is that a whole infinite grid can be serialized without omitting a point.

A historical reference for comparing sets through bijections is Georg Cantor's paper Ein Beitrag zur Mannigfaltigkeitslehre, published in 1878 in volume 84 of Journal für die reine und angewandte Mathematik, pages 242–258.^[ref]

 Read the full article

Past exam of the computer science course of the University of Oxford / 2015 / Computer graphics / 1 / a / Solution 2026-10-31

 View more

Three advantages of integer arithmetic for raster graphics are:

Cheap, fast operations. Historically, integer addition, subtraction and shifts needed much less hardware and execution time than floating-point arithmetic, especially on machines without a floating-point unit. An incremental Bresenham's line algorithm needs only a few additions and a decision per pixel.
Exact decisions. Provided integer overflow is excluded, the scaled midpoint line decision is exact. There is no accumulated rounding error that can move a decision across a half-pixel boundary; equality has a reproducible tie rule.
Simple storage and hardware. Integer coordinates already identify pixels. An implementation can keep coordinates and decision variables in ordinary registers and use a small add/compare datapath, avoiding conversion between floating-point numbers and integer screen addresses.

These are historical and structural advantages, not a claim that every modern processor executes every integer operation faster. An integer implementation still needs sufficient width for its intermediate results. The compact answer is

fast elementary operations; exact decisions; simple hardware and storage .

(1)

 Read the full article

Past exam of the computer science course of the University of Oxford / 2015 / Computer graphics / 3 / c / Solution 2026-10-31

 View more

Let the oriented vertices be

a, b, c

, in an order such that

n = (b - a) \times (c - a)

points out of the front face. The given nondegeneracy guarantees

n \neq = 0

. A ray with nonzero direction

r

approaches the front face only if

n \cdot r < 0.

(1)

Reject

n \cdot r \geq 0

. This combines back-face culling with rejection of parallel rays, including the coplanar case for this ordinary surface-hit test. The ray-plane intersection parameter is

λ = \frac{n \cdot ( a - o )}{n \cdot r} .

(2)

Reject

λ < 0

; for rays launched from existing surfaces use a positive lower cutoff to avoid counting the launch surface. Form

p = o + λ r

Finally determine whether

p

lies inside the triangle. For

u = b - a

v = c - a

w = p - a

, compute

D = (u \cdot u) (v \cdot v) - (u \cdot v)^{2} > 0,

(3)

β = \frac{( v \cdot v ) ( w \cdot u ) - ( u \cdot v ) ( w \cdot v )}{D}, γ = \frac{( u \cdot u ) ( w \cdot v ) - ( u \cdot v ) ( w \cdot u )}{D} .

(4)

These solve

p = a + β u + γ v

. The barycentric coordinates are

(1 - β - γ, β, γ)

, so the complete acceptance criterion is

n \cdot r < 0, λ \geq 0, β \geq 0, γ \geq 0, β + γ \leq 1; p = o + λ r .

(5)

The inequalities include the edges and vertices. A renderer sharing edges between triangles should use a consistent numerical edge convention; small floating-point arithmetic errors near a shared boundary should not create cracks. Normalizing

n

is unnecessary for this intersection calculation because its scale cancels.

 Read the full article

Past exam of the computer science course of the University of Oxford / 2015 / Computer graphics / 3 / e / Solution 2026-10-31

 View more

At each primary hit

p

, send a shadow ray toward each point light. An especially clear parameterization is

s_{j} (t) = p + t (ℓ_{j} - p), 0 < t < 1.

(1)

If an opaque scene surface intersects this open segment, the source is hidden. The finite interval is essential: an object beyond the light must not cast a shadow at

p

. Exclude the launch surface with a small geometrically appropriate origin offset or positive minimum parameter, so floating-point arithmetic does not create a false self-intersection.

Introduce a visibility factor

V_{j} (p)

equal to zero when blocked and one otherwise. The modified equation is

I (p) = k_{a} I_{a} + j \sum V_{j} (p) I_{j} (p) [k_{d} max (0, \hat{N} \cdot \hat{L}_{j}) + k_{s} χ_{j} max (0, \hat{R}_{j} \cdot \hat{V})^{n}] .

(2)

A blocked light loses both its diffuse and specular contributions. The simple ambient approximation remains. Opaque point lights produce hard shadow boundaries. Use a consistent sidedness policy for occluders: an opaque closed solid blocks from either crossing side, so a two-sided blocker test is appropriate even when primary visibility uses back-face culling. An any-hit query can stop as soon as it finds one blocker before the light.

 Read the full article

Past exam of the computer science course of the University of Oxford / 2015 / Geometric modelling / 1 / e / Solution 2026-10-31

 View more

The tight axis-aligned bounding box is

Q = [0, 2] \times [0, 2]

. Bisect vertically at

x = 1

, giving a right child

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

. For a scalar implicit function representing the same constructive solid geometry solid, use

Φ = max {f_{5}, f_{6}, f_{4}, f_{1}, min (f_{2}, f_{3})}, P = {Φ \leq 0} .

(1)

With the convention

f_{i} \leq 0

inside, set intersection corresponds to the maximum, and set union to the minimum. These operations encode membership;

Φ

is not a signed distance function.

Since every

f_{i}

is an affine function, its extreme values on

R

occur at corners. Exact real-arithmetic interval arithmetic gives:

Function	Range on $R$	$f_{1}$	$[- 1 - 3, 3]$	$f_{2}$	$[- 1, 1]$	$f_{3}$	$[0, 1]$	$f_{4}$	$[- 2, 0]$	$f_{5}$	$[- 2, - 1]$	$f_{6}$	$[- 2, 0]$	$min (f_{2}, f_{3})$ , propagated	$[- 1, 1]$	$Φ$ , propagated	$[- 1, 3]$

Each leaf range is exact, but combining independent intervals forgets that the expressions share

x, y

. The final enclosure need not be the exact range of

Φ

. For example, equality

Φ = - 1

would require both

f_{6} \leq - 1

, hence

y \geq 1

, and

min (f_{2}, f_{3}) \leq - 1

, hence

y \leq 0

R

, which is impossible. The enclosure is still sufficient for a reliable three-way decision:

sup Φ (R) < 0 : inside; in f Φ (R) > 0 : outside; 0 \in Φ (R) : unresolved, so refine .

(2)

Here the result is unresolved. This is a real mixed box, rather than only a pessimistic test:

(5/4, 3/4)

is strictly inside, whereas

(3/2, 3/2)

is outside. A zero endpoint alone can also mean contact with the boundary, so equality must not be treated as proof of strict interior or strict exterior. In floating-point arithmetic, certified bounds require outward rounding.

Continue splitting unresolved boxes, normally alternating axes or splitting the longest side. Propagate the interval arithmetic results through the constructive solid geometry tree, discard exterior boxes and stop refining strict interior boxes. In sufficiently small boundary boxes, locate crossings of the active implicit equations with box edges and connect them into boundary line segments. Retain corners and resolve ambiguous multiple crossings by further refinement; here every genuine boundary piece is straight, so its supporting implicit equation allows exact intersection rather than a curved approximation. A maximum boundary-box diagonal

ε

gives spatial localization at that scale, provided no boundary components are discarded; corner-sign sampling alone provides no such guarantee.

The tight box saves empty-space work, but it is awkward for extracting facets. Entire boundary line segments, notably

EO

O A

and

D E

, lie on the initial box boundary, and

C D

lies on the first splitting straight line. A sign-change-only method may have no exterior sample on the other side and lose such pieces; exact zeros also persist through refinement. A slightly padded initial axis-aligned bounding box, together with an explicit rule for boundary-coincident pieces, makes the faceting more robust. The tight box is therefore economical for enclosure but not the best unqualified choice for a simple faceting implementation.

 Read the full article

Past exam of the computer science course of the University of Oxford / 2015 / Object oriented programming / 2 / f / Solution 2026-10-31

 View more

For each sample, include both its marker and its bar. A bar's lower vertical limit is

min (0, y_{i})

, including negative bars below the horizontal axis. Its horizontal half-width is

w_{i} /2

, rather than

w_{i}

. Its marker extends by

C

in both directions. Therefore that sample's axis-aligned bounding box is

[x_{i} - max (C, w_{i} /2), x_{i} + max (C, w_{i} /2)] \times [min (0, y_{i} - C), max (0, y_{i} + C)] .

(1)

Taking successive bounding unions includes every sample. Finally unite with

[- M, M] \times [- M, M]

to include the stipulated axes. The printed upper-corner formula has one extra closing parenthesis; both the lower-corner formula and the picture confirm the coordinatewise interpretation used here.^[ref]

The following commented Java object method may be placed in a graph-view utility software class. It uses Rect from part (a) and one iterator snapshot from part (e).

static Rect barBounds(DataList data, double C, double M) {
    if (!Double.isFinite(C) || !Double.isFinite(M) || C <= 0 || M <= 0)
        throw new IllegalArgumentException("positive finite margins required");
    List<Point> p = new ArrayList<Point>();
    for (Point point : data) p.add(point);
    int n = p.size();
    if (n < 2) throw new IllegalArgumentException("at least two points required");

    // This central square supplies the required minimum axis extents.
    Rect bounds = new Rect(new Point(-M, -M), new Point(M, M));
    for (int i = 0; i < n; ++i) {
        double x = p.get(i).getX(), y = p.get(i).getY();
        double width;
        if (i == 0) {
            width = p.get(1).getX() - x;
        } else if (i == n - 1) {
            width = x - p.get(i - 1).getX();
        } else {
            width = Math.min(x - p.get(i - 1).getX(),
                             p.get(i + 1).getX() - x);
        }
        // Include whichever is wider: marker or half of the bar.
        double half = Math.max(C, width / 2.0);
        // The bar reaches the baseline; the marker reaches y +/- C.
        Rect item = new Rect(new Point(x - half, Math.min(0.0, y - C)),
                             new Point(x + half, Math.max(0.0, y + C)));
        bounds = bounds.Union(item);
    }
    return bounds;
}

This uses the usual assumption that coordinate differences and expanded bounds are representable as finite double values. Rect rejects overflowed bounds; numerical range extension is not part of the requested geometry. As with all floating-point arithmetic, the calculation approximates the corresponding exact real number formulas.

Every side produced by the coordinatewise extrema is required by some bar, marker enclosure or axis extent. The result is therefore the smallest guaranteed rectangle for the specified marker enclosures and axes. If

X_{l}, Y_{l}, X_{h}, Y_{h}

denote the graph-only bounds, the final answer is

[min (- M, X_{l}), max (M, X_{h})] \times [min (- M, Y_{l}), max (M, Y_{h})] .

(2)

 Read the full article

Past exam of the mathematical and Computational Biology course of the University of Cambridge / 2024 / Paper 2 / A1 / a / Solution 2026-10-20

 View more

The expression is the Binet formula for the Fibonacci numbers. Compute the two roots and round the result to the nearest integer:

from math import sqrt


def binet(n):
    phi = (1 + sqrt(5)) / 2
    psi = (1 - sqrt(5)) / 2
    return round((phi**n - psi**n) / sqrt(5))

The returned object is an integer. This floating-point arithmetic implementation correctly rounds the requested range through

n = 50

; it is not an arbitrary-precision method for indefinitely large

n

. An integer recurrence is a useful independent numerical check.

 Read the full article

Past exam of the year 1 of the computer science course of the University of Oxford / 2015 / Imperative programming / 4 / a / Solution 2026-10-31

 View more

Represent the zero polynomial by degree = -1, with an allocated one-element array whose element is zero. Its coefficient prefix is empty; the spare cell is merely convenient storage. For every nonzero polynomial, require degree >= 0, sufficient allocated capacity, and a nonzero leading polynomial coefficient.

PROCEDURE Zero(VAR p: Poly);
BEGIN
  p.degree := -1; NEW(p.coeff, 1); p.coeff[0] := 0.0
END Zero;

PROCEDURE Normalize(VAR p: Poly);
BEGIN
  WHILE p.degree >= 0 DO
    IF p.coeff[p.degree] # 0.0 THEN RETURN END;
    DEC(p.degree)
  END
END Normalize;

Thus

de g_{stored} 0 = - 1; p \neq = 0 ⟹ p_{p . degree} \neq = 0.

(1)

The sentinel

- 1

is a representation convention, not an assertion that the mathematical degree of the zero polynomial is an ordinary integer. In parts (b)–(e), algebraic correctness uses exact arithmetic on real numbers; ordinary floating-point arithmetic introduces rounding and needs an additional numerical policy if exact zero tests are to be replaced.

 Read the full article

Past exam of the year 1 of the computer science course of the University of Oxford / 2015 / Imperative programming / 8 / a / Solution 2026-10-31

 View more

For Uniform, require a nonempty integer interval. Compute its size in Long before converting to Double, avoiding integer overflow in u-l+1. For Flip, require

0 \leq p \leq 1

, and omit zero-probability outcomes from the reported support.

class Uniform(l: Int, u: Int) extends FiniteDist[Int] {
  require(l <= u)
  private val size = u.toLong - l.toLong + 1L
  def getDom: List[Int] = {
    val result = List.newBuilder[Int]
    var x = l.toLong
    while (x <= u.toLong) { result += x.toInt; x += 1L }
    result.result()
  }
  def getProb(x: Int): Double =
    if (l <= x && x <= u) 1.0 / size.toDouble else 0.0
}

class Flip(p: Double) extends FiniteDist[Boolean] {
  require(0.0 <= p && p <= 1.0)
  def getDom: List[Boolean] =
    if (p == 0.0) List(false)
    else if (p == 1.0) List(true)
    else List(true, false)
  def getProb(x: Boolean): Double = if (x) p else 1.0-p
}

The first is a discrete uniform distribution; the second is a Bernoulli distribution on Boolean values. Their exact mathematical probability mass functions are

p_{U} (x) = ⎩ ⎨ ⎧ \frac{1}{u - l + 1} 0 l \leq x \leq u, otherwise, p_{F} (true) = p, p_{F} (false) = 1 - p .

(1)

Each support contains exactly the positive-probability outcomes once, and the probabilities sum to one. Double approximates these exact values by floating-point arithmetic. Materializing a huge support is limited by available memory; the Long calculation keeps its cardinality and endpoints arithmetically correct.

 Read the full article