For a prime power $q$ the group ${\rm GL}_2({\bf F}_q)$ has $(q^2-1)(q^2-q)$ elements. This happens to be a triangular number for $q=2$ (being $6 = 1+2+3$), and $-$ more notably, especially this year :-) $-$ for $q=7$ (being $1+2+3+\cdots+63 = 2016$). Indeed $2016$ is surely the largest triangular number of the form $(q^2-1)(q^2-q)$ for $q \in \bf Z$, even without the assumption that $q$ be a prime power. (The others are $300 = 1+2+3+\cdots+24$, for $q=-4$, and the trivial $0$ for $q=0$ and $q=\pm1$.)
As far as I know this question has no serious mathematical significance or application. But it seems surprisingly hard to prove that $2016$ is the last such number. I know how to do it, but not easily; this Diophantine problem might be just a bit beyond the available tools for an elementary proof, and thus a good test case for extending those tools.
How accessible / elementary a proof can we give for the fact that $-4,-1,0,1,2,7$ are the only integers $q$ for which $(q^2-1)(q^2-q)$ is triangular (equivalently: for which $8(q^2-1)(q^2-q)+1$ is a square)?
Here's what I know about this.
For starters, there are only finitely many such $q$ by Siegel's theorem on integral points, because the equation $(q^2-1)(q^2-q) = (s^2-s)/2$ gives a curve $C$ of positive genus $1$. But not only is Siegel's theorem far from elementary, it's also ineffective: one needs several large solutions to prove that no further solution exists, and usually there isn't even one solution that's large enough. (Or is there a proof of Siegel's theorem for which the $q = -4$ and $q=7$ solutions of our equation are enough to prove that the list of solutions is complete?)
At this point we're tempted to throw $C$ into a standard routine for finding integral points on an elliptic curve via Baker's effective bounds. But those bounds are even farther from elementary than Siegel's theorem. Moreover, they do not immediately apply here, because our curve $C$ is not in Weierstrass form $y^2 = x^3 + ax^2 + bx + c$.
Now Magma does have a function IntegralQuarticPoints for finding all integer points on a curve $y^2 = Q(x) = ax^4+bx^3+cx^2+dx+e$ given also a rational point. Here we have $x=q$ and $Q(x) = 8(q^2-1)(q^2-q)+1$, and we already know some rational points such as $(x,y)=(0,\pm 1)$; and indeed the input
time IntegralQuarticPoints([8,-8,-8,8,1],[0,1])
yields output ending with
[
[ 7, 127 ],
[ 2, -7 ],
[ -4, 49 ],
[ 0, -1 ],
[ -1, 1 ],
[ 1, 1 ]
]
Time: 0.460
(the Time is reported in seconds). But I don't know what Magma did to check that there are no other $q \in \bf Z$, and thus don't know how accessible the resulting proof might be.
Since $C$ is a genus-$1$ curve with rational points, we can transform $C$ birationally into Weierstrass form; for example, starting from the rational point $(q,s)=(0,0)$, where $s = -2q + O(q^2)$, we may write $s = -2q + (x+1)q^2$ to obtain $$ (x^2+2x-1)q^2 - (4x+2) q + 5-x = 0, $$ and this quadratic in $q$ has discriminant $4(x^3+x^2-7x+6)$ $-$ so we have identified $C$ with the curve $$E : y^2 = x^3 + x^2 - 7x + 6.$$ But this transformation does not take integral points of $E$ to integral points of $C$. (Geometrically, integrality depends on both the curve and the point(s) at infinity, and here we have the same curve but different distinguished points.)
We can still use this Weierstrass form $E$ to compute the rational points of $C$. If the group of rational points were finite we would be done; if it had rank $1$ we could usually, with more effort, use $p$-adic techniques to find all the integral points on $C$. (It may be a stretch to call these techniques elementary $-$ and even computing the rank of $E$ requires at least the arithmetic of a cubic number field, because the polynomial $x^3 + x^2 - 7x + 6$ is irreducible $-$ but this still feels like a much smaller hammer than the Baker bounds.) Unfortunately mwrank finds that $E$ has rank $2$, with generators $(x,y) = (1,1)$ and $(-3,3)$. Magma reported this too en route to computing its list of integral points, which must be why the input had to include some rational point of $C$. (This curve has conductor $2^3 331 = 2648$, and is listed in the LMDFB as elliptic curve 2648.a1.)
That leaves me with two backup options:
First, we can apply Baker over the quadratic field $K = {\bf Q}(\sqrt 2)$. The two points at infinity of $C$ are defined over $K$, so we can choose one of them as the origin of the group law to identify $C$ with $E$ over $K$, and then integral points of $C$ do go to integral points in $E(K)$. [There can be yet more integral points, because we're ignoring the other point at infinity of $C$; but those spurious can be eliminated at the end, or ignored using the condition that our point in $E(K)$ must come from $C({\bf Q})$. The group $E(K)$ has rank $3$, since it contains with finite index the direct sum of $E({\bf Q})$ and $E_2({\bf Q})$ where $E_2$ is the quadratic twist $2y^2 = x^3 + x^2 - 7x + 6$, isomorphic with the rank-$1$ LMFDB curve 21184.d1 of conductor $2^6 331$.]
Second $-$ and possibly requiring only a smaller "hammer" $-$ we can reduce to some quartic Thue equations over $K$. In the course of computing $E({\bf Q})$ and the relevant part of $E(K)$, we usually do a $2$-descent, and thus find a finite list of rational functions $f_i = N_i/D_i \in K(X)$ of degree $4$ such that every $(x,y) \in C({\bf Q})$ has $x = f_i(x')$ for some $i$ and $x' \in K$. At this point we don't really need to finish computing $E(K)$ by searching for rational points, because it's enough to find all integral values of each $f_i$ $-$ and this amounts to solving a finite list of quartic Thue equations $\Delta_i(t,u) = d$, where $\Delta_i$ is a homogeneous quartic corresponding to $D_i$, and the variables $t,u$ and target value $d$ are in $O_K = {\bf Z}[\sqrt 2]$. With luck all of those equations might be solvable with Skolem's $p$-adic method instead of invoking Baker. But that is still a long way from elementary, and will take orders of magnitude more than $0.460$ seconds...
Is there a simpler proof?
