Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

When working on this question regarding the divisibility of the sum of factorials, I decided to write some code to test "small values" of the problem using the following code.

f[p_] := Total[Mod[#!, p] & /@ Range[p - 1]];
Table[Mod[f@Prime@i, Prime@i], {i, 1, 500}]

Basically, what the code does is sum up all the factorials $$1!+2!+3!+\dots+(p-1)!$$

and find the remainder modulo $p$, for prime $p$.

Unfortunately, my code as written takes a very long time to run. Checking the first 500 primes takes 88.280966 seconds on my computer, but checking the first 2000 primes took me about 4 hours.

Is there any way to improve the code, or is it already the best we can do?

As for optimizations not involving the code, I used Wilson's Theorem, which states that for all primes $p$,

$$(p-1)!\equiv-1 \bmod p$$

Using the above theorem, we can modify the code as follows.

h[p_] := Total@Flatten[{Mod[#!, p], PowerMod[(# - 1)!*(-1)^(#), -1, p]} & /@ Range[(p - 1)/2]];
Table[Mod[h@Prime@i, Prime@i], {i, 1, 500}]

This is considerably faster than the previous code, since checking the first 500 primes takes only 25.896166 seconds. However, checking the first 2000 primes still takes an inordinately long time.

share|improve this question
up vote 13 down vote accepted

This is bit faster:

toPrime = 500;
sums = Accumulate@FoldList[Times, 1, Range[2, Prime@toPrime - 1]];
primes = Prime[Range[toPrime]]; 
Mod[sums[[primes - 1]], primes]

Precompute factorial sums and primes. Mod is fast on lists.

share|improve this answer
you are way too humble! calculating your sum for even the first 2000 primes takes less than a second. However, is there a way to get around storing large numbers in sums in memory? It keeps crashing my computer when I try toPrime=5000. – Vincent Tjeng Mar 26 '13 at 6:03
+1 (that's freaking fast!) Can you explain why Accumulate@FoldList[#1 #2 &, 1, Range[n] is so much quicker to Accumulate@Array[#! &, n] + 1? I really don't get it. – gpap Mar 26 '13 at 11:23
@gpap Calculating factorial so many times costs a lot. Since we know we want all the factorials up to Prime[toPrime]-1, we ultimately gain a lot keeping the intermediate results with FoldList. – Michael E2 Mar 26 '13 at 12:04
@MichaelE2 Yes, worked it out myself in the meantime - you multiply the previous result and don't calculate a factorial at every step. Thanks – gpap Mar 26 '13 at 12:11
Is there any reason why you used #1 #2 & instead of Times? – J. M. Jun 15 '15 at 12:53

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.