Exact arithmetic, and what it costs instead

A fraction recovered from one remainder

A solution over the rationals can be computed modulo a prime power and then recovered — the residue determines the fraction uniquely, but only once the modulus is twice the square of the fraction's longer part. Below that there is no partial credit: the algorithm returns a different fraction with the same residue, and it is a perfectly good one.

Worth reading first: An answer with no error in it · Elimination is a sequence of choices · The condition number is an amplifier.

The Chinese remainder theorem reconstructs an integer from its residues. A determinant is an integer, so the modular route finishes there. A solution vector is not: every entry is a ratio of two determinants, and the reconstruction has one more step in it.

That step is remarkable and it is not obvious it should be possible at all. Given a single residue u = n/d modulo m — one integer — recover n and d. There are m fractions with that residue and the claim is that exactly one of them is small enough to be the answer.

Recovering 355/113 from its residue: nothing below 2^18, everything aboveThe fraction 355/113 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 18 the reconstruction returns a different fraction or refuses; at and above it, it returns 355/113 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^17.94, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^18. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^18modulus, as a power of twoa lattice with one short vectornumerator355denominator1132·max(n, d)², bits18first recovered at18moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 1 Recovering 355/113 from its residue modulo a prime just above 2ᵉ, for e from 3 to 44. Below 2¹⁸ the reconstruction returns something else; at and above it, 355/113 exactly, every time.

There is a second surprise in it, which is that the same trick has been in use for two thousand years under a different description. A continued-fraction expansion of a real number produces the best rational approximations to it with bounded denominators, and it is computed by the Euclidean algorithm; recovering n/d from its residue is that expansion run on a pair of integers instead of on a real, stopped by a size condition instead of by an accuracy one. The convergents of π are the canonical example and 355/113 is the third of them, which is why it is the fraction the figure recovers.

The algorithm is one that is already there

Run the extended Euclidean algorithm on (m, u) and stop it early — the first time the remainder falls below √(m/2). At that point the remainder is the numerator and the accumulated coefficient is the denominator, up to sign.

That is the whole method. It is the same sequence of divisions that computes a greatest common divisor, run on a pair it will not be allowed to finish, and its intermediate state is a sequence of pairs (rᵢ, tᵢ) with rᵢ ≡ tᵢ·u (mod m) at every step. The stopping rule chooses which of those pairs to return.

Two things about that are worth noticing.

The cost is a gcd, which the field has already priced. The first essay in this field measured the price of reducing fractions during an elimination and found it to be a gcd per operation. Here the same routine appears once per entry of the answer rather than once per arithmetic operation, which is n² times less often, and it is the reason the modular route reaches an exact solution without ever forming a fraction in the middle.

And it needs no search. There are m candidate fractions and the algorithm looks at about log m of them, because the Euclidean sequence’s remainders fall geometrically. That is not a heuristic that usually works: the pairs it skips are provably too large, so the one it stops at is the only one it could return.

The implementation has one place where it can quietly do the wrong thing, and it is the return value on failure. When no pair meets the size condition — because the modulus was too small for this fraction — there is still a pair to return, and returning it is the natural thing for the code to do. The routine here returns nothing instead, and says so, because a reconstruction that cannot say the modulus was too small converts a recoverable failure into a silent one. That is the same discipline the interval-arithmetic essays apply to a verdict of nothing can be said: it is an outcome, it is useful, and a routine that has no way to express it will express something else.

Why exactly one fraction is small enough

The uniqueness is a statement about a lattice, and it is short.

The set of integer pairs (a, b) with a ≡ b·u (mod m) is a lattice in the plane — closed under addition and integer multiples — and its determinant is m. The answer (n, d) is one of its points. Any other point (n′, d′) with the same residue is also in it, and if both are shorter than √(m/2) in each coordinate then the parallelogram they span has area at most 2·(m/2) = m, with equality only in degenerate cases; a lattice of determinant m has no two independent points spanning less than that. So the short point is unique, and finding it is finding the short vector in a two-dimensional lattice — which is what the Euclidean algorithm does.

That is the same object the lattice essay in this field is about, in the smallest dimension where it has anything to say. Reconstruction is lattice reduction in two dimensions, and in two dimensions the reduction is exact, terminating, and older than the problem it is being used on.

The bound is about the larger part, and that is the mistake to avoid

The natural way to write the bound down is wrong, and the measurement is what said so.

