Diophantine equation for 2016: triangular $|{\rm GL}_2({\bf F}_q)|$ 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?

 A: This solution works iff $q$ is a prime power, so it's not a full solution, but I guess it's better than nothing, as it is completely elementary.
If $q$ is a prime power, then $q$ must fully divide either $s$ or $s-1$. We may WLOG assume it divides $s$ (since otherwise we may take $t=1-s$ as $s$). Then we have that, letting $s=mq$,
$$2(q-1)^2(q+1)=m(mq-1).$$
Taking $\bmod q$, we have that $m\equiv -2\bmod q$, so we may set $m=(x+1)q-2$, giving $s=-2q+(x+1)q^2$ and 
$$(x^2+2x-1)q^2 - (4x+2) q + (5-x) = 0,$$
as noted in the question. From this, we have $x\equiv 5\bmod q$. Let $x=nq+5$. Then
$$s=-2q+(nq+6)q^2=nq^3+6q^2-2q=6q^2+q(nq^2-2).$$

We will prove the following lemma:
Lemma: If $q$ is a nonzero integer and 
$$(2s-1)^2 = 8(q-1)^2q(q+1)+1,$$
then $|s|<3q^2.$
Proof: First, note that
$$3x^4-(x-1)^2x(x+1)=2x^4+x^3+x^2-x=(2x^2-x)(x^2+x+1)$$
is positive for all $x$ outside of $[0,1/2]$, and specifically is positive if $x$ is a nonzero integer. Thus,
$$(x-1)^2x(x+1)<3x^4.$$
$$(2s-1)^2=8(q-1)^2q(q+1)+1<24q^4+1\leq 25q^4.$$
So,
$$|2s-1|<5q^2.$$
$$|2s-1|+|1|<5q^2+1\leq 6q^2.$$
$$2|s|< 6q^2,$$
implying the result.

If $n=0$, then $s=-2q+6q^2$, and we have
$$(6q^2-2q)(6q^2-2q-1)=2(q-1)^2q(q+1).$$
This quartic has a triple root at $q=0$ and a single root at $q=11/17$, so we can safely ignore this case. Thus, $|n|\geq 1$. If $|q|\geq 11$, then $|nq|\geq 11$. So, 
$$|nq^2|\geq 11|q|.$$
$$|nq^2|\geq 9|q|+2.$$
$$|nq^2-2| \geq 9|q|.$$
But
$$s=6q^2+q(nq^2-2).$$
So, since $|nq^2-2|\geq 9|q|$, we have that $|s|\geq 3q^2,$ which contradicts our lemma. Thus, we have reduced the problem to a finite case check on integers $-10\leq q\leq 10$, which can easily be done.
