I am a newcomer here. If p >3 is congruent to 3 mod 4, there is an answer which involves only $p\pmod 8$ and $h\pmod 4$, where $h$ is the class number of $Q(\sqrt{-p})$ .
Namely one has $(\frac{p-1}{2})!\equiv 1 \pmod p$ if an only if either (i) $p\equiv 3 \pmod 8$ and $h\equiv 1 \pmod 4$ or (ii) $p\equiv 7\pmod 8$ and $h\equiv 3\pmod 4$.
The proof may not be original: since $p\equiv 3 \pmod 4$, one has to determine the Legendre symbol
$${{(\frac{p-1}{2})!}\overwithdelims (){p}}
=\prod_{x=1}^{(p-1)/2}{x\overwithdelims (){p}}=\prod_{x=1}^{(p-1)/2}(({x\overwithdelims (){p}}-1)+1).$$ It is enough to know this modulo 4 since it is 1 or -1. By developping, one gets $(p+1)/2+S \pmod 4$, where
$$S=\sum_{x=1}^{(p-1)/2}\Bigl({x\over p}\Bigr).$$ By the class number formula, one has $(2-(2/p))h=S$ (I just looked up Borevich-Shafarevich, Number Theory), hence the result, since $\Bigl({2\over p}\Bigr)$ depends only on $p \pmod 8$.
Edit: For the correct answer see KConrad's post or Mordell's article.
I just wanted to remark that if $p$ is a prime such that
$\ell$ splits in $F = \mathbb{Q}(\sqrt{-p})$ for all $\ell \le N$,
then one may prove the existence of $N$ consecutive integers which
are norms of integers in $\mathcal{O}_F$, providing one is willing
to assume a standard hypothesis about prime numbers, namely,
Schinzel's Hypothesis H.
First, note the following:
Lemma 1: If $C$ is an abelian group of odd order, then there exists a finite
(ordered) set $S = \{c_i\}$ of elements of $C$ such that every element in $C$ can be
written in the form
$\displaystyle{\sum \epsilon_i \cdot c_i}$
where $\epsilon_i = \pm 1$.
Proof: If $C = A \oplus B$, take $S_C = S_A \cup S_B$. If $C
= \mathbb{Z}/n \mathbb{Z}$ then take $S = \{1,1,1,\ldots,1\}$ with
$|S| = 2n$.
Let $C$ be the class group of $F$. It has odd order, because
$2$ splits in $F$ and thus $\Delta_F = -p$.
Let $S$ be a set as in the lemma. Let $A$ denote an ordered
set of distinct primes $\{p_i\}$ which split in $\mathcal{O}_F$ such
that one can write $p_i = \mathfrak{p}_i \mathfrak{p}'_i$ with
$[\mathfrak{p}_i] = c_i \in C$, where $c_i$ denotes a set of elements
whose existence was shown in Lemma 1.
Lemma 2: If $n$ is the norm of some ideal $\mathfrak{n} \in \mathcal{O}_F$,
and $n$ is not divisible by any prime $p_i$ in $A$, then
$$n \cdot \prod_{A} p_i$$
is the norm of an algebraic integer in $\mathcal{O}_F$.
Proof: We may choose
$\epsilon_i = \pm 1$ such that
$\displaystyle{[\mathfrak{n}] + \sum \epsilon_i \cdot c_i = 0 \in C}$.
By assumption, $[\mathfrak{p}_i] = c_i \in C$ and thus
$[\mathfrak{p}'_i] = -c_i \in C$. Hence the ideal
$$\mathfrak{n}
\prod_{\epsilon_i = 1}
{ \mathfrak{p}}
\prod_{\epsilon_i = -1}
\mathfrak{p}'$$
is principal, and has the desired norm.
By the Chebotarev density theorem (applied to the Hilbert class field
of $F$), there exists a set $A$ of primes as above which avoids any fixed finite set of
primes. In particular, we may
find $N$ such sets which are pairwise distinct and which contain
no primes $\le N$. Denote these sets by $A_1, \ldots, A_N$.
By the Chinese remainder theorem, the set of integers $m$ such that
$$m \equiv 0 \mod p \cdot (N!)^2$$
$$m + j \equiv 0 \mod \prod_{p_i \in A_j} p_i, \qquad 1 \le j \le N$$
is of the form $m = d M + k$ where $0 \le k < M$, $d$ is arbitrary, and $M$
is the product of the moduli.
Lemma 3: Assuming Schinzel's Hypothesis H, there exists infinitely many integers
$d$ such that
$$ P_{dj}:= \frac{dM + k + j}{j \cdot \prod_{p_i \in A_j} p_i}$$
are simultaneouly prime for all $j = 1,\ldots,N$.
Proof: By construction, all these numbers are coprime to $M$ (easy check).
Hence, as $d$ varies, the greatest common divisor of the product of these
numbers is $1$, so Schinzel's Hypothesis H applies.
Let $\chi$ denote the quadratic character of $F$. Note that
$dM + k + j = j \mod p$, and so
$\chi(dM + k + j) = \chi(j) = 1$ (as all primes less than $N$ split in $F$).
Moreover, $\chi(p_i) = 1$ for all
primes $p_i$ in $A_j$
by construction. Hence
$\chi(P_{dj}) = 1$. In particular, if $P_{dj}$ is prime, then
$P_{dj}$ and $j \cdot P_{dj}$ are norms of (not necessarily principal) ideals in
the ring of integers of $F$. By Lemma 2, this implies that
$$dM + k + j = j \cdot P_{dj} \prod_{p_i \in A_j} p_i$$
is the norm of some element of $\mathcal{O}_F$ for all $j = 1,\ldots, N$.
One reason to think that current sieving technology will not be sufficient
to answer this problem is the
following: when Sieving produces a non-trivial lower bound, it usually
produces a pretty good lower bound. However, there are no good (lower) bounds known for the following problem: count the number of integers $n$
such that $n$, $n+1$, and $n+2$ are all sums of two squares. Even for the
problem of estimating the number of $n$ such that $n$ and $n+1$ are both sums of squares is tricky - Hooley implies that the natural sieve does not give lower bounds
(for reasons analogous to the parity problem). Instead,
he relates the problem to sums of the form
$\displaystyle{\sum_{n < x} a_n a_{n+1}}$ where
$\sum a_n q^n = \theta^2$ is a modular form. In particular, he
implicitly uses
automorphic methods which won't work
with three or more terms.
Best Answer
Here's a massive generalization, with references.
Consider a sequence $\{a_n \}_{n\ge 1}$ of integers. We'll say it satisfies the necklace congruences if $a_{n} \equiv a_{n/p} \mod n$ whenever $p \mid n$ ($p$ prime).
Here are some equivalent formulations:
So any such sequence gives rise to an integer sequence $\lambda_{\tilde{a}}(n) := \frac{1}{n}\sum_{d \mid n} a_d \mu(n/d)$ which satisfies $a_n = \sum_{d \mid n} d \lambda_{\tilde{a}}(n)$ by Mobius inversion. One can then define $\sigma_{\tilde{a}}(n):=\sum_{d \mid n, 1<d<n} d\lambda_{\tilde{a}}(d)$ and obtain the following generalization of Fermat's little theorem:
$a_n \equiv a_1 + \sigma_{\tilde{a}}(n) \mod n$
Now, let's go back to necklace congruences. A less obvious equivalence is the following:
Theorem: $\{a_n\}_{n \ge 1}$ satisfies the necklace congruences iff $\zeta_{a}(x):=\exp(\sum_{n\ge1} \frac{a_n}{n}x^n)$ has integer coefficients.
Proof: Write, $\zeta_a$ formally as $\prod_{n\ge 1} (1-x^n)^{-\frac{b_n}{n}}$. The $b_n$'s are uniquely determined, and in fact (by taking logarithmic derivative) it can be seen that they are $b_n = \sum_{d \mid n} a_d \mu(n/d)$. Now one direction is immediate and the other requires an inductive argument. $\blacksquare$
(Reference: Exercise 5.2 (and its solution) in Richard Stanley's book "Enumerative Combinatorics, vol. 2")
As Sergei remarked, for any integer square matrix $A$, the sequence $a_n := \text{Tr}(A^n)$ satisfies the necklace congruences. This is immediate from the last theorem, as the corresponding "zeta" function $\zeta_a$ is just:
$\exp( \sum_{n \ge 1} \frac{\text{Tr}((Ax)^n)}{n}) = \exp (\text{Tr} (- \ln (I -Ax)))=\det(I-Ax)^{-1} \in \mathbb{Z}[x]$
In particular, by taking a 1 on 1 matrix, we recover your observation.
Another generalization is $a_n = [x^n]f^n(x)$, where $f \in \mathbb{Z}[[x]]$. By taking the linear polynomial $f(x)=1+ax$ we recover your example, by taking $f(x)=(1+x)^m$ we get $a_n = \binom{nm}{n}$.
The really interesting feature is that this notion has a local version. Specifically, theorems of Dieudonne-Dwork and Hazewinkel give the following result:
Theorem: Given $\{ a_n \}_{n\ge 1} \subseteq \mathbb{Q}_{p}$ ($p$-adic numbers), the following are equivalent:
Applying this simultaneously for all $p$, we recover the previous theorem. A reference for this theorem is "A Course in p-adic Analysis" by Alain M. Robert (the theorems are stated there in a much more general context).