Since there are about m fractions and the answer has parts n and d, the tempting statement is that the modulus must exceed 2·|n|·d — the product. It reads well and it is false. The condition the algorithm actually needs is that both parts fall below √(m/2) separately, which requires

m  >  2·max(|n|, d)²

and those two are different whenever the parts are of different sizes. For 355/113 the product bound predicts a crossover at 2¹², and the measured first success is 2¹⁸ — six bits out, in the direction that would have made a routine budget too little modulus and return the wrong fraction with no warning.

Recovering 355/7 from its residue: nothing below 2^18, everything aboveThe fraction 355/7 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 18 the reconstruction returns a different fraction or refuses; at and above it, it returns 355/7 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^17.94, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^18. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^18modulus, as a power of twoa lattice with one short vectornumerator355denominator72·max(n, d)², bits18first recovered at18moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 2 355/7, whose product is fifty times smaller and whose crossover is in exactly the same place. The bound does not know about the denominator until the denominator is the larger part.

The drag range makes the correction visible in one motion: the cliff does not move at all while the denominator is below the numerator, and then moves two bits for every bit the denominator gains. A bound proportional to the product would have moved throughout.

There is no partial credit

The failure below the bound is the thing this essay is really about, and it is unlike every other failure on this site.

Below the crossover the algorithm does not return an approximation of 355/113. It returns a different fraction — a genuinely correct one, in the sense that its residue modulo m is the same u, and one with smaller parts. There is no sense in which it is nearly right. 355/113 and the fraction returned at 2¹² are not close, they are not related, and no measurement of the returned value would suggest anything was wrong.

Compare that with every truncation elsewhere on this site. A truncated model is wrong by an amount the discarded singular values bound. A restarted Krylov method loses accuracy in a way its own residual estimate tracks. Even a rank chosen at the wrong threshold is wrong by one, and the singular value it crossed says by how much. All of those degrade.

This does not degrade, it substitutes — which puts it in the small class of failures the graph field opened with, where a rounding does not make an answer slightly wrong but makes it a different answer. The only defence is arithmetic: know how large the parts can be, and supply a modulus twice the square of that.

The measurement, and what it took to state it

The figure sweeps a fraction against forty-two moduli and records, for each, whether the reconstruction returned that fraction. It is a binary outcome per modulus, which makes the picture a pair of blocks rather than a curve — and the sharpness of the boundary between them is the finding.

Between 2³ and 2¹⁷ inclusive, every modulus fails. At 2¹⁸ and every modulus above it, up to 2⁴⁴, the reconstruction is exact. There is no modulus at which it sometimes works, none at which it returns a fraction close to the answer, and none at which the failure is marginal. The assertion in the figure checks that the first success is where the lattice bound puts it, within two bits, and it is: 2¹⁸ against a bound of 2¹⁷·⁹.

That is a stronger claim than the algorithm works when the modulus is large enough and it is worth the extra apparatus, because the weaker claim is compatible with a routine that is unreliable near the boundary — which is what a reader would reasonably assume, since almost everything else here is.

Knowing how large the parts can be

Which is answerable, and by the same theorem that budgeted the primes.

By Cramer’s rule every entry of the solution is det(Aᵢ)/det(A), so the denominator is bounded by Hadamard’s bound on A and the numerator by Hadamard’s bound on the column-replaced matrix. Both are computable in n² operations before anything is solved. Take the larger, square it, double it, and that is the modulus the lifting has to reach.

So the whole modular solve is budgeted in advance from end to end: how many p-adic digits to lift, because the modulus must reach 2·max(bound)²; and then one Euclidean run per entry. The bounds are generous — Hadamard’s is slack by whatever cancellation the determinant had in it, and the answer’s parts are usually far below their ceilings — so a routine sized this way does noticeably more lifting than it needed. That is the same trade the prime budget made and for the same reason: the only alternative is to watch the reconstruction and stop when it settles, and a reconstruction that settles has coincided rather than converged. Nothing in it adapts, nothing watches itself, and as the previous essay noted that is the opposite of how every iterative method here terminates.

