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).
Best Answer
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.