The value $\pm 1$ for the square root of Wilson's theorem, ((p-1)/2)! mod p While teaching number theory this quarter, I have come across a phenomenon which was already addressed in another MO posting, but I have new questions.  Let $p$ be a prime congruent to 3 mod 4.  Then an elementary refinement of Wilson's theorem says that $\frac{p-1}2!$ is congruent to $\pm 1$ mod $p$.  In fact, a published result of Mordell says that it is congruent to $(-1)^{(h+1)/2}$, where $h$ is the class number of the imaginary quadratic number field $\mathbb{Q}(\sqrt{-p})$.  Mordell's identity is taken to mean that (1) you cannot expect $\frac{p-1}2! \in \mathbb{Z}/p$ to be governed by any very simple property of $p$ such as a congruence, and (2) we can conjecture but not prove that the density of both values is $\frac12$.  So I wonder about the following:
(1) What is fastest known algorithm to compute $\frac{p-1}2! \in \mathbb{Z}/p$?  I read somewhere that you can compute $h$ in time $O(p^{1/5})$, but for this question (actually for the interest of the class) I only care about $h$ mod 4.
(2) Never mind the density, is it known that both values occur infinitely often?
(3) If it is known that both values occur infinitely often, then is it known, for each $\text{gcd}(n,a) = 1$, that they occur infinitely often when $p = kn+a$?  (It is of course Dirichlet's theorem that there are infinitely many such primes.)

Update:  There as an subexponential time algorithm of McCurley to compute the class number $h$.  I can believe that computing $h$ mod 4 wouldn't be any faster.
 A: As suggested in the comment by Greg Martin, the determination of the parity of $(h + 1)/2$ can be performed more efficiently than by computing the factorial expression in Lagrange’s refinement of Wilson’s theorem.
Dirichlet, “Question d’analyse indéterminée,” Journal für die Reine und Angewandte Mathematik  3 (1828): 407–408; reprinted in Werke 1:107–108, noted that its parity depends upon the number of quadratic nonresidues of $p$ lying between $1$ and $p/2$ (OEIS sequence no. A178151). This number is usually designated $m$ in the literature.
Skipping over many partial results, the value of $m$ was investigated in detail by Louis C. Karpinski in his doctoral dissertation (Mathematischen und Naturwissenschaftlichen Facultät der Kaiser Wilhelms-Universität zu Strassburg, 1903), published as “Über die Verteilung der quadratischen Reste,” Journal für die Reine und Angewandte Mathematik 127 (1904): 1–19. Karpinski proved a set of formulae involving sums over Legendre symbols, and showed that the most concise sums possible contain only $\lfloor p/6 \rfloor$ terms:
\begin{equation}
m = \frac{p-1}{4} - \frac{1}{2} \sum_{k=1}^{(p-1)/2} \left(\frac{k}{p}\right) \quad (p \equiv 3 \bmod{4}; p > 3);
\end{equation}
\begin{equation}
m = \frac{p-1}{4} - \frac{2-\left( \frac{2}{p} \right)}{3-\left( \frac{3}{p} \right)} \sum_{k=1}^{\lfloor p/3 \rfloor} \left(\frac{k}{p}\right) \quad (p \equiv 3 \bmod{4}; p > 3);
\end{equation}
\begin{equation}
m = \frac{p-1}{4} - \frac{1}{2} \sum_{k=\lfloor p/4 \rfloor +1}^{(p-1)/2} \left(\frac{k}{p}\right) \quad (p \equiv 3 \bmod{8}; p > 3);
\end{equation}
\begin{equation}
m = \frac{p-1}{4} - \frac{1}{2} \quad \sum_{k=1}^{\lfloor p/4 \rfloor} \quad \left(\frac{k}{p}\right) \quad (p \equiv 7 \bmod{8}; p > 3);
\end{equation}
\begin{equation}
m = \frac{p-1}{4} - \frac{2-\left( \frac{2}{p} \right)}{1 + \left(\frac{2}{p}\right) + \left( \frac{3}{p}\right) - \left( \frac{6}{p} \right)} \sum_{k=1}^{\lfloor p/6 \rfloor} \left(\frac{k}{p}\right) \quad (p \equiv 7, 11, 23 \bmod{24});
\end{equation}
\begin{equation}
m = \frac{-5p-1}{4} + \frac{3}{2} \sum_{k=1}^{\lfloor p/6 \rfloor} \left(\frac{k}{p}\right) \quad (p \equiv 19 \bmod{24}).
\end{equation}
A: Regarding your question (1), here are two articles I know dealing with the problem of computing factorials :
Crandall, Dilcher, Pomerance. A search for Wieferich and Wilson primes. Math. Comp.  66  (1997),  no. 217, 433--449. (MR1372002)
Costa, Gerbicz, Harvey. A search for Wilson primes. 
(Recall that Wilson primes are those primes $p$ such that $(p-1)! \equiv -1 \mod{p^2}$.)
It turns out that there exists an algorithm to compute the product of $N$ terms in arithmetic progression in no more than $O(N^{\alpha+\epsilon})$ multiplications, with $\alpha = \frac{\sqrt{5}-1}{2} = 0.618...$ (see p. 441 in the first article). They also mention an algorithm to compute $(p-1)! \mod{p^2}$ in time $O(p^{1/2+\epsilon})$.
The second article seems to use a somewhat different method and achieves an average polynomial time algorithm to compute this quantity, making use of the fact that there are redundant products when considering many values of $p$.
Both articles always work mod $p^2$. It is certainly faster to work only mod $p$, but I'm not sure whether this has an impact on the theoretical running time of the algorithms.