The exact solution of a random system is longer than the systemAx = b with A a random integer matrix and b chosen so the answer is a genuine vector of fractions. At n = 11 the question occupies 484 bits — 11² entries of 4 bits each — and the answer occupies 729, a ratio of 1.51. Every entry is a ratio of two determinants by Cramer's rule, checked here by multiplying each one back by det A and comparing it with the determinant of the column-replaced matrix as an integer — so every denominator divides det A and none of them is longer than its 33 bits, while a numerator can be, at 37 here. The answer IS that long, whatever route produced it. That is a floor no exact solver can get under, and it is why the modular route reconstructs the answer rather than eliminating its way to it.3456789101110²10³nbitsthe answerthe question · det Abits in, bits outn11question484answer729det A33widest part37the answer is the floorand no route writes it down more cheaply
Fig. 3 The parts being bounded: the answer’s own length against the question’s, where every denominator divides det A and the numerators can exceed it.

The consequence for cost is a factor of two that is easy to miss. The modulus must be twice the square of the longer part, so it has about twice as many bits as the answer does. A p-adic solve that produces k correct digits recovers fractions with about k/2 digits in each part, and a routine sized against the answer’s own length rather than against twice it will lift exactly half as far as it needs to and reconstruct confidently wrong fractions.

Where the ratio of the parts matters and where it does not

The bound being about the larger part has a consequence for how a solve should be posed, and it is the kind of thing that only shows up once the bound is written correctly.

An entry of the solution whose parts are 10²⁰ and 3 costs the same modulus as one whose parts are 10²⁰ and 10²⁰. So a system whose answer happens to have small denominators is not cheaper to reconstruct, and the common intuition — nice answers are cheap — is exactly half true: a nice answer is cheap only if both parts are small.

Where it bites is a system with a mixed answer, some entries large and some tiny. The modulus is sized against the worst entry and every entry is reconstructed at that modulus, so the cost is set by one component. That is the same shape as a stopping test governed by whichever quantity is slowest, except that here there is nothing to race: the requirement is a maximum over the components and it is computed rather than observed.

The repair, where one is available, is to reconstruct entry by entry and stop lifting for the entries that have already been recovered and verified. It works, and it changes the routine’s character from a fixed budget into an adaptive one — which reintroduces exactly the thing the budget existed to avoid, since an entry that reconstructs to the same fraction at two successive moduli has not been verified, it has coincided.

What the failure looks like from outside

It is worth being concrete about what a caller sees, because the answer is: nothing.

The reconstruction returns a fraction. That fraction, substituted into the original system, does not satisfy it — so the failure is detectable, by exactly one test: form A·y − d·b in integers and check that it is the zero vector. That is cheap, it is exact, and it is the check the essay on the answer’s length argues should be the return type’s natural form anyway.

Without it there is no signal. The fraction has plausible-sized parts, its residue is the right residue, and there is no residual to be small or large because nothing was rounded. A routine that lifts to a modulus it guessed at and does not verify has a silent failure mode with no observable symptom, and the observable symptom exists and costs one matrix–vector product.

The whole route, end to end

It is worth assembling the three modular essays into the single algorithm they describe, because each of them is a component and none of them is the method.

Given an integer matrix A, an integer right-hand side b, and a request for the exact rational solution:

  1. Bound the answer. Hadamard’s bound on A gives a ceiling on every denominator; Hadamard’s bound on each column-replaced matrix gives one on every numerator. Take the larger, call it B.
  2. Choose a prime and check it. A prime p is usable when A is invertible modulo p, which one elimination decides. An unlucky prime is rejected here rather than discovered later.
  3. Lift. Solve modulo p, correct, solve again — doubling the number of correct p-adic digits each round — until the modulus passes 2B². Every operation in this stage is on machine words.
  4. Reconstruct. One truncated Euclidean run per entry recovers the fraction from its residue.
  5. Verify. Form A·y − d·b in integers and check it is the zero vector.

Step 1 is the budget, step 2 is the unlucky prime, steps 3 and 4 are this essay, and step 5 is the only part that is not an algorithm at all — it is the check that turns the answer into a certificate, and it costs one matrix–vector product in exact integers.

What is notable about the list is that nothing in it adapts. Every quantity that decides how much work happens is computed before any work happens, from bounds rather than from observations. A routine built this way either returns the exact answer or fails a check, and it never returns something approximately right — which is the whole difference between this field and the rest of the site, stated as a control flow rather than as a philosophy.

At other settings

