How many primes the answer needs

The first three essays in this field all had the same shape of problem: exact arithmetic on integers that grow, with the argument being about how much they grow and whose fault it is. There is a way to remove the growth entirely, and it is not a cleverer elimination.

Work in a ring where the numbers cannot grow. Choose a prime p below 2²⁶, reduce every entry of the matrix modulo p, and eliminate. Every quantity that appears is a residue in [0, p), a product of two of them fits exactly in a double, and there is no bit-length curve to draw because the bit length is 26 from the first step to the last, whatever the matrix was.

The answer is then wrong, in a specific and recoverable way: what has been computed is the determinant of A over the field with p elements, which is det A reduced modulo p. One residue is not the determinant. Several of them are.

That is a completely different bargain from the two the field has already made. Fraction-free elimination keeps the intermediates near the answer’s own length and pays for it with a long division per entry; reduced-rational elimination keeps them there and pays with a greatest common divisor. Both are still doing arithmetic on numbers whose length is decided by the matrix. This one fixes the length at twenty-six bits and moves the entire cost into how many times the elimination is run.

Putting it back together #

The Chinese remainder theorem says that residues modulo pairwise coprime moduli determine a unique value modulo their product. So computing det A mod p₁, det A mod p₂, … and combining them recovers det A modulo p₁p₂…, and once that product exceeds the range det A can lie in, the value is det A.

Two details in that sentence carry all the weight.

The range is signed. A determinant can be negative and the residues cannot, so the reconstruction has to be read in the symmetric range [−M/2, M/2) rather than in [0, M). Reading it in the wrong one reports a determinant of −6 as a twenty-digit integer, which is not a small error and is not detectable from inside: the arithmetic is right, the interpretation is not.

“Once the product exceeds the range” has to be decided in advance. There is no way to look at a partial reconstruction and tell whether it has settled. Adding a prime always changes the value — that is what a larger modulus does — and a value that stops changing has not converged, it has coincided. The number of primes is not something the computation can discover about itself.

That second point is the interesting one, and it is the opposite of how almost every other computation on this site terminates. An iterative solver stops when its residual is small enough; the essay on what a method reports is about the gap between that number and the real one, and the essay on stopping tests is about two quantities racing to decide it. Every one of those is a routine watching itself and deciding when to stop. A modular reconstruction cannot watch itself at all. It has to be told, before it starts, how much work it is going to do.

The bound, which is a theorem #

Hadamard’s inequality supplies the missing number: for any real matrix,

|det A|  ≤  Π ‖row i‖₂

the product of the row two-norms. That is computable from the entries in n² operations, before any elimination, and it says how many bits det A can occupy. Divide by 26 and round up and that is how many word-sized primes are needed.

Across sizes of 5, 7, 9, 10, 11, 12 and 13 the bound reads 18, 27, 37, 42, 47, 52 and 56 bits, the determinant actually needs 15, 24, 25, 37, 33, 40 and 44, and the number of primes budgeted is 1, 2, 2, 2, 2, 3 and 3.

The gap between the two columns is not a constant and cannot be, which is the practical content of the whole page. It is 3, 3, 12, 5, 14, 12 and 12 bits — the bound over-estimates by between 13% and 48% — because the bound is a function of the entries’ sizes alone and the answer’s length is a property of the particular matrix. A determinant that happens to come out small is a determinant the bound could not have known about.

And the primes are budgeted from the bound rather than from the answer, which is what makes the gap cost something rather than merely being untidy. At n = 9 and n = 11 the code takes two primes for answers of 25 and 33 bits; at n = 12 and 13 it takes three for 40 and 44. Every one of those is one modulus more than the answer needed, and the extra modulus is a full elimination over a finite field. There is no way to spend less without knowing the answer first, which is the position every a-priori bound in this collection puts its user in.

On the random family with one-digit entries the bound at n = 12 is 52 bits and the determinant occupies 40, so three primes are budgeted where two would have done. On the wide family with three-digit entries the bound at n = 12 is 133 bits and the answer 121: six primes budgeted, five needed. The waste is small and it is systematic — the bound is a product of row norms and the determinant is a signed sum with cancellation in it, so the bound exceeds the answer by whatever the cancellation was.

Why a loose bound would still be the right thing to use #

There is an obvious objection: a bound that is not tight is a bound that makes the routine do more work than it needs to. It is worth answering carefully, because the answer is not “the waste is small”.

The alternative to budgeting is checking, and checking is not available. A routine that computed residues until the reconstruction stopped changing would be reading a coincidence as a convergence. A routine that computed one extra prime and compared would detect a disagreement rather than a sufficiency — which is a genuinely useful thing to do, and it is what the next essay is about, but it tests for a different failure.

