Given a prime $p$, an integer $0\le n\le p-1$, and a divisor $q\mid(1+p+p^2)$, we compute $n!\bmod p$ in expected bit complexity $\widetilde{O}\left(q^c+\frac{\sqrt{p}}{q^{1/4}}\right)$ for some absolute constant $c\ge1$. More generally, the construction applies when $q\mid\Phi_r(p)$, where $\Phi_r$ is the $r$-th cyclotomic polynomial and $r$ is any fixed odd prime power. Combining these constructions, for every fixed $\epsilon\in(0,1/2)$, we obtain an expected bit complexity of $p^{1/2-\delta_\epsilon+o(1)}$, with positive $\delta_\epsilon$, for at least a $1-\epsilon$ fraction of primes up to $X$, for all sufficiently large $X$. Under the same divisor conditions, the method extends to moduli $p^k$ with an additional factor polynomial in $k$.

The method recovers ratios of factorials modulo $p$ from Jacobi sums, character sums over finite fields. We use Lenstra and Silverberg's algorithm to reconstruct these sums up to a root of unity from their ideal factorizations and products with their complex conjugates. We determine this root using van Wamelen's criterion. Lattice rounding selects nearly equal parts while keeping the remaining factorials small, so each recursive step uses only one smaller factorial and short interval products. The only randomized steps are Las Vegas constructions of the finite-field representations and multiplicative characters.

Related constructions use suitable divisors of $p-1$ for moduli $p$ and $p^2$, and of $p+1$ for modulus $p$. A separate deterministic algorithm uses modular halving and simultaneous evaluation to improve the Bostan–Gaudry–Schost bound modulo $p$ and $p^2$ by a factor of $\sqrt{\log p/\log\log p}$, without any divisor assumption.