How a question about matrix powers turns out to be solvable in polynomial time — even though the same question about a single matrix entry has been open for decades.
Take a square matrix $A$ with rational entries and start computing its powers: $A, A^2, A^3, \dots$ The entries jump around — some positive, some negative — but for certain matrices, something curious happens: after some point, every entry of every power is nonnegative, and stays that way forever. The matrix might be full of negative numbers, its early powers might flip signs chaotically, and then, past some threshold, all that turbulence dies out.
Call such a matrix eventually nonnegative. This post is about the following question:
Given $A$, can a computer decide whether it is eventually nonnegative? How fast?
The answer turns out to be yes, in polynomial time — and the road there passes through coin problems, resonating eigenvalues, a rigidity theorem from the theory of weighted automata, and a Galois-theoretic trick for taming numbers that are too big to write down. Everything you need is a solid course in linear algebra and a first course in abstract algebra. Let’s go.
A puzzle about matrix powers
Here is a concrete matrix to keep in mind throughout:
$$ A_\star=\frac{1}{30} \begin{pmatrix} 21 & -9 & 21\\ -9 & 21 & 21\\ 21 & 21 & -9 \end{pmatrix}. $$
Three of its nine entries are negative. Compute $A_\star^2$: all entries positive. Compute $A_\star^3$: negative entries again. This flip-flopping continues for a while — even powers look fine, odd powers don’t — and then, from the eighth power on, it stops for good:
Why? Because $A_\star$ is symmetric, we can diagonalize it. Its eigenvalues are $\tfrac{11}{10}$, $1$, and $-1$, with orthogonal projections $P_\rho, P_+, P_-$ onto the corresponding eigenlines, and
$$ A_\star^n=\left(\tfrac{11}{10}\right)^n P_\rho + P_+ + (-1)^n P_-. $$
The $(-1)^n$ term is what makes odd and even powers behave differently. But its size is fixed, while the $\left(\tfrac{11}{10}\right)^n$ term grows. A short computation shows the smallest entry of $A_\star^n$ is $\bigl((\tfrac{11}{10})^n-1\bigr)/3$ for even $n$ and $\bigl((\tfrac{11}{10})^n-2\bigr)/3$ for odd $n$; once $(\tfrac{11}{10})^n>2$ — that is, once $n\geq 8$ — every entry is strictly positive, forever.
So this instance we can settle by hand. The eigenvalue $\tfrac{11}{10}$ eventually drowns out the oscillation between $1$ and $-1$. The general question is whether a computer can always settle it: for every rational matrix, decide whether the sign turbulence eventually dies out.
The surprise: one entry is an open problem
This problem has a natural variant that sounds simpler but is actually significantly more complicated.
Fix a single position $(i,j)$ and look at the sequence of numbers $u_n = (A^n)_{ij}$. By the Cayley–Hamilton theorem, $A$ satisfies its own characteristic polynomial, and therefore so does every entry sequence: $u_n$ satisfies a linear recurrence whose coefficients come from the characteristic polynomial of $A$. For example, with
$$ A=\begin{pmatrix}1&1\\1&0\end{pmatrix}, \qquad A^n=\begin{pmatrix}F_{n+1}&F_n\\F_n&F_{n-1}\end{pmatrix}, $$
the entries are Fibonacci numbers, satisfying $u_{n+2}=u_{n+1}+u_n$. The correspondence goes both ways: every linear recurrence sequence (LRS) over $\mathbb{Q}$ shows up, up to a harmless index shift, as an entry of the powers of some rational matrix. Matrix entries and linear recurrences are the same objects in different clothing.
Deciding whether a given LRS is eventually nonnegative is a famous open problem — the Ultimate Positivity Problem. Despite sustained attacks, nobody knows an algorithm that decides it for all linear recurrences; positive results exist only for special classes (low order, or no repeated eigenvalues, and similar restrictions). It is one of those problems that looks like it should have been settled in the 1930s and simply wasn’t.
So we are in the following bizarre situation. A $d \times d$ matrix has $d^2$ entry sequences. Asking whether one prescribed entry sequence is eventually nonnegative: open. Asking whether all $d^2$ of them simultaneously are eventually nonnegative: decidable in polynomial time — as this post will explain.
The resolution of the paradox is that the $d^2$ sequences are not independent — they are cut from a single multiplicative orbit $A, A^2, A^3,\dots$ — and demanding that all of them behave forces an enormous amount of hidden algebraic structure. A single sequence can misbehave in ways that a full matrix of sequences cannot. Most of this post is about uncovering that structure.
Why single sequences are hard: rotations and Diophantine trouble
To appreciate the difficulty, let’s look at the complex unit circle.
Take the rotation-like matrix
$$ R=\frac15\begin{pmatrix}3&-4\\4&3\end{pmatrix}, \qquad R^n=\begin{pmatrix}\cos n\theta&-\sin n\theta\\\sin n\theta&\cos n\theta\end{pmatrix}, $$
where $\theta$ is the angle with $\cos\theta=\tfrac35$, $\sin\theta=\tfrac45$. Its eigenvalues are $z=\tfrac{3+4i}{5}$ and $\bar z$, sitting exactly on the unit circle.
Is $\theta$ a rational multiple of $\pi$? No — and there is a lovely classical reason. Niven’s theorem says that if $\theta$ is a rational multiple of $\pi$ and $\cos\theta$ is rational, then $\cos\theta\in\{0,\pm\tfrac12,\pm1\}$. Since $\cos\theta=\tfrac35$ is none of these, $\theta/\pi$ is irrational. Equivalently: $z$ is a point on the unit circle that is not a root of unity, so the sequence $z, z^2, z^3, \dots$ never cycles; by Weyl’s equidistribution theorem it wanders around the circle forever, filling it densely.
Consequently every entry of $R^n$ — say $\cos n\theta$ — changes sign infinitely often, in a pattern governed by how the multiples of an irrational angle fall around the circle. Questions like “is $u_n \geq 0$ for all large $n$?” for sequences built from such terms hinge on how well $\theta/\pi$ can be approximated by rationals — a Diophantine approximation problem. For specific constants, such approximation questions are notoriously out of reach, and this is precisely the wall the Ultimate Positivity Problem runs into: a sequence like $u_n = 1-\cos n\theta$ is nonnegative forever, but a tiny perturbation of it dips below zero depending on number-theoretic facts nobody knows how to compute.
Notice also what $R$ itself does as a matrix: no power of $R$ is ever entrywise nonnegative. The off-diagonal entries are $-\sin n\theta$ and $+\sin n\theta$; for both to be $\geq0$ we would need $\sin n\theta = 0$, impossible since $\theta/\pi$ is irrational. The matrix question dodges the Diophantine trap here — the answer is just “no”, and for a robust structural reason. That is a first glimpse of the rigidity we will exploit.
The general shape of the enemy is now visible. Write down the closed form of an entry sequence (we assume $A$ is diagonalizable; repeated eigenvalues add polynomial-in-$n$ factors, just like resonant terms in differential equations):
$$ (A^n)_{ij}=\sum_{\lambda} c_\lambda\, \lambda^n , $$
a sum over eigenvalues with fixed coefficients. Trouble comes from eigenvalues of equal largest modulus whose contributions oscillate against each other — like $1$ and $-1$ did for $A_\star$, or $z$ and $\bar z$ for $R$ — and the worst trouble comes from ratios $\lambda/\gamma$ that are not roots of unity, where the oscillation never settles into a cycle.
The proof strategy has three main ideas. We’ll look at them one at a time.
Idea 1: The chicken nugget reduction
Consider the set of “good exponents”
$$ S=\{n\geq1 : A^n\geq0\}. $$
Here’s the key point, which is surprisingly simple: the product of two nonnegative matrices is nonnegative. So if $A^m\geq0$ and $A^n\geq0$, then $A^{m+n}=A^mA^n\geq0$. The set $S$ is closed under addition — it is a sub-semigroup of the natural numbers.
Additively closed sets of integers are the subject of a classic piece of mathematical folklore, the Frobenius coin problem (also known as the chicken nugget theorem): if a set of natural numbers is closed under addition and contains two coprime elements $a$ and $b$, then it contains every positive integer greater than $ab-a-b$. (With nuggets sold in boxes of $a=9$ and $b=13$: every order size beyond $9\cdot13-9-13=95$ is exactly purchasable.)
Now suppose we could show, for some modulus $D$, that
$$ A^{Dk+1}\geq0 \quad\text{for all large } k. $$
Then $S$ contains the two consecutive terms $a=DK+1$ and $b=D(K+1)+1$ of that progression, and
$$ \gcd(a,b)=\gcd(DK+1,\,D)=\gcd(1,D)=1 $$
— consecutive terms of the progression $Dk+1$ are automatically coprime, because their difference is $D$ and each is $\equiv 1 \pmod D$. The coin theorem then fills in every sufficiently large exponent. Conclusion:
$A$ is eventually nonnegative $\iff$ $A^{Dk+1}\geq0$ for all sufficiently large $k$ — for any modulus $D$ we like.
This matters because, as we will see, the natural analysis of matrix powers involves splitting the exponents into residue classes modulo a certain integer $D$ — and that $D$, while it has polynomially many digits, can be exponentially large as a number. Checking all $D$ residue classes one by one is hopeless. The semigroup trick means that we don’t need to.
You can watch it work on $A_\star$: we will see that its natural modulus is $D=4$, and indeed $A_\star^{4k+1}\geq0$ from $k=2$ onward. The exponents $9$ and $13$ are coprime, so nonnegativity propagates to every exponent beyond $95$. (The true threshold is $8$ — the coin bound is deliberately crude, but crude is fine.)
One warning: this argument is genuinely about matrices. The scalar analogue — “if $u_{Dk+1}\geq0$ for large $k$ then $u_n \geq 0$ for large $n$” — is false: the sequence $u_n=1-2(-1)^n$ equals $3$ at every odd index and $-1$ at every even index. The set $\{n : u_n\geq0\}$ of a scalar sequence is not closed under addition; the good set of a matrix is. This is a foretaste of why the matrix version of the problem is easier than the scalar version.
Idea 2: Merging the resonant eigenvalues
Which modulus $D$ should we choose? The one that kills all the rational resonances among eigenvalues.
Call two nonzero eigenvalues $\lambda,\gamma$ torsion-equivalent if $\lambda/\gamma$ is a root of unity. (The name is from group theory: the roots of unity are exactly the torsion elements — elements of finite order — of the multiplicative group $\mathbb{C}^\times$.) This chops the eigenvalues into torsion classes. Within one class, eigenvalues differ only by roots of unity: they have the same modulus, and their powers cycle around each other with a common period.
Take, for instance, a matrix with eigenvalues $\{2,\,-2,\,3i,\,-3i,\,5\}$. Then $-2/2=-1$ and $3i/(-3i)=-1$ are roots of unity, while any ratio between the groups $\{2,-2\}$, $\{3i,-3i\}$, $\{5\}$ has modulus $\neq1$ and can’t be one. Three torsion classes. Every order that appeared was $2$, so let $D=4$ — twice the least common multiple of the orders, the factor of two being a convenience that makes real numbers’ $D$-th powers positive. Now raise everything to the $D$-th power:
$$ 2^4=(-2)^4=16,\qquad (3i)^4=(-3i)^4=81,\qquad 5^4=625. $$
Each class collapses to a single positive number — its class root — and the three class roots have no root-of-unity ratios left among them:
This is where the residue classes come from. Along the progression $n=Dk+1$, the closed form of an entry reorganizes itself beautifully. Each eigenvalue contributes $c_\lambda\lambda^{Dk+1} = c_\lambda\lambda\cdot(\lambda^D)^k$, and eigenvalues in the same torsion class share the base $\mu_C := \lambda^D$. Grouping them:
$$ (A^{Dk+1})_{ij}\;=\;\sum_{C}\; b_C^{ij}\,\mu_C^{\,k}, \qquad b_C^{ij}\;=\;\sum_{\lambda\in C} \lambda\,c_\lambda^{ij}, $$
a sum over torsion classes, with one base $\mu_C$ per class and one grouped coefficient $b_C^{ij}$ per class. (If $A$ has repeated eigenvalues with nontrivial Jordan structure, the coefficients become polynomials in $k$ rather than constants — the statements below survive verbatim with “leading coefficient” in place of “coefficient”, so we’ll suppress this.) Two facts make this expansion the right object:
- No resonances remain. Distinct class roots never have a root-of-unity ratio — we chose $D$ precisely to guarantee it. In LRS language, the sequence $k \mapsto (A^{Dk+1})_{ij}$ is nondegenerate, the well-behaved kind of recurrence.
- The expansion is unique. A standard lemma says an exponential sum with distinct bases and nonzero coefficients represents a sequence in only one way. So the numbers $b_C^{ij}$ are intrinsic to the matrix — they are the quantities our algorithm should examine.
For $A_\star$: eigenvalues $\tfrac{11}{10}, 1, -1$; the quotient $1/(-1)=-1$ has order $2$, and no quotient involving $\tfrac{11}{10}$ has modulus $1$; classes $\{\tfrac{11}{10}\}$ and $\{1,-1\}$; modulus $D=4$. Along $n=4k+1$:
$$ A_\star^{4k+1} =\underbrace{\tfrac{11}{10}\left(\tfrac{14641}{10000}\right)^{k}P_\rho}_{\text{class }\{11/10\}} +\underbrace{(P_+-P_-)}_{\text{class }\{1,-1\},\ \mu_C=1}. $$
The oscillating pair $1,-1$ has merged into a single constant term with base $1$. All sign turbulence within a class disappears along the progression; what remains is a clean competition between classes.
Idea 3: Nonnegativity is rigid
Along our progression, each entry is a tame exponential sum $\sum_C b_C^{ij}\mu_C^k$. When is such a sum eventually nonnegative? The obvious sufficient condition: if among the classes with $b_C^{ij}\neq0$ (active classes) there is a unique one of largest modulus, and its coefficient is positive, then that term eventually dominates everything else and the entry is eventually positive. If no class is active at all, the entry is eventually zero. Either way, fine.
The dangerous case is a tie: two active classes of equal largest modulus, whose contributions could cancel each other in perpetuity, nonnegativity hanging on a knife’s edge. This is exactly the configuration that makes Ultimate Positivity impenetrable for a single sequence.
The third idea is that for a whole eventually nonnegative matrix, ties simply cannot happen. The tool is a rigidity theorem with a Perron–Frobenius flavor.
The Perron–Frobenius theorem says that an entrywise nonnegative matrix $M$ has its spectral radius $\rho$ as an actual eigenvalue. Also, its eigenvalues of maximal modulus are exactly of the form $\rho\zeta$ with $\zeta$ a root of unity — the peripheral spectrum is a scaled cyclic pattern, never an irrational-angle mess. Nonnegativity of the matrix disciplines the spectrum. A theorem of Berstel from the theory of rational series pushes this to sequences: if a sequence can be realized with nonnegative data, as
$$ v_m=x^\top M^m y \qquad\text{with } x, M, y \text{ entrywise nonnegative}, $$
and is nondegenerate and not eventually zero, then its closed form has a unique base of maximum modulus.
Where does this theorem come from? Some intuition, in an optional aside — click to expand.
Two remarks on what this theorem does and does not say. It does not say nonnegative systems cannot oscillate — they can, but only periodically. The permutation matrix $\begin{pmatrix}0&1\\1&0\end{pmatrix}$ realizes the sequence $1,0,1,0,\dots$ with entirely nonnegative data, and that sequence has two equally large dominant bases, $1$ and $-1$. What saves the theorem is its nondegeneracy hypothesis: $-1/1$ is a root of unity, so this sequence is degenerate — and root-of-unity ties are exactly what our choice of $D$ eliminated in advance. Berstel and nondegeneracy split the work: nonnegativity forces dominant roots to differ only by roots of unity, and nondegeneracy says no two of ours do.
What nonnegativity truly forbids is aperiodic oscillation — a dominant root $\rho e^{i\theta}$ with $\theta/\pi$ irrational. Here is the mechanism, and it is worth internalizing. Think of a nonnegative $M$ as a weighted graph on states. If the graph mixes — some power $M^k$ is strictly positive — then that power averages: it sends every nonzero nonnegative vector to a strictly positive one, and iterated averaging collapses all directions onto a single dominant ray, the Perron eigenvector. An undamped rotation at the top scale would be a direction that refuses to collapse; contradiction. So a mixing nonnegative system has no peripheral oscillation at all. The only escape from mixing is combinatorial: the zero pattern of $M$ may split the states into classes $1\to2\to\cdots\to h\to1$ that the dynamics cycles through like a clock, with mixing inside each class. Then all top-scale oscillation is the clock itself — eigenvalues $\rho\zeta$ with $\zeta^h=1$ — and a finite clock can only tick at rational angles. A nonnegative machine’s only way to remember phase is “which class am I in,” and that memory is finite. Sustaining an irrational rotation would require amplitudes that interfere destructively at a never-repeating phase — negative or complex weights. That is exactly what our rotation matrix $R$ used, and notice it had negative entries; the irrational-angle mess is a luxury of mixed signs.
Now watch the matrix structure hand us exactly such realizations. Suppose $A$ is eventually nonnegative, and pick $t$ large enough that $M_0=A^{Dt}$ and $M_1=A^{Dt+1}$ are both nonnegative. Sample the entry sequence along a sub-progression:
$$ v_m=\bigl(A^{Dt+1}\cdot (A^{Dt})^m\bigr)_{ij}=e_i^\top M_1M_0^m\,e_j $$
— nonnegative vector, times powers of a nonnegative matrix, times a nonnegative vector. Berstel’s theorem applies, forces a unique dominant base for $v$, and (after checking that sampling merges no classes — this is where nondegeneracy earns its keep) a unique dominant class for the entry itself. A limiting argument then pins down the sign: the dominant coefficient must be positive. We have arrived at an exact characterization:
Criterion. $A$ is eventually nonnegative $\iff$ for every entry $(i,j)$: either no class is active (the entry sequence is eventually zero), or there is a unique dominant active class and its grouped coefficient satisfies $b^{ij}>0$.
At this point the mathematical characterization is complete, and with it decidability. What is not at all clear is that any of it can be computed efficiently. That is the second half of the story, and it’s where the classical Galois theory gets involved.
The computational cliff
To run the criterion we must, for each entry, (a) find the torsion classes, (b) compute the grouped coefficients $b_C^{ij}$, (c) find the dominant active class, (d) test the sign of its coefficient. Each step conceals a trap of exponential size.
Trap 1: $D$ is astronomically large. Each root-of-unity order among eigenvalue ratios is at most polynomial (an order-$r$ root of unity has algebraic degree $\varphi(r)\geq\sqrt{r/2}$, and ratios of eigenvalues of a $d\times d$ matrix have degree at most $d^2$, so $r\leq 2d^4$). But $D$ is built from the least common multiple of many such orders and can genuinely have magnitude exponential in $d$. So we can never afford to compute $A^D$, nor any $\lambda^D$: the number $\lambda^D$ has about $D$ digits. The class roots $\mu_C=\lambda^D$ — the very bases of our expansion — are too large to ever write down. The algorithm must reason about them purely symbolically, comparing $|\mu_C|$’s by comparing $|\lambda|$’s of representatives, which is fine since torsion-equivalent eigenvalues share a modulus.
Trap 2: the splitting field is huge. The coefficient $b_C^{ij}=\sum_{\lambda\in C}\lambda\,c^{ij}_\lambda$ sums contributions of different eigenvalues, each living in its own number field $\mathbb{Q}(\lambda)$ of degree at most $d$. To add them “properly” you would work in a field containing all of them — and the splitting field of a degree-$d$ polynomial can have degree $d!$ (the Galois group can be the full symmetric group $S_d$, shuffling the roots every which way). For $d=20$, that’s a field of degree $2{,}432{,}902{,}008{,}176{,}640{,}000$. Exact arithmetic there is not happening.
Each individual contribution is harmless: eigenvalues are roots of the factors of the characteristic polynomial, each factor has degree $\leq d$, and classical polynomial-time algorithms (polynomial factorization over $\mathbb{Q}$, linear algebra over $\mathbb{Q}(\lambda)$) compute each $\lambda\,c^{ij}_\lambda$ exactly, with polynomially many bits, inside its own small field. The problem is purely the sum across fields — deciding, say, whether contributions from $\mathbb{Q}(\sqrt[3]{2})$ and $\mathbb{Q}(\omega\sqrt[3]{2})$ cancel exactly, without ever building a field containing both.
Galois symmetry to the rescue
Recall from your algebra course: an algebraic number $\alpha$ has a minimal polynomial over $\mathbb{Q}$, and the other roots of that polynomial are $\alpha$’s conjugates — its algebraically indistinguishable siblings. ($\sqrt2$ and $-\sqrt2$; the three cube roots of $2$.) The degree of $\alpha$ equals its number of conjugates, and the conjugates are exactly the images of $\alpha$ under the field symmetries (automorphisms) of any large enough field. The one fact we need is a mirror of this:
If every symmetry $\sigma$ maps $\alpha$ into a fixed finite set $B$ (with $\alpha\in B$), then the degree of $\alpha$ is at most $|B|$ — because all of $\alpha$’s conjugates lie in $B$.
(In the extreme case $|B|=1$, this is the familiar fact that a number fixed by all symmetries is rational.)
Now look at a grouped coefficient $b_C=\sum_{\lambda\in C}\lambda\,c_\lambda$ through the eyes of a symmetry $\sigma$. Two observations:
- $\sigma$ permutes the eigenvalues of $A$ (they are roots of the rational polynomial $\chi_A$), and it respects torsion classes: if $\lambda/\gamma$ is a root of unity, so is $\sigma(\lambda)/\sigma(\gamma)$ (symmetries map roots of unity to roots of unity). So $\sigma$ doesn’t scramble eigenvalues arbitrarily — it permutes whole classes.
- The coefficient $c_\lambda$ is built from $\lambda$ and the (rational!) matrix $A$ by field operations — it is a “formula in $\lambda$”. Symmetries commute with formulas: $\sigma(c_\lambda)=c_{\sigma(\lambda)}$.
Put together: applying $\sigma$ to the class sum $b_C$ just transports it to another class sum,
$$ \sigma\Bigl(\sum_{\lambda\in C}\lambda\,c_\lambda\Bigr)=\sum_{\lambda\in C}\sigma(\lambda)c_{\sigma(\lambda)}=\sum_{\gamma\in \sigma(C)}\gamma\,c_\gamma=b_{\sigma(C)}. $$
Every conjugate of $b_C$ is itself one of the class sums — and there are at most $d$ classes! So:
Orbit bound. Each grouped coefficient $b_C^{ij}$ is an algebraic number of degree at most $d$.
Not $d!$, not $d^{\,\text{anything}}$ — just $d$. The individual summands may live in incompatible fields whose compositum is monstrous; their sum is constrained by symmetry to be tame enough to touch. We never find out which small field $b_C$ lives in, and we never need to: the degree bound alone, plus a routine (if fiddly) bookkeeping argument showing $b_C$ also has polynomially many bits’ worth of “height”, is all the leverage required.
Let’s watch the whole mechanism on one worked example. Take the $6\times6$ matrix (with blank blocks all zero)
$$ M=\begin{pmatrix} 0&0&0&4&&\\ 1&0&0&0&&\\ 0&1&0&0&&\\ 0&0&1&0&&\\ &&&&1&1\\ &&&&1&0 \end{pmatrix}. $$
The first block is the companion matrix of $X^4-4$: a matrix built to have exactly that characteristic polynomial. Its eigenvalues are the fourth roots of $4$, namely $\pm\sqrt2$ and $\pm i\sqrt2$; by Cayley–Hamilton, this block — call it $B$ — satisfies $B^4=4I$, which we will use below. The second block is the Fibonacci matrix from earlier, with eigenvalues the golden ratio $\varphi=(1+\sqrt5)/2\approx1.618$ and its conjugate $\psi=(1-\sqrt5)/2\approx-0.618$. The torsion classes are
$$ C_0=\{\sqrt2,\,-\sqrt2,\,i\sqrt2,\,-i\sqrt2\},\qquad C_\varphi=\{\varphi\},\qquad C_\psi=\{\psi\}. $$
The four eigenvalues of the first block merge into one class: their pairwise ratios are $\pm1$ and $\pm i$, all roots of unity. The golden pair stays split: $|\varphi/\psi|=\varphi^2\approx2.6$, and a ratio of modulus $\neq1$ cannot be a root of unity. Every other ratio also has modulus $\neq1$, so there is no further merging. The root-of-unity orders that appeared are $1$, $2$, and $4$, so the modulus is $D=2\cdot\mathrm{lcm}(1,2,4)=8$, and the three class roots $(\sqrt2)^8=16$, $\varphi^8$, $\psi^8$ have no root-of-unity ratios left.
Notice that torsion classes are not conjugacy families. These are two different relations between the eigenvalues, and here they cut across each other.
The class $C_0$ glues together two different conjugate pairs. The pair $\pm\sqrt2$ consists of the roots of $X^2-2$; the pair $\pm i\sqrt2$ consists of the roots of $X^2+2$. The two pairs share no minimal polynomial, so they are not conjugate. They do not even live in the same field: $\mathbb{Q}(\sqrt2)$ is real, $\mathbb{Q}(i\sqrt2)$ is not. Torsion merges them anyway, because their ratio $i$ is a root of unity.
In the other direction, $\varphi$ and $\psi$ are conjugates, roots of the same $X^2-X-1$, yet they land in different torsion classes, because they have different magnitudes.
Now the symmetries. The splitting field of all six eigenvalues is $\mathbb{Q}(\sqrt2,i,\sqrt5)$, and its symmetries flip the signs of $\sqrt2$, $i$, and $\sqrt5$ independently. Every symmetry maps $C_0$ to itself, shuffling its four members internally. The symmetry $\sqrt5\mapsto-\sqrt5$ swaps $C_\varphi$ and $C_\psi$ wholesale. Note what that means: it carries the dominant class ($\varphi\approx1.62$) to the subdominant one ($|\psi|\approx0.62$). Symmetries do not respect magnitude. That is why the orbit argument uses them only to bound degrees, and never to decide which class dominates.
The split conjugate family. The $(5,5)$ entry of $M^n$ is the Fibonacci number $F_{n+1}$ — we computed the powers of the Fibonacci matrix earlier. The Fibonacci numbers have a closed form, Binet’s formula:
$$ F_{n+1}=\frac{\varphi^{n+1}}{\sqrt5}-\frac{\psi^{n+1}}{\sqrt5}. $$
The class sums of this entry can be read off directly: substitute $n=8k+1$ and collect each eigenvalue’s leftover powers into its coefficient,
$$ F_{8k+2} =\underbrace{\frac{\varphi^{2}}{\sqrt5}}_{=\,b_{C_\varphi}}\bigl(\varphi^{8}\bigr)^{k} +\underbrace{\Bigl(-\frac{\psi^{2}}{\sqrt5}\Bigr)}_{=\,b_{C_\psi}}\bigl(\psi^{8}\bigr)^{k}. $$
In exact form, $b_{C_\varphi}=(5+3\sqrt5)/10\approx1.171$ and $b_{C_\psi}=(5-3\sqrt5)/10\approx-0.171$.
Now check the general theory against these two numbers. The symmetry $\sqrt5\mapsto-\sqrt5$ swaps the two classes, so it must swap the two class sums — and visibly it does. Each class sum is an irrational number of degree $2$; the orbit bound caps the degree at the number of classes, which is $3$. The two class sums add up to $1$: the rational entry $M_{55}$, split into two conjugate irrational pieces. Finally, the verdict. The dominant class is $C_\varphi$, and its coefficient is positive. So the criterion accepts this entry. The negative coefficient on the subdominant class does no harm.
The merged class, and the cancellation. Every symmetry maps $C_0$ to itself, so every symmetry fixes the class sums of $C_0$, and they must be rational numbers. The $(1,1)$ entry shows what this looks like in practice. Since $B^4=4I$, and since $B$, $B^2$ and $B^3$ all have $(1,1)$ entry zero, the sequence $(B^n)_{11}$ is easy to describe: it equals $4^{n/4}$ when $4$ divides $n$, and $0$ otherwise. Written as a sum over the four eigenvalues, this says
$$ (B^n)_{11}=\frac{(\sqrt2)^n}{4}+\frac{(-\sqrt2)^n}{4}+\frac{(i\sqrt2)^n}{4}+\frac{(-i\sqrt2)^n}{4} $$
— check it: the four $n$-th powers add up to $4\cdot4^{n/4}$ when $4$ divides $n$, and cancel to $0$ otherwise. So each eigenvalue $\theta\in C_0$ appears with coefficient $c_\theta=\tfrac14$, and contributes $\theta/4$ to the class sum. The four contributions are
$$ \frac{\sqrt2}{4},\qquad-\frac{\sqrt2}{4},\qquad\frac{i\sqrt2}{4},\qquad-\frac{i\sqrt2}{4}. $$
Each one is nonzero. The first two live in the field $\mathbb{Q}(\sqrt2)$; the last two live in the different field $\mathbb{Q}(i\sqrt2)$. And yet they cancel exactly: the class sum is $0$, and along the progression the entry is identically zero (consistent with the description above, since $4$ never divides $8k+1$).
This is the situation the orbit bound exists for. Here a single coefficient is a sum of terms from incompatible fields, so no common-field arithmetic can simplify it. The algorithm instead evaluates each term numerically in its own field, with certified error bars, and adds the results. It sees a tiny number, something like $0.0000000003$. It must then decide: is the coefficient exactly zero — an eventually-zero entry, the first branch of the criterion — or merely small, in which case its sign decides the entry?
Which brings us to the last tool.
Knowing when zero is zero
Numerical computation can never prove a quantity is exactly $0$ — unless you know in advance that the quantity, if it is nonzero at all, cannot be too small. For algebraic numbers of bounded degree and bounded description length, that is a classical fact.
Suppose $\alpha\neq0$ has degree $e$, with minimal polynomial $a_eX^e+\dots+a_0\in\mathbb{Z}[X]$. The product of all the roots is $\pm a_0/a_e$, and $a_0\neq0$ (else $X$ divides a minimal polynomial), so
$$ |\alpha|\cdot\prod_{\substack{\text{conjugates }\beta\neq\alpha}}|\beta| \;=\;\Bigl|\frac{a_0}{a_e}\Bigr|\;\geq\;\frac{1}{|a_e|}. $$
If each conjugate has absolute value at most $T$, then $|\alpha|\geq 1/(|a_e|T^{\,e-1})$. (We can easily bound $T$ in terms of the coefficients of the minimal polynomial.)
For our grouped coefficients we have exactly the needed inputs: degree $\leq d$ by the orbit bound, height polynomial by bookkeeping. So there is an explicit $\delta=2^{-\mathrm{poly}}$ with: either $b_C^{ij}=0$, or $|b_C^{ij}|\geq\delta$. Now we can safely use numerical calculation. Using ball arithmetic — every quantity carried as (center, error radius), with rigorously tracked error bars — evaluate each eigenvalue’s contribution in its own field, precisely enough that the total error of the sum stays below $\delta/2$, and add the balls. If the resulting ball excludes $0$, the coefficient is nonzero and we read off its sign; if it contains $0$, then since a nonzero coefficient of size at least $\delta$ would have kept the whole ball away from $0$, the coefficient is exactly zero. Polynomially many bits of precision suffice. Exact decisions, from finite-precision computation.
The algorithm, assembled
Everything now clicks together into five polynomial-time steps:
- Factor and isolate. Factor $\chi_A$ over $\mathbb{Q}$; represent each distinct eigenvalue by its minimal polynomial plus an isolating disc in $\mathbb{C}$.
- Find the classes. For each pair of eigenvalues, test whether $\lambda^r=\gamma^r$ for some $r\leq 2d^4$ (certified zero tests again); this reveals the torsion classes and the modulus $D$ — kept in binary, never as a value.
- Compute grouped coefficients. For each entry and class, evaluate each eigenvalue’s exact contribution inside its own $\mathbb{Q}(\lambda)$, sum with ball arithmetic, and decide zero/sign via the separation bound.
- Compare classes. Order the active classes by $|\lambda|$ of representatives (again a certified comparison of small algebraic numbers).
- Apply the criterion. Accept the entry if it is eventually zero, or if a unique dominant class survives with positive coefficient; accept $A$ iff every entry is accepted — correctness courtesy of the coin trick, the grouped expansion, and Berstel rigidity.
No $A^D$, no $\lambda^D$, no splitting field — the three unaffordable objects are precisely the three never constructed.
What’s still open
The strangest part of the story is the gap it leaves behind. We now know:
- Eventual nonnegativity of one prescribed entry: open (this is Ultimate Positivity).
- Eventual nonnegativity of all $d^2$ entries: polynomial time.
What about a prescribed subset of entries — a row? the diagonal? two corners? Nothing is known between the extremes. Somewhere between “one sequence” and “all sequences from one matrix orbit,” rigidity kicks in, and nobody knows where.
Two more directions beckon. For families of several commuting matrices, the analogous questions are only known to be decidable assuming Schanuel’s conjecture, a deep open problem in transcendence theory; whether the coin-theorem collapse has a multi-dimensional analogue is open. And while the algorithm certifies a threshold beyond which all powers are nonnegative, computing the least such threshold — or deciding “$A^n\geq0$ for all $n\geq N$” for a given binary $N$ — remains unresolved.
The moral: hard problems about sequences can become easy problems about the systems that generate them, because a system must be globally consistent in ways a single observable need not. Sometimes the way past an impossible question about one number is to ask it about all of them at once.