A curious identity related to finite fields To three elements $a_1$, $a_2$, $a_3$ in the finite field $\mathbb F_q$
of $q$ elements we associate the number $N(a_1,a_2,a_3)$
of elements $a_0\in \mathbb F_q$ such that the polynomial
$x^4+a_3x^3+a_2x^2+a_1x+a_0$ splits into four distinct linear factors over
$\mathbb F_q$. The number 
$$\sum_{(a_1,a_2,a_3)\in\mathbb F_q^3}{N(a_1,a_2,a_3)\choose 2}$$
counts then the number of minimal pairs in the Craig lattice $C_{q-1,3}$ 
(for $q$ a prime number). Experimentally (for all primes up to 2000)
this number seems to be given by
$$\frac{1}{1152}q(q-1)(q^3-21q^2+171q-c_q)$$
(the factors $q$ and $(q-1)$ are easy to explain through the action of the affine group) where 
$$c_q=\left\lbrace\begin{array}{ll}
455\qquad&q\equiv 5\pmod{24}\\
511&q\equiv 7\pmod{24}\\
583&q\equiv 13\pmod{24}\\
383&q\equiv 23\pmod{24}\end{array}\right.$$
and no nice formula seems to exist for the remaining cases (which after exclusion of powers of $2$ and $3$ form a subgroup of the multiplicative group
$(\mathbb Z/24\mathbb Z)^\ast$).
Is there an explanation for these identities?
Remark: One can of course define similarly $N(a_1,a_2)$ or $N(a_1,a_2,a_3,a_4)$.
Nice formulae for $\sum{N(\ast)\choose 2}$ exist for all primes in the first case (and are easy to prove using quadratic reciprocity). I could see nothing in the second case (computations become however quite heavy and I could not get very far).
 A: For prime $q \geq 5$ write the count as