Recovering 355/9973 from its residue: nothing below 2^28, everything aboveThe fraction 355/9973 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 28 the reconstruction returns a different fraction or refuses; at and above it, it returns 355/9973 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^27.57, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^28. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^28modulus, as a power of twoa lattice with one short vectornumerator355denominator99732·max(n, d)², bits28first recovered at28moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 4 A denominator larger than the numerator, where the cliff has moved right by the difference.
Recovering 1/99991 from its residue: nothing below 2^35, everything aboveThe fraction 1/99991 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 35 the reconstruction returns a different fraction or refuses; at and above it, it returns 1/99991 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^34.22, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^35. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^35modulus, as a power of twoa lattice with one short vectornumerator1denominator10·10⁴2·max(n, d)², bits34first recovered at35moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 5 And the extreme case, a unit numerator, where the bound is entirely about the denominator.
Recovering 65537/29 from its residue: nothing below 2^34, everything aboveThe fraction 65537/29 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 34 the reconstruction returns a different fraction or refuses; at and above it, it returns 65537/29 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^33.00, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^34. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^34modulus, as a power of twoa lattice with one short vectornumerator6.6·10⁴denominator292·max(n, d)², bits33first recovered at34moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 6 A large numerator and a small denominator, which the product bound would have priced as the same problem as the last one and is not.
Recovering 355/113 from its residue: nothing below 2^18, everything aboveThe fraction 355/113 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 18 the reconstruction returns a different fraction or refuses; at and above it, it returns 355/113 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^17.94, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^18. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^18modulus, as a power of twoa lattice with one short vectornumerator355denominator1132·max(n, d)², bits18first recovered at18moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 7 The classical approximation to π, at every modulus from 2³ to 2⁴⁴.
The exact solution of a wide system is longer than the systemAx = b with A a wide integer matrix and b chosen so the answer is a genuine vector of fractions. At n = 10 the question occupies 1000 bits — 10² entries of 10 bits each — and the answer occupies 2058, a ratio of 2.06. Every entry is a ratio of two determinants by Cramer's rule, checked here by multiplying each one back by det A and comparing it with the determinant of the column-replaced matrix as an integer — so every denominator divides det A and none of them is longer than its 103 bits, while a numerator can be, at 105 here. The answer IS that long, whatever route produced it. That is a floor no exact solver can get under, and it is why the modular route reconstructs the answer rather than eliminating its way to it.34567891010²10³nbitsthe answerthe question · det Abits in, bits outn10question1000answer2058det A103widest part105the answer is the floorand no route writes it down more cheaply
Fig. 8 Three-digit entries, where both parts of every entry are longer and the modulus needed doubles with them.
One lattice, two bases: an orthogonality defect of 12.04 reduced to 1The dots are the lattice generated by (1, 0) and (12, 1) — every one of them, out to ±4 in each coefficient. The long pair of arrows is that basis; the short pair is what LLL returns after 1 steps, in exact rational arithmetic with δ = 3/4. The two bases generate the same set of points, which is checked rather than drawn: the lattice determinant is 1 before and 1 after, and every step of the reduction is a unimodular operation, which is what makes that so. What the reduction changes is the orthogonality defect — the product of the basis vectors' lengths over the determinant — from 12.042 to 1. One is the floor, and it is attained only by an orthogonal basis.the same lattice, twicewhat the reduction may not changedet, before1det, after1defect, before12defect, after1LLL steps1the determinant is the invariantand the defect is what is being reduced
Fig. 9 The two-dimensional lattice reduction this reconstruction is, in the field’s own picture of it.
One lattice, two bases: an orthogonality defect of 40.01 reduced to 1The dots are the lattice generated by (1, 0) and (40, 1) — every one of them, out to ±4 in each coefficient. The long pair of arrows is that basis; the short pair is what LLL returns after 1 steps, in exact rational arithmetic with δ = 3/4. The two bases generate the same set of points, which is checked rather than drawn: the lattice determinant is 1 before and 1 after, and every step of the reduction is a unimodular operation, which is what makes that so. What the reduction changes is the orthogonality defect — the product of the basis vectors' lengths over the determinant — from 40.012 to 1. One is the floor, and it is attained only by an orthogonal basis.the same lattice, twicewhat the reduction may not changedet, before1det, after1defect, before40defect, after1LLL steps1the determinant is the invariantand the defect is what is being reduced
Fig. 10 And at a skew of forty, where the short vector is very much shorter than the basis given.
The determinant's length, and the bound that decides how many primes to use before any of them is chosenSolid: Hadamard's bound on |det A|, the product of the row 2-norms, in bits. Dotted: the determinant's actual length. The bound is what a modular determinant budgets against — at n = 11 the random family needs 47 bits, which is 2 word-sized primes at 26 bits each, and that count is fixed before the first residue is computed. On a Sylvester Hadamard matrix of order 8 the bound is 13 bits and the determinant has 13, so the bound is attained rather than merely satisfied. On the unit triangular family the determinant is 1 at every size while the bound grows, which is the gap the budget pays for: a matrix can have a determinant of one and rows that are long.34567891011110¹10²nbitsthe budget and what it buysrandom, bound47random, actual33primes needed2Hadamard n = 8, bound13Hadamard n = 8, actual13the count is decided by a theorembefore any arithmetic happens
Fig. 11 Hadamard’s bound, which supplies the numbers this reconstruction has to be sized against.
det A = 30030, and its residue at 25 primes — 5 of which divide itEach dot is det A mod p, drawn as a fraction of p so the primes are comparable. A dot on the floor is a prime at which the elimination reports a singular matrix — and the matrix is not singular: its determinant is 30030. 5 of 25 primes here are unlucky, a rate of 20 per cent, which is what it looks like when the primes are small. At the 26-bit primes a modular determinant actually uses, a prime divides a 15-bit determinant with probability of order 2.24·10⁻⁷ — but the consequence is unchanged, and there is no way to see it from one residue. What sees it is the budget: computing more residues than Hadamard's bound requires, and refusing to reconstruct when the extra ones disagree.det mod p, as a fraction of p3713a prime that divides the answerdet A3·10⁴primes swept25unlucky5rate0.2det, in bits15singular mod p is not singularand one residue cannot tell them apart
Fig. 12 The other way a modular computation returns a correct answer to the wrong question.
10 × 10: three exact determinants agree to the digit, and the float one does notdet A = -77985307738 on a random matrix of order 10, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 2 word-sized primes put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns -7.79853077·10¹⁰, a relative error of 3.913·10⁻¹⁶, which is what a determinant of 37 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 43 bits are needed and 2 primes supply them.Bareiss37 bitsrationals36 bitsmodular26 bitsthe answer37 bitsHadamard43 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 2 primes0floating-point LU3.9·10⁻¹⁶det A = -77985307738three routes, one integerand no tolerance anywhere in the comparison
Fig. 13 The determinant every denominator divides, computed three ways.
An exact solve of a measured 8 × 8 system, against the error the measurement put thereThe rising curve is the distance from the exact solution of a perturbed system to the exact solution of the unperturbed one, as the perturbation grows from a relative 10⁻¹⁶ to 1. Every point on it was computed with no rounding anywhere: the elimination is over the rationals, the residual of each solve is the zero vector, and the error is entirely the data's. At a measurement accurate to 10^-12 the answer is wrong by 9.988·10⁻⁴, which is 10·10⁸ times the data's own error — the matrix's condition number, doing the only job left for it to do. The straight line of slope one is the whole finding: removing the arithmetic error changes nothing about this quantity, because this quantity was never about the arithmetic.10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁷10⁻⁵10⁻³10⁻¹10¹relative error in the datarelative error in the answerthe answer's errorthe data's errornothing here is the arithmetic'sdata error10⁻¹²answer error10·10⁻⁴amplification10·10⁸rounding committed0residual0the residual is exactly zeroand the answer is wrong anyway
Fig. 14 And the question of whether the exact answer was worth reconstructing at all.
Iterative refinement of the inverse-and-multiply solve at κ = 10^14The backward error starts at 4.5·10⁻⁶ and falls by about κu a step — 3.2·10⁻¹⁰, 1.4·10⁻¹³, 10⁻¹⁷, 2·10⁻¹⁷, 1.1·10⁻¹⁷ — reaching the LU route's 4.2·10⁻¹⁷ after 5 corrections, each costing 2n² flops against the 2n³ the inversion cost. The forward error, drawn above, does not improve: 0.0043 to 2.3·10⁻⁴, against the LU route's 2.8·10⁻⁴. What refinement at working precision buys is stability, and the accuracy floor belongs to the problem.01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 4.2·10⁻¹⁷LU route: forward 2.8·10⁻⁴forward errorbackward errorwhat a correction buysη before refinement4.5·10⁻⁶η after four steps1.1·10⁻¹⁷forward, unchanged2.3·10⁻⁴cost of a step, flops128the residual is repairableand the accuracy floor is the problem's
Fig. 15 Iterative refinement, whose p-adic lifting is the analogue and whose failure mode is not.
Iterative refinement from a 24-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 16 Buying accuracy back in floating point, where the recovery is gradual rather than a cliff.
Recovering 22/7 from its residue: nothing below 2^10, everything aboveThe fraction 22/7 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 10 the reconstruction returns a different fraction or refuses; at and above it, it returns 22/7 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^9.92, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^10. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^10modulus, as a power of twoa lattice with one short vectornumerator22denominator72·max(n, d)², bits9.9first recovered at10moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 17 A small fraction, recovered at 2^9 rather than 2^18.
Recovering 9973/113 from its residue: nothing below 2^28, everything aboveThe fraction 9973/113 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 28 the reconstruction returns a different fraction or refuses; at and above it, it returns 9973/113 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^27.57, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^28. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^28modulus, as a power of twoa lattice with one short vectornumerator9973denominator1132·max(n, d)², bits28first recovered at28moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 18 A four-digit numerator, where the cliff has moved by eight bits.
Recovering 355/1013 from its residue: nothing below 2^21, everything aboveThe fraction 355/1013 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 21 the reconstruction returns a different fraction or refuses; at and above it, it returns 355/1013 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^20.97, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^21. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^21modulus, as a power of twoa lattice with one short vectornumerator355denominator10132·max(n, d)², bits21first recovered at21moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 19 And a denominator larger than the numerator, where it moves again.
The exact solution of a random system is longer than the systemAx = b with A a random integer matrix and b chosen so the answer is a genuine vector of fractions. At n = 9 the question occupies 324 bits — 9² entries of 4 bits each — and the answer occupies 454, a ratio of 1.4. Every entry is a ratio of two determinants by Cramer's rule, checked here by multiplying each one back by det A and comparing it with the determinant of the column-replaced matrix as an integer — so every denominator divides det A and none of them is longer than its 25 bits, while a numerator can be, at 28 here. The answer IS that long, whatever route produced it. That is a floor no exact solver can get under, and it is why the modular route reconstructs the answer rather than eliminating its way to it.345678910²nbitsthe answerthe question · det Abits in, bits outn9question324answer454det A25widest part28the answer is the floorand no route writes it down more cheaply
Fig. 20 The parts being recovered, and how long they are.
The determinant's length, and the bound that decides how many primes to use before any of them is chosenSolid: Hadamard's bound on |det A|, the product of the row 2-norms, in bits. Dotted: the determinant's actual length. The bound is what a modular determinant budgets against — at n = 10 the random family needs 42 bits, which is 2 word-sized primes at 26 bits each, and that count is fixed before the first residue is computed. On a Sylvester Hadamard matrix of order 8 the bound is 13 bits and the determinant has 13, so the bound is attained rather than merely satisfied. On the unit triangular family the determinant is 1 at every size while the bound grows, which is the gap the budget pays for: a matrix can have a determinant of one and rows that are long.345678910110¹10²nbitsthe budget and what it buysrandom, bound42random, actual37primes needed2Hadamard n = 8, bound13Hadamard n = 8, actual13the count is decided by a theorembefore any arithmetic happens
Fig. 21 The bound that supplies the modulus this reconstruction is sized against.
8 × 8: three exact determinants agree to the digit, and the float one does notdet A = -8650618 on a random matrix of order 8, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 2 word-sized primes put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns -8.650618·10⁶, a relative error of 2.153·10⁻¹⁶, which is what a determinant of 24 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 33 bits are needed and 2 primes supply them.Bareiss24 bitsrationals23 bitsmodular26 bitsthe answer24 bitsHadamard33 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 2 primes0floating-point LU2.2·10⁻¹⁶det A = -8650618three routes, one integerand no tolerance anywhere in the comparison
Fig. 22 The determinant every denominator divides.
det A = 30030, and its residue at 36 primes — 5 of which divide itEach dot is det A mod p, drawn as a fraction of p so the primes are comparable. A dot on the floor is a prime at which the elimination reports a singular matrix — and the matrix is not singular: its determinant is 30030. 5 of 36 primes here are unlucky, a rate of 13.9 per cent, which is what it looks like when the primes are small. At the 26-bit primes a modular determinant actually uses, a prime divides a 15-bit determinant with probability of order 2.24·10⁻⁷ — but the consequence is unchanged, and there is no way to see it from one residue. What sees it is the budget: computing more residues than Hadamard's bound requires, and refusing to reconstruct when the extra ones disagree.det mod p, as a fraction of p311a prime that divides the answerdet A3·10⁴primes swept36unlucky5rate0.14det, in bits15singular mod p is not singularand one residue cannot tell them apart
Fig. 23 The prime that has to be checked before any lifting starts.
One lattice, two bases: an orthogonality defect of 20.02 reduced to 1The dots are the lattice generated by (1, 0) and (20, 1) — every one of them, out to ±4 in each coefficient. The long pair of arrows is that basis; the short pair is what LLL returns after 1 steps, in exact rational arithmetic with δ = 3/4. The two bases generate the same set of points, which is checked rather than drawn: the lattice determinant is 1 before and 1 after, and every step of the reduction is a unimodular operation, which is what makes that so. What the reduction changes is the orthogonality defect — the product of the basis vectors' lengths over the determinant — from 20.025 to 1. One is the floor, and it is attained only by an orthogonal basis.the same lattice, twicewhat the reduction may not changedet, before1det, after1defect, before20defect, after1LLL steps1the determinant is the invariantand the defect is what is being reduced
Fig. 24 The two-dimensional lattice this reconstruction reduces, at a larger skew.
4 × 4: the invariant factors, computed twice by routes with no algorithm in commonThe Smith normal form of a random matrix of order 4 is the diagonal [1, 1, 1, 4386], reached by unimodular row and column operations, with each entry dividing the next — a property the elimination does not give and which has to be enforced and asserted. The same numbers computed from the definition — the gcd of all k × k minors, over the gcd of all (k−1) × (k−1) minors — are [1, 1, 1, 4386], and the two agree entry for entry. Their product is 4386, which is |det A| = 4386. The Hermite form of the same matrix reaches 13 bits against the Smith form's 13, and its transform has determinant -1, so the lattice the rows generate is unchanged.s10 bitss20 bitss30 bitss412 bitsinvariant factors, in bitstwo routes, and a conserved quantitydet A-4386Π invariants4386SNF widest13HNF widest13det of the transform-1an algorithm and a definitionagreeing as integers
Fig. 25 The invariant factors that decide which primes are unusable.
Three exact eliminations of the same wide matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 9 the answer itself is 91 bits and Hadamard's bound allows 98. Fraction-free elimination never forms a number wider than 91 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 89 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 182755 bits — 2010 times the width of the fraction-free route and 2010 times the width of the answer. All three are correct.3456789110¹10²10³10⁴10⁵nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn9answer91Hadamard bound98fraction-free91reduced rationals89unreduced1.8·10⁵the error is zero on every curvethe cost is the length of the numbers
Fig. 26 The growth the p-adic route exists to avoid.
Recovering 1013/29 from its residue: nothing below 2^21, everything aboveThe fraction 1013/29 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 21 the reconstruction returns a different fraction or refuses; at and above it, it returns 1013/29 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^20.97, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^21. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^21modulus, as a power of twoa lattice with one short vectornumerator1013denominator292·max(n, d)², bits21first recovered at21moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 27 A large numerator with a small denominator.
Recovering 65537/113 from its residue: nothing below 2^34, everything aboveThe fraction 65537/113 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 34 the reconstruction returns a different fraction or refuses; at and above it, it returns 65537/113 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^33.00, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^34. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^34modulus, as a power of twoa lattice with one short vectornumerator6.6·10⁴denominator1132·max(n, d)², bits33first recovered at34moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 28 A five-digit numerator, at 2^35.
The exact solution of a hilbert system is longer than the systemAx = b with A a hilbert integer matrix and b chosen so the answer is a genuine vector of fractions. At n = 7 the question occupies 931 bits — 7² entries of 19 bits each — and the answer occupies 135, a ratio of 0.145. Every entry is a ratio of two determinants by Cramer's rule, checked here by multiplying each one back by det A and comparing it with the determinant of the column-replaced matrix as an integer — so every denominator divides det A and none of them is longer than its 49 bits, while a numerator can be, at 16 here. The answer IS that long, whatever route produced it. That is a floor no exact solver can get under, and it is why the modular route reconstructs the answer rather than eliminating its way to it.3456710²nbitsthe answerthe question · det Abits in, bits outn7question931answer135det A49widest part16the answer is the floorand no route writes it down more cheaply
Fig. 29 The parts of an ill-conditioned system’s answer.
8 × 8: three exact determinants agree to the digit, and the float one does notdet A = 3458508270104059228258152 on a wide matrix of order 8, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 4 word-sized primes put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns 3.45850827·10²⁴, a relative error of 3.105·10⁻¹⁶, which is what a determinant of 82 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 87 bits are needed and 4 primes supply them.Bareiss82 bitsrationals77 bitsmodular26 bitsthe answer82 bitsHadamard87 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 4 primes0floating-point LU3.1·10⁻¹⁶det A = 3458508270104059228258152three routes, one integerand no tolerance anywhere in the comparison
Fig. 30 The determinant every denominator divides.
The determinant's length, and the bound that decides how many primes to use before any of them is chosenSolid: Hadamard's bound on |det A|, the product of the row 2-norms, in bits. Dotted: the determinant's actual length. The bound is what a modular determinant budgets against — at n = 12 the random family needs 52 bits, which is 3 word-sized primes at 26 bits each, and that count is fixed before the first residue is computed. On a Sylvester Hadamard matrix of order 8 the bound is 13 bits and the determinant has 13, so the bound is attained rather than merely satisfied. On the unit triangular family the determinant is 1 at every size while the bound grows, which is the gap the budget pays for: a matrix can have a determinant of one and rows that are long.3456789101112110¹10²nbitsthe budget and what it buysrandom, bound52random, actual40primes needed3Hadamard n = 8, bound13Hadamard n = 8, actual13the count is decided by a theorembefore any arithmetic happens
Fig. 31 Hadamard, bounding the parts this reconstruction is sized against.
One lattice, two bases: an orthogonality defect of 5.099 reduced to 1The dots are the lattice generated by (1, 0) and (5, 1) — every one of them, out to ±4 in each coefficient. The long pair of arrows is that basis; the short pair is what LLL returns after 1 steps, in exact rational arithmetic with δ = 3/4. The two bases generate the same set of points, which is checked rather than drawn: the lattice determinant is 1 before and 1 after, and every step of the reduction is a unimodular operation, which is what makes that so. What the reduction changes is the orthogonality defect — the product of the basis vectors' lengths over the determinant — from 5.099 to 1. One is the floor, and it is attained only by an orthogonal basis.the same lattice, twicewhat the reduction may not changedet, before1det, after1defect, before5.1defect, after1LLL steps1the determinant is the invariantand the defect is what is being reduced
Fig. 32 The lattice, at a small skew.
5 × 5: the invariant factors, computed twice by routes with no algorithm in commonThe Smith normal form of a wide matrix of order 5 is the diagonal [1, 1, 1, 1, 2580458337283640], reached by unimodular row and column operations, with each entry dividing the next — a property the elimination does not give and which has to be enforced and asserted. The same numbers computed from the definition — the gcd of all k × k minors, over the gcd of all (k−1) × (k−1) minors — are [1, 1, 1, 1, 2580458337283640], and the two agree entry for entry. Their product is 2580458337283640, which is |det A| = 2580458337283640. The Hermite form of the same matrix reaches 106 bits against the Smith form's 52, and its transform has determinant -1, so the lattice the rows generate is unchanged.s10 bitss20 bitss30 bitss40 bitss551 bitsinvariant factors, in bitstwo routes, and a conserved quantitydet A-2.6·10¹⁵Π invariants2.6·10¹⁵SNF widest52HNF widest106det of the transform-1an algorithm and a definitionagreeing as integers
Fig. 33 The invariant factors that decide which prime may be lifted from.
det A = 30030, and its residue at 44 primes — 5 of which divide itEach dot is det A mod p, drawn as a fraction of p so the primes are comparable. A dot on the floor is a prime at which the elimination reports a singular matrix — and the matrix is not singular: its determinant is 30030. 5 of 44 primes here are unlucky, a rate of 11.4 per cent, which is what it looks like when the primes are small. At the 26-bit primes a modular determinant actually uses, a prime divides a 15-bit determinant with probability of order 2.24·10⁻⁷ — but the consequence is unchanged, and there is no way to see it from one residue. What sees it is the budget: computing more residues than Hadamard's bound requires, and refusing to reconstruct when the extra ones disagree.det mod p, as a fraction of p311a prime that divides the answerdet A3·10⁴primes swept44unlucky5rate0.11det, in bits15singular mod p is not singularand one residue cannot tell them apart
Fig. 34 The prime that must be rejected before lifting starts.

What links here

Computed from the collection, not written here: the essays that point at this one.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Bit lengthChinese remainder theoremCramers ruleExact arithmeticLatticeModular arithmeticRational reconstruction