Second part of the factorial sum divisibility question Which primes $p$ divide the sum of factorials $1! + 2! + 3! + 4! + 5! + \cdots + (p-1)!$? This is related to my previous question.
 A: Gottfried's answer and Steven's replies bring up the point that there's a choice to be made in a direct numerical search: either you can compute $(p-1)!$ $(\textrm{mod }p)$ separately for each $p$, or you can store $(p-1)!$ and then compute $p! = p(p-1)!$ to test the next $p$.  This same tradeoff comes up in other similar problems, such as the search for Brown Numbers.
In order to check every $p<N$, the first approach requires $O(N^2)$ multiplies of numbers of order $p$; the second approach requires $N$ multiplies, but each multiply is very expensive, being between a number of order $p$ and one of order $p!$.  Another consideration is that the second solution, storing $(p-1)!$ un-reduced, requires memory linear in $p$.  I don't really want to see what happens when $(p-1)!$ no longer fits in cache.  Furthermore, it's not obvious how to parallelize the second algorithm.
My impression is that if you want to check all $p<N$ and not just the primes $p$ it's fastest to store $(p-1)!$ for small $N \lesssim 10^6$, but for larger $N$ it's faster to compute $(p-1)!$ $(\textrm{mod }p)$ separately for each $p$.  If you want to only check primes $p$, then the "first approach" gets a $O(\log(n))$ speedup and the "second approach" hardly benefits at all, meaning that the first approach is faster: see the "Update:" section below the code for an explanation.
Here is a C code that uses the GMP to implement this "second approach"; I implemented the "first approach" in my first answer above.
#include <gmp.h>
#include <time.h>

#define NMAX 100000

int main(){

    // Print begin time & date
    time_t beginTime;
    struct tm * timeinfo;
    time ( &beginTime );
    timeinfo = localtime ( &beginTime );
    gmp_printf("# Begin time/date is: %s", asctime (timeinfo) ); 

    // Stores n! and sum(n!)
    mpz_t cur_factorial;
    mpz_t sum_factorial;
    mpz_t sum_factorial_modn;
    mpz_init_set_ui(cur_factorial, 1);
    mpz_init_set_ui(sum_factorial, 1);
    mpz_init(sum_factorial_modn);

    for(unsigned long int n=2; n<=NMAX; n++){

      // Check to see if n divides sum(m=1..n-1) m!
      if( mpz_divisible_ui_p(sum_factorial, n) )
    gmp_printf("%u\n", n);

      // Update n!
      mpz_mul_ui(cur_factorial, cur_factorial, n);

      // Update sum(n!)
      mpz_add(sum_factorial, sum_factorial, cur_factorial);

    }

    // Print end time & date; total run time
    time_t endTime;
    time ( &endTime );
    timeinfo = localtime ( &endTime );
    gmp_printf("# End time/date is:   %s", asctime (timeinfo) ); 
    gmp_printf("# Total runtime of %.0f seconds.\n", difftime(endTime, beginTime) ); 


    // Release and return
    mpz_clear(cur_factorial);
    mpz_clear(sum_factorial);
}


Update:
This code implementing the "second approach" is strictly slower than the one in my previous answer. The above code takes ~9s to for an upper limit to $10^5$, which only takes ~5s using the code in my previous answer.  A search up to $10^6$ takes ~19 min with this code, compared to ~8 min with the other one.  Part of the reason for this is that, with the other code, you get a huge speedup if you only test primes $p$.  This code is set up to check all $n$, even the pointless even ones.  This might not seem like a fair comparison, but modifying the above code to only check odd $n$ doesn't give any speed improvement; the point is that you anyway have to calculate $(n-1)!$ for every $n$ with this method, so you've already paid the cost; the divisibility check doesn't take any time compared with the multiplies.
A: In addition to @Douglas post, here is a very short (and I think: efficient) implementation in Pari/GP - however, I didn't run it to high N myself:
[Update] I added an improved computation method at the end[/update]
MaxN=101   \\ set some upper limit for the computation
S=[1,1]    \\ contains current factorial and current sum
for(k=3,MaxN, S *= [k-1,k-1;0,1]; if(isprime(k),print([k,S[1] % k,S[2]%k])))

Result :     
p | (p-1)! %p | S %p
--+-----------+-----
[3, 2, 0]
[5, 4, 3]
[7, 6, 5]
[11, 10, 0]
[13, 12, 9]
[17, 16, 12]
[19, 18, 8]
[23, 22, 20]
[29, 28, 16]
[31, 30, 1]
[37, 36, 4]
[41, 40, 3]
...    

Up to the 1000'th prime MaxN=7919      
MaxN = prime(1000)
...
...
[7919, 7918, 2882]

this needed 281 msec.
A Version using the modulo-computation is the following
gettime()
{for(j=1,1000,p=prime(j); 
   S=[1,1];
   for(k=3,p,
        S*=[k-1,k-1;0,1];
         S=S % p
       );
   print([p,S[1] % p,S[2]%p])
   )}