So the budget is not a conservative choice among several. It is the only mechanism available, and its being an over-estimate by a few per cent is the price of the whole route working at all.

It is worth putting that beside the site’s usual complaint about bounds. The growth factor in Gaussian elimination has a bound of 2ⁿ⁻¹ that no matrix anybody has ever assembled comes near, and the finding there is that the bound is useless for deciding anything. The a-priori bound in model reduction is the other extreme: attained exactly, at every order, to four decimals. Hadamard’s sits between them and is used differently from either — not to predict what will happen, and not to reassure, but as an input to the algorithm. A routine that budgets against a bound is a routine whose correctness depends on the bound being true rather than on its being tight.

And the bound is attained #

A bound that could never be reached would be a different objection, and this one can be.

Hadamard’s inequality holds with equality exactly when the rows are mutually orthogonal. A Hadamard matrix — entries ±1, rows orthogonal — attains it: at order eight, every row has norm √8, the product of the norms is 8⁴ = 4096, and the determinant is 4096. Twenty-five bits budgeted, twenty-five bits used, nothing wasted.

That is the refusal this essay carries. A claim that the bound is loose by more than a byte on every matrix is fed the Sylvester construction at order eight and has to reject it — and the assertion is run rather than described, so the claim in the paragraph above is one the build would catch if it stopped being true.

The families where it is enormously slack are equally instructive. The unit triangular family has a determinant of exactly 1 and rows of growing length, so the bound climbs linearly while the answer stays at one bit. The Pascal family is worse: entries reaching twenty bits, a bound of 149 bits at n = 12, and a determinant of 1. Both are correct behaviours of the inequality — orthogonal rows are what make it tight, and neither family has any — and both are the standing reminder that the size of an answer is not predicted by the size of the entries.

What this buys, measured against the alternatives #

The modular route’s advantage is not speed at these sizes. It is that the cost is a known number of word-sized eliminations rather than one elimination on operands of unknown length.

Fraction-free elimination performs n³/3 operations on integers that grow to about the determinant’s length — so the real cost is n³/3 multiplied by the cost of arithmetic at that length, which is quadratic in it for schoolbook multiplication. The modular route performs k·n³/3 operations on machine words, where k is the number of primes, and k is itself about the determinant’s length over 26. So the two are the same order and the constants are entirely different, and the second’s inner loop is a loop over doubles.

The comparison is worth making precisely, because the two routes are the same order and it is easy to conclude from that that they are the same. It is also easy to make it wrong, and the obvious version of it is wrong. Fraction-free elimination performs its n³/3 operations on operands that reach the determinant’s length L, at a cost quadratic in L for schoolbook arithmetic, so the total is about n³L²/3. The modular route performs k = L/26 eliminations of n³/3 word operations each, so the total is about n³L/78. Dividing one by the other gives 26L, which is the field’s usual answer to everything — the length of the numbers is the cost, and a route that never lets them get long never pays it. It is also about fifty times too large, and the section began by saying so.

The two expressions are not in the same unit #

The division above is not a comparison. The n³L²/3 counts bit operations, because L² is the schoolbook cost of multiplying two L-bit numbers measured in bit operations. The n³L/78 counts word operations, because the whole point of the modular route is that its inner loop runs on machine doubles. Dividing a bit count by a word count produces a number whose unit is bits per word, and 26 of the factor is exactly that: the word size, smuggled in as though it were an advantage.

Put both in word operations. Fraction-free elimination performs n³/3 operations on operands of L/26 words, at a schoolbook cost of (L/26)² word multiplications each, for n³/3 × (L/26)² in total. The modular route performs L/26 eliminations of n³/3 word operations, for (L/26) × n³/3. The ratio is L/26, which is not the length of the determinant but the number of primes — the same k the rest of the essay is about. That is a much better answer than 26L for a reason beyond the arithmetic: the modular route’s whole structure is k independent copies of one word-sized elimination, so a ratio that is k says the route costs one elimination against fraction-free’s k, which is a sentence about the algorithm rather than about a unit conversion.

What the ratio actually is, counted #

Even L/26 is an overestimate, because it charges every fraction-free operation at the final length. The Bareiss intermediates do not start there. At step k of the elimination the entries carry roughly the length of a k×k minor, so the operands grow from one word to L/26 words across the whole sweep, and the operations that dominate the count by number are the early ones on short operands. The only honest way to settle it is to count the words as they actually occur.

assertTheModularAdvantageIsAFactorOfTwoNotL does that: it runs the Bareiss elimination and, at every one of its multiplications and divisions, adds the product of its two operands’ lengths in 26-bit words. Against k·n³/3 for the modular route, on random matrices with two-digit entries:

