My code

… was written in C++11 and can be compiled with G++, Clang++, Visual C++. You can download it, too. Or just jump to my GitHub repository.

       
#include <iostream>
       
#include <vector>
       
 
       
// ---------- standard prime sieve from my toolbox ----------

       
 
       
// odd prime numbers are marked as "true" in a bitvector

       
std::vector<bool> sieve;
       
 
       
// return true, if x is a prime number

       
bool isPrime(unsigned int x)
       
{
       
  // handle even numbers
       
  if ((x & 1) == 0)
       
    return x == 2;
       
 
       
  // lookup for odd numbers
       
  return sieve[x >> 1];
       
}
       
 
       
// find all prime numbers from 2 to size

       
void fillSieve(unsigned int size)
       
{
       
  // store only odd numbers
       
  const unsigned int half = (size >> 1) + 1;
       
 
       
  // allocate memory
       
  sieve.resize(half, true);
       
  // 1 is not a prime number
       
  sieve[0] = false;
       
 
       
  // process all relevant prime factors
       
  for (unsigned int i = 1; 2*i*i < half; i++)
       
    // do we have a prime factor ?
       
    if (sieve[i])
       
    {
       
      // mark all its multiples as false
       
      unsigned int current = 3*i+1;
       
      while (current < half)
       
      {
       
        sieve[current] = false;
       
        current += 2*i+1;
       
      }
       
    }
       
}
       
 
       
// ---------- problem specific code ----------

       
 
       
// compute n! % modulo (slow and obvious algorithm)

       
unsigned int naive(unsigned int n, unsigned int modulo)
       
{
       
  unsigned long long result = 1;
       
  while (n > 1)
       
  {
       
    result *= n;
       
    result %= modulo;
       
    n--;
       
  }
       
  return (unsigned int)result;
       
}
       
 
       
// return (base^exponent) % modulo for 32-bit values, no need for mulmod
       
// copied from my toolbox, not needed anymore because replaced by the next function

       
unsigned int powmod(unsigned int base, unsigned int exponent, unsigned int modulo)
       
{
       
  unsigned int result = 1;
       
  while (exponent > 0)
       
  {
       
    // fast exponentation:
       
    // odd exponent ? a^b = a*a^(b-1)
       
    if (exponent & 1)
       
      result = (result * (unsigned long long)base) % modulo;
       
 
       
    // even exponent ? a^b = (a*a)^(b/2)
       
    base = (base * (unsigned long long)base) % modulo;
       
    exponent >>= 1;
       
  }
       
  return result;
       
}
       
 
       
// extended Euclidean algorithm, see also en.wikipedia.org/wiki/Extended_Euclidean_algorithm
       
// simplified version of ExtendedGcd from problem 134

       
unsigned int modularInverse(unsigned int a, unsigned int modulo)
       
{
       
  // pretty much the same code can be found on https://rosettacode.org/wiki/Modular_inverse
       
 
       
  auto originalModulo = modulo;
       
  // note: s and t can be negative inside the loop
       
  int s = 0;
       
  int t = 1;
       
  while (a > 1)
       
  {
       
    auto tmp      =     modulo;
       
    auto quotient = a / modulo;
       
    modulo        = a % modulo;
       
    a = tmp;
       
 
       
    auto tmp2 = s;
       
    s = t - quotient * s;
       
    t = tmp2;
       
  }
       
 
       
  // avoid negative result
       
  return t < 0 ? t + originalModulo : t;
       
}
       
 
       
// return n! % modulo (only if modulo is prime)
       
// en.wikipedia.org/wiki/Wilson's_theorem

       
unsigned int facmod(unsigned int n, unsigned int modulo)
       
{
       
  // if n >= modulo then n! is a multiple of modulo
       
  if (n >= modulo)
       
    return 0;
       
 
       
  // for small n the naive algorithm can be faster
       
 
       
  // Wilson's Theorem is: (n-1)! == -1 (mod modulo)
       
  // avoid negative numbers:  -1 == modulo - 1 (mod modulo)
       
  unsigned long long result = modulo - 1;
       
 
       
  // (p - i)! = p! / (p-1)   / (p-2)   / (p-3)    / ... / (p - i + 1)
       
  //          = p! * (p-1)^1 * (p-2)^1 * (p-3)^-1 * ... * (p - i + 1)^-1
       
  // with modular arithmetic the inverse of a is a^-1
       
  // where a^-1 can be found with Fermat's Little Theorem
       
  for (auto i = modulo - 1; i > n; i--)
       
  {
       
    // my first approach
       
    // Fermat's Little Theorem: a^p     == a mod p
       
    // divide both sides by a:  a^(p-1) == 1 mod p
       
    // once more:               a^(p-2) == a^-1 mod p
       
    // where a^-1 is the inverse of a => equal to a^(p-2) % p
       
    // (in my code: i=a and modulo=p)
       
    // https://en.wikipedia.org/wiki/Fermat%27s_little_theorem
       
    //result *= powmod(i, modulo - 2, modulo);
       
    //result %= modulo;
       
 
       
    // my second approach (4x faster):
       
    // Extended Euclidean algorithm https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
       
    // solve gcd(a,b) = ax + by
       
    // if a and b are coprime then gcd(a,b) == 1 which means:
       
    // 1 = ax + by
       
    // when I set b=p and apply the modulo to both sides:
       
    // 1 mod p = (ax + py) mod p
       
    // 1 mod p == 1 and py is a multiple of p:
       
    // 1 = ax mod p
       
    // therefore x must be the inverse of a modulo p
       
    result *= modularInverse(i, modulo);
       
    result %= modulo;
       
  }
       
 
       
  return (unsigned int)result;
       
}
       
 
       
int main()
       
{
       
  unsigned int limit = 100000000;
       
  std::cin >> limit;
       
 
       
  // generate prime numbers
       
  fillSieve(limit);
       
 
       
  // iterate over all prime numbers
       
  unsigned long long sum = 0;
       
  for (unsigned int i = 5; i < limit; i++)
       
    if (isPrime(i))
       
    {
       
      unsigned long long minus5 = facmod(i - 5, i);
       
      // facmod is "expensive", therefore re-use the result:
       
      // facmod(i - 4, i) = facmod(i - 5, i) * (i - 4) % i
       
      unsigned long long minus4 = (minus5 * (i - 4)) % i; // same as facmod(i - 4, i)
       
      // and the same trick for the remaining values
       
      unsigned long long minus3 = (minus4 * (i - 3)) % i; // same as facmod(i - 3, i)
       
      unsigned long long minus2 = (minus3 * (i - 2)) % i; // same as facmod(i - 2, i)
       
      unsigned long long minus1 = i - 1; // again Wilson's theorem
       
 
       
      sum += (minus1 + minus2 + minus3 + minus4 + minus5) % i;
       
 
       
      // note:
       
      // of course I could customize facmod() to emit the relevant values minus1, minus2, ..., minus5
       
      // (because they are generated internally) but I wanted it to keep the code simple in case I need it for other problems
       
    }
       
 
       
  // display result
       
  std::cout << sum << std::endl;
       
  return 0;
       
}

This solution contains 27 empty lines, 62 comments and 2 preprocessor commands.