gettime()

%1124 = 9219   \\ this is msecs


[update] The following is only relevant, if we want a routine, which produces the list over all primes up to some limit.
It has then some advantage to do a compromise between the first method where we need all factorials and sums only once but which needs many bits for the representation of the numbers and the complete modulo-computation, which operates on small numbers but must do the whole run of computing factorials (residuals) for each prime such that we have a double loop.
If we compute and store the required factorials and partial sums in chunks of, say, $t=10$ terms, such that we have a list of precomputed values 
$$ \small \begin{array} {r|l|l|} 
       & f   & S \\
\hline \\
  c_1=& 10! & \sum_{k=1}^{10} k! \\
  c_2=& {20!\over 10!} & {\sum_{k=11}^{20} k! \over 10!} \\
  c_3=& {30!\over 20!} & {\sum_{k=21}^{30} k! \over 20!} \\
  \ldots
 \end{array} $$
then we can compose our results for each prime by modular multiplication of that constants and a short remainder. Finetuning the parameter $t$ for higher allows then to find a compromise between storage/size of the stored numbers and number of steps in the inner loop, which contains the modular computation per prime.
Using some empirical data I have done an estimation with a quadratic curve in Excel, which says, for the list of primes up to 1000000 I need 90 min with blocksize $t=20$ , 30 min with $t=100$ and 15 min with $t=500$ .
(If it is wished I can append the Pari/GP code here)
A: This is what I thought about, it should be some comment, but this exceeds the character limit for comment.
If you take Wilson Theorem that 
$(p-1)!\equiv -1 \mod{p}$, then
$(p-2)!\equiv 1 \mod{p}$.
If we take inverse of $a$ in $F_p$ as $\dfrac{1}{a}$, we know that $\dfrac{1}{a} + \dfrac{1}{b} = \dfrac{a+b}{ab}$ in $F_p$, it can taken as usual addition.
Then 
$(p-3)! \equiv -\dfrac{1}{2}\mod{p}$
$(p-4)! \equiv -\dfrac{1}{2}\cdot\left(-\dfrac{1}{3}\right)\mod{p}$ 
$\cdots$
$1! \equiv (-1)^{p-3}\dfrac{1}{(p-2)!} \mod{p}$
Thus 
$1!+2!+\cdots(p-3)! \equiv -\sum_{k=0}^{p-2}(-1)^k\dfrac{1}{k!}\mod{p}$
And we state without proof.
Lemma:
Consider the fraction $\dfrac{p}{q}$ with $q = m!$, then the closest number to $\exp(-1)$ with form $\dfrac{p}{q}$ should satisfy
$$\dfrac{p}{q} = \sum_{j=0}^m (-1)^j\dfrac{1}{j!}$$
Then 
$-\sum_{k=0}^{p-2}(-1)^k\dfrac{1}{k!} \equiv - \mathrm{floor}\left[\dfrac{(p-2)!}{\exp(1)}\right]\mod{p}$
However, this cannot be checked for large $p$, since the constant $\exp(1)$ is not accurate enough. I simply checked for smaller primes.
However, I felt that there should be other solutions than 3 and 11, otherwise, this can be a very interesting property for the constant.
A: Note that heuristically, if these numbers were equidistributed mod $p$ we would expect each one to be zero (i.e., expect $p$ to divide the result) with probability $\approx 1/p$; presuming that all the values are independent (which seems a reasonable assumption), we shoud expect to have $\sum_{p\lt n}\frac{1}{p}\approx M+\ln \ln n$ of them less than $n$, where $M$ is the Meissel-Mertens Constant — this sum is approximately 2.887 for $n=10^6$, so while we would heuristically 'expect' another one in that range it's not surprising that there aren't any, and despite the apparent lack of any more solutions, because of the divergence of the series the 'expected' number of solutions is still infinite!
Added: I put together my own software implementation in C++, confirming the results (that only $p=3$ and $p=11$ work) for $p\lt 5\times 10^5$ but also looking at the distribution of results $\mod p$ to see if any patterns arose.  I computed the value of $\dfrac{\sum_{1}^{p-1}i!\bmod p}{p}$ (the division by $p$ serving to 'normalize' values into the range $[0, 1)$ ) for all $p$ in that range and then grouped them in bins; these are the plots with 100 bins and 256 bins.  (These bin counts should be small enough relative to the prime values in consideration that there shouldn't be any substantial aliasing effects in binning the data.)  The results are a bit scattershot but there doesn't appear to be any skew towards particular values (e.g. $\frac{p-1}2$) which would show up as a spike in the bins.  If I have more time I may take a look at the square of the sum and particularly even the inverse of the sum $\pmod p$, but since the latter would actually require me writing a quick GCD algorithm to find the inverse it'll have to wait for later.
(Note: the sharp drop at the last element is an artifact of my adding a 0 value to the end of the bin data so my spreadsheet didn't automatically normalize the range; it's not an actual value.)


A: I have a few small analytic comments and a relatively aggressive numerical search.
Firstly, I checked all odd $p$ (not necessarily prime) up to $10^6$ and found only the same five as Yoni above, namely 3, 9, 11, 33, and 99.  This took ~1 hour using 1 core of my MacBook Pro, with an early unoptimized version of my code.
Secondly, of course we have:
\begin{equation}
p\nmid\sum_{n=1}^{p-1}n!\quad \forall p\in\mathbb{N}, p \textrm{ even}
\end{equation}
because all of the terms in the sum are even except the first one, so the sum is always odd.
Thirdly, a numerical search is sped-up slightly using the identity:
\begin{equation}
\sum_{n=1}^{p-1} n! = 1+2[1+3[1+4[...[1+(p-1)]]]]
\end{equation}
I've written software to check primes $p$ for $p|\sum_{n=1}^{p-1} n!$ using the above identity, and set it up to use multiple threads in an embarrassingly-parallel fashion.  Using this software I've shown:
\begin{equation}
\textrm{3 and 11 are the only primes }p<10^7\textrm{ that divide }\sum_{n=1}^{p-1} n!
\end{equation}
If you want to compile it yourself then you'll have to link with the primesieve library.
#include <stdio.h>
#include <stdlib.h>
#include <primesieve/soe/PrimeSieve.h>


// To compile: g++ -O3 factorialSum.cpp -lprimesieve -o factorialSum
//         or: g++ -O3 -static factorialSum.cpp -lprimesieve -o factorialSum


// We perform arithmetic at u1 precision; p is restricted to u2 precision
typedef uint_fast64_t u1;
typedef uint_fast32_t u2;


// This evaluates $\sum_{n=1}^{p-1} n!$ based on the definition.
// We don't actually use this in our program below.
u2 factorialSum_modp_method1(u2 p){
  u1 factorial_modp=1;
  u1 factorialSum_modp=0;
  for(u1 n=1; n<=p-1; n++){
    factorial_modp*=n;
    factorial_modp%=p;
    factorialSum_modp+=factorial_modp;
    factorialSum_modp%=p;
  }
  return factorialSum_modp;
}


// This evaluates $\sum_{n=1}^{p-1} n! = 1+2[1+3[1+4[...[1+(p-1)]]]]$ via the RHS of the expression.
// This is somewhat faster than 'method1' above.
u2 factorialSum_modp_method2(u2 p){
  u1 factorialSum_modp=1;
  for(u1 n=p-1; n>1; n--)
    factorialSum_modp=(factorialSum_modp*n+1)%p;
  return factorialSum_modp;
}


// This prints out p iff p divides $\sum_{n=1}^{p-1} n!$.
// We pass this as an arument to PrimeSieve.generatePrimes().
void check_factorialSum_modp(u2 p){
  if(factorialSum_modp_method2(p)==0)
    printf("%u\n", (unsigned)p);
}


void error_incorrectUsage(char *args[]){
  printf("Usage: %s thread numThreads blockSize pMin pMax\n", args[0]);
  printf("       thread is a positive integer 0, ..., numThreads-1\n");
  printf("       numThreads is a natural number, 1, 2, ...\n");
  printf("       blockSize is the interleave size between threads\n");
  printf("       pMin and pMax are integers, the upper and lower limits on candidate primes p\n\n");
  printf("e.g.:  %s 0 1 10000 0 100000\n", args[0]);
  printf("       will check [0, 100000] for primes p that divide $\\sum_{n=1}^{p-1} n!$ \n");
  printf("       (expression written LaTeX).\n\n");
  printf("In order to use multiple threads, run multiple instances with different values of 'thread'.\n");
  exit(1);
}


int main(int nargs, char *args[]){

  // Parse command line arguments and check for proper usage
  if( !(nargs==6) ) error_incorrectUsage(args);
  char *error1, *error2, *error3, *error4, *error5;
  u2 thread     = strtoul(args[1], &error1, 10);
  u2 numThreads = strtoul(args[2], &error2, 10);
  u2 blockSize  = strtoul(args[3], &error3, 10);
  u2 pMin       = strtoul(args[4], &error4, 10);
  u2 pMax       = strtoul(args[5], &error5, 10);
  if(*error1 || *error2 || *error3 || *error4 || *error5) error_incorrectUsage(args);

  // Sieve intervals [pLower, pUpper] that belong to this thread
  PrimeSieve ps;
  for(u2 pLower=pMin+thread*blockSize; pLower<=pMax-blockSize; pLower+=numThreads*blockSize){
    u2 pUpper=pLower+blockSize;
    ps.generatePrimes(pLower, pUpper, check_factorialSum_modp);
    printf("interval [%u,%u] completed on thread %u\n", pLower, pUpper, thread);
  }

}

