11
$\begingroup$

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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
15
$\begingroup$

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.

Improve this answer
$\endgroup$
9
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Mar 26, 2013 at 6:03
  • $\begingroup$ +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. $\endgroup$ Commented Mar 26, 2013 at 11:23
  • $\begingroup$ @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. $\endgroup$ Commented Mar 26, 2013 at 12:04
  • $\begingroup$ @MichaelE2 Yes, worked it out myself in the meantime - you multiply the previous result and don't calculate a factorial at every step. Thanks $\endgroup$ Commented Mar 26, 2013 at 12:11
  • 2
    $\begingroup$ Is there any reason why you used #1 #2 & instead of Times? $\endgroup$ Commented Jun 15, 2015 at 12:53
6
$\begingroup$

Let $x \equiv r_1 \bmod p$ and $y \equiv r_2 \bmod p$. Then, $x y \equiv r_1 r_2 \bmod p$. So, we can compute the sum of the factorials mod p using:

f[p_] := Mod[ Total @ FoldList[ Mod[Times[##], p]&, Range[p-1]], p]

Let's compare this to the naive implementation:

t[p_] := Mod[Sum[k!, {k, p-1}], p]

For example:

f[Prime[500]] //AbsoluteTiming t[Prime[500]] //AbsoluteTiming

{0.00058, 1813}

{0.085628, 1813}

The nice thing about using Mod[Times[##], p] as the FoldList function is that the output should be a machine number unless you are working with very large primes. That means that f can be compiled:

fc = Compile[{{p, _Integer}}, Mod[ Total @ FoldList[ Mod[Times[##], p]&, Range[p-1]], p], RuntimeAttributes->{Listable} ];

Let's compare for a large prime:

f[Prime[10^6]] //AbsoluteTiming fc[Prime[10^6]] //AbsoluteTiming

{1.83041, 9308538}

{1.05449, 9308538}

Faster, but the real difference is that fc is Listable. Hence, comparing timings on a list shows a much larger difference. For example:

r1 = f /@ Prime[Range[5000, 5500]]; //AbsoluteTiming r2 = fc[Prime[Range[5000, 5500]]]; //AbsoluteTiming r1 === r2

{2.06691, Null}

{0.395402, Null}

True

Finally, since the list of all factorials is not stored anywhere, the memory used is much more manageable. Here is the memory and timing for the first 5000 primes:

r1 = fc[Prime[Range[5000]]]; //MaxMemoryUsed //AbsoluteTiming

{3.03463, 462848}

This is quite a bit smaller than @MichaelE2's answer:

MaxMemoryUsed[ toPrime = 5000; sums = Accumulate@FoldList[Times, 1, Range[2, Prime@toPrime - 1]]; primes = Prime[Range[toPrime]]; r2 = Mod[sums[[primes - 1]], primes] ] //AbsoluteTiming r1 === r2

{2.40441, 4064607008}

True

The compiled answer uses .462KB while the approach where all the factorials are precomputed takes 4GB, and the timing is not too different.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.