A factor of about two. Not L, which runs from 37 to 112 down that column, and not quite L/26 either, which is the right shape but still an overcount by the growth argument above. The two rows that look out of order — n = 8 at 1.66 and n = 12 at 1.95, both below their neighbours — are the number of primes stepping up one before the determinant’s length has caught up with it, which makes the modular route pay for a prime it barely needs. That is the granularity the essay’s own bound section is about: k is a ceiling, and a ceiling is loose by up to a whole prime.

The ratio does grow, from 2.00 at n = 6 to 2.32 at n = 16, and it grows more slowly than L does over the same range. So the modular route does win, and it wins by more on a larger matrix, and both of those are worth saying — but the win at these sizes is a factor of two on a well-implemented schoolbook baseline, and any real implementation’s constant factors are larger than that. Which is precisely the claim this section opened with, and which the quoted ratio of L had spent a paragraph contradicting: the modular route’s advantage here is not speed.

That shape has a name by now. A formula is written down, it is arithmetically fine as a formula, and it is quoted as a measurement of something nobody measured — the same failure as quoting a tolerance the arithmetic never reaches or reporting a bound every answer satisfies. The tell here was internal: two paragraphs of one section disagreed by a factor of fifty, and the prose one was right. The prime that actually divides the answer is a number to be found; the ratio of two routes is a number to be counted; neither is a number to be derived and then trusted. An essay that is exact about a measured problem has no business being approximate about its own cost.

There is a second advantage that matters more in practice and is invisible in an operation count: the k eliminations are completely independent. Nothing in the modulo-p₁ computation is needed by the modulo-p₂ one. The essays about what a machine actually costs spend their time on data movement and dependence, and here there is none of either — which is a different kind of statement from a flop count and is the reason this route is the one that scales.

Where the residues are actually wrong #

One thing has been assumed throughout and it is not always true: that the elimination modulo p succeeds. It can fail, and it fails in a way that looks exactly like a correct answer.

If p divides det A then the matrix is singular over the field with p elements. The elimination reaches a column of zeros, reports a determinant of zero, and is right: that determinant genuinely is zero. What it is not is a residue of the determinant that was wanted — or rather it is, and it is the uninformative one, and there is no way to tell from inside a single residue whether the answer is zero because the matrix is singular or because the prime was unfortunate.

The same substitution reaches the rank. Rank over the field with p elements is a perfectly good quantity, computed by the same elimination with no threshold anywhere, and it is not the rank over the rationals — a second and quite separate way for an exact answer to answer a question nobody asked. Rank as a threshold decision is the floating-point version of the same trouble reached from the opposite direction, and the two have nothing in common but the symptom.

That is the whole cost of having no growth, and it is the next essay. It is a failure mode with no analogue in floating point: nothing in the arithmetic is inaccurate, no quantity is small enough to be mistaken for zero, and the essay about deciding that a zero has arrived is about a judgement that here is not being made. The zero is exact. It is the question that was substituted.

Choosing the primes, which is not arbitrary #

Three properties are wanted of the moduli and only one of them is obvious.

They must be pairwise coprime, which for distinct primes is free. Using a composite modulus is the failure this reaches for: the modular inverse of a nonzero entry may not exist, the elimination reports a singular matrix, and the report is about the modulus rather than about the matrix. A list of primes typed by hand is a place for a composite to hide, so the list here is generated by trial division rather than written down.

They must be small enough that a product of two residues is exact in the arithmetic doing the work. A double holds integers exactly to 2⁵³, so a modulus below 2²⁶·⁵ keeps every product inside it. That is the whole reason for the twenty-six: it is not a word size, it is half a significand. The same routine written for 64-bit integers would use 31-bit primes and the same argument.

And they must be large enough that not many are needed, which is the counter-pressure. Every prime is a full elimination, so halving the modulus doubles the work. Twenty-six bits is where those two meet on a machine whose fast arithmetic is floating point — which is a slightly startling place for an exact algorithm to end up, and it is the same observation the low-precision essays make from the other direction: the fastest arithmetic on a modern machine is the one with the shortest operands, and an algorithm that wants speed goes to meet it.

What is left over #

Two things about this route are not settled by the budget and both get their own essay.

The first is the unlucky prime: a modulus dividing the determinant, at which the elimination is correct and uninformative. The defence is to compute more residues than the bound requires and check that the extra ones agree — which is a genuinely different test from the sufficiency the budget provides, and which turns an undetectable wrong answer into a detected one.

The second is that everything above is about a determinant, which is a single integer. A solution is a vector of fractions, and the Chinese remainder theorem reconstructs integers. Getting a fraction back out of a residue is a separate mechanism with a separate bound, and the bound is about a lattice rather than about a product of primes.