$$
\frac1{1152} q (q-1) (q^3 - 21q^2 + 171 q - c_q)
$$
where
$$
c_q = 483 + 36 \left(\frac{-1}{q}\right) + 64 \left(\frac{-3}{q}\right)
+ \delta_q.
$$
Then for $(\frac{-2}{q}) = -1$ Ronald Bacher's calculations
indicate $\delta_q=0$.  If $(\frac{-2}{q}) = +1$ then
$q$ can be written as $m^2 + 2n^2$, uniquely up to changing
$(m,n)$ to $(\pm m, \pm n)$, and we have
$$
\delta_q = 24(m^2 - 2n^2) + 192 + 72 \left(\frac{-1}{q}\right).
$$
The explanation is as follows.  Start as did Will Sawin
by considering the variety of $(s_1,s_2,s_3,s_4,t_1,t_2,t_3,t_4)$
such that for $i=1,2,3$ the $i$-th elementary symmetric function
of the $s$'s equals the $i$-th elem.sym.fn. of the $t$'s.
We may apply any $aX+b$ transformation to all $8$ variables,
which explains the $q(q-1)$ factor.  (The factor $1152 = 2 \cdot 4!^2$
is from coordinate permutations that respect the partition of the
$8$ variables into two sets of $4$.)  In odd characteristic, there's
a unique representative with $\sum_{i=1}^4 s_i = \sum_{i=1}^4 t_i = 0$;
this takes care of the translations, and then we mod out by scalars
by going to projective space.  We end up with the complete intersection
of a quadric and a sextic in ${\bf P}^5$.  This threefold, call it
${\cal M}$, turns out to be rational.  (This has probably been known
for some time, because ${\cal M}$ classifies perfect multigrades of order $4$,
and such things have been studied since the mid-19th century, see the
Prouhet-Tarry-Escott problem; I outline a proof below.)
However, by requiring that all coordinates be distinct
we're removing some divisor ${\cal D}$ on this threefold,
so the final count decreases by the outcome of an inclusion-exclusion formula
whose terms are point counts over some subvarieties of ${\cal M}$.
Most of these sub varieties are rational curves, or points that may be defined
over ${\bf Q}(i)$ or ${\bf Q}(\sqrt{-3})$, the latter explaining
the appearance of Legendre symbols $(\frac{-1}{p})$, $(\frac{-3}{p})$
in the counting formula.  But the two-dimensional components of ${\cal D}$
are isomorphic K3 surfaces, arising as a complete intersection of
a quadric and a cubic in ${\bf P}^4$; and those components make a more complicated
contribution.  Fortunately these K3 surfaces are "singular" (i.e. their
Picard number attains the maximum of $20$ for a K3 surface in
characteristic zero) $-$ I computed that they're birational with
the universal elliptic curve over $X_1(8)$ $-$ and it is known that
the point-count of this singular K3 surface can be given by a formula
that involves $m^2-2n^2$ when $(\frac{-2}{q}) = +1$.
To show that $\cal M$ is rational, it is convenient to apply a linear
change of variables from the "$A_3$" coordinates $s_i,t_i$ to 
"$D_3$" coordinates, say $a,b,c$ and $d,e,f$, with
$$
s_i = a+b+c, \phantom+ a-b-c, \phantom+ -a+b-c, \phantom+ -a-b+c
$$
and likewise $t_i = d+e+f, \phantom. d-e-f, \phantom. -d+e-f, \phantom. -d-e+f$.
Then $\sum_{i=1}^4 s_i = \sum_{i=1}^4 t_i = 0$ holds automatically,
and the quadric and cubic become simply
$$
a^2+b^2+c^2 = d^2+e^2+f^2,
\phantom\infty
abc = def.
$$
Let $d=pa$ and $e=qb$.  Then $f=(pq)^{-1}c$, and the quadric becomes
a conic in the $(a:b:c)$ plane with coefficients depending on $p,q$:
$$
(p^2-1)a^2 + (q^2-1)b^2 + ((pq)^{-2}-1) c^2 = 0.
$$
So $\cal M$ is birational to a conic bundle over the $(p,q)$ plane,
and this conic bundle has a section $(a:b:c:d:e:f) = (1:p:pq:p:pq:1)$
which lets us birationally identify $\cal M$ with the product of the
$(p,q)$ plane with ${\bf P}^1$.  This is a rational threefold, QED.
A: Let $t_1,t_2,t_3,t_4,s_1,s_2,s_3,s_4$ be eight distinct elements of $\mathbb F_q$ such that the middle three coefficients of $(x-t_1)(x-t_2)(x-t_3)(x-t_4)$ are the same as the middle three coefficients of $(x-s_1)(x-s_2)(x-s_3)(x-s_4)$
Then $(x-t_1)(x-t_2)(x-t_3)(x-t_4)$ and $(x-s_1)(x-s_2)(x-s_3)(x-s_4)$ are two polynomials satisfying your conditions. Moreover, each pair of polynomials arises from this construction in exactly $4! \cdot 4! \cdot 2=1152$ ways, since any such pair must have distinct roots.
This is an affine variety with $8$ variables satisfying $3$ equations, hence a $5$-dimensional algebraic variety. So your polynomial is just counting points on this algebraic variety. The group of affine transformations of $\mathbb F_q$ acts faithfully on this, which explains the factor of $q(q-1)$. The remaining polynomial is the number of points on the quotient variety.
By the Lefschetz trace formula, the number of points on this variety is equal to the (alternating) trace of the action of Frobenius on its cohomology. Unfortunately, based on your computations, it seems the cohomology groups are quite large in dimension, since each power of $q$ should correspond to a single eigenvalue. So a direct computation of the cohomology would be difficult and unenlightening, unless there is some trick that I am missing.
Instead, I just want to point out that the formula you get is exactly what you would expect from a Galois representation. The various powers of $q$ come from copies of the cyclotomic character. The periodic term comes from a weight $0$ Galois representation. The dependence on $q$ mod $24$ comes from the action of the Galois group of $\mathbb Q(\mu_{24})$. The formula that is more complicated on a multiplicative subgroup and less complicated elsewhere likely indicates that part of the representation is an induced representation from that subgroup, so has trace $0$ outside it. The rest is some representation of the group that produces the desired traces on the other elements. For $q$ in the multiplicative subgroup, the point count might depend on the conjugacy class of $Frob_q$ in some nonabelian number field.
