Is there a prime number $p$ dividing $1+2!^2+3!^2+\cdots (p-1)!^2$?

Question

Is there a prime number $p$ with $p \mid \sum_{j=1}^{p-1} j!^2$?

I checked the primes upto $600\ 000$ without finding a solution.

Heuristic : If we can assume that the probability that $p$ is a solution is $\frac{1}{p}$ , then there should be infinite many solutions and therefore the desired prime should exist. But I guess this is not the case and maybe the special expression can easily be determined to be or not to be divisible by $p$.

Motivation : Such a prime number would prove that $1+2!^2+3!^2+\cdots n!^2$ can be prime only for finite many positive integers $n$ and also give an upper bound for the possible $n$. If the answer to the question is no , there is a chance that there are infinite many such primes.

Answer 1

Yes. By running a script which computes this sum mod $p$ through the first $100\,000$ primes, I found one solution: $p = 1\,248\,829$ ($96\,379$th prime).

UPDATE. This is also the only solution for the first $1\,000\,000$ primes.

UPDATE 2. If the probability that $p$ is the solution is indeed $1/p$, it follows from the Mertens' third theorem, that the probability to find the solution between $10^{m-1}$th prime and $10^m$th prime is asymptotically $1/m$. This could explain our serendipitous finding, but also discourages to look any further.

Answer 2

Let's implement $S$ in Python to make large searches possible (for whoever is interested).

Let $S(m) = \sum_{j=1}^{m-1} (j!)^2$. We want to focus our study on prime $m$, but let's have $S$ available for any positive integer $m \geq 2$.

(There is an extensive comment in the Markdown here. It details an attempt to move directly from $S(m-1) = q_{m-1}(m-1) + r_{m-1}$ to $S(m) = q_m m + r_m$. The result was no faster, so it is not presented. Feel free to edit and read if you are somehow interested.)

The function sPlain() is a direct translation of the definition of $S$. The function sReuse is an attempt to implement Efim Mazhnik's description of his code. Sconsecutive() evaluates $S(m) \pmod{m}$ for $m =2, 3, 4, \dots$. SconsecutivePrimes() evaluates $S(p) \pmod{p}$ for prime integers, $p = 2, 3, 5, \dots$.

(I find that SconsecutivePrimes() runs faster as (equivalent) Mathematica code. Most likely, Mathematica tries harder to cache sieving partial results between NextPrime[] calls than does sympy for nextprime(), but I haven't really looked into it.)

Sconsecutive() has been run up to 3.1 million finding only $m \mid S(m)$ for $m = 1\,248\,829$. SconsecutivePrimes() has been run up to 10 million and found no new successes.

from math import factorial

def Splain(m):
    retVal = 0
    for j in range(1, m):    # iterates over [1,m).
        retVal += factorial(j)**2
    return retVal

def Sreuse(m):
    retVal = 0
    jFactorialSquared = 1
    for j in range(1, m):
        jFactorialSquared *= j**2
        retVal += jFactorialSquared
    return retVal

def Sconsecutive():
    # Has been run to m = 3.1 million (and a little further); only found 1248829|S(1248829).
    m = 2  #  least m for which S(m) is not empty
    j = 1  #  new index in S(2) compared to S(1)
    jFactorialSquared = factorial(j)**2
    s = jFactorialSquared  #  S(m) = S(2) = 1!^2 = 1
    print("Runs forever.  Press <ctrl>-c (or whatever interrupts execution in your environment) to stop.")
    while True:
        if s % m == 0:
            print(" m | S(m) for m = ",  m)
        m+=1
        j += 1
        jFactorialSquared *= j**2
        s += jFactorialSquared
        if m % 100000 == 0:  #  Depending on your hardware, you might like feedback more frequently that every hundred thousand.
            print("About to test m = ",  m)

from sympy import nextprime
def SconsecutivePrimes():
    # Has been run to p = 10 million (and a little further); only found 1248829|S(1248829).
    p = 2  #  S(2) is first nonempty sum
    s = 1
    jTerm = 1    # = factorial(1)**2
    count = 1
    print("Runs forever.  Press <ctrl>-c (or whatever interrupts execution in your environment) to stop.")
    while True:
        if s % p == 0:
            print(" p | S(p) for p = ",  p)
        pNext = nextprime(p)
        # Note: sympy's nextprime() promises to deliver primes for p less 
        # than about 2^32.  For larger p, nextprime() seems to produce strong 
        # pseudoprimes to several bases.  It's fine that we use such ps, 
        # but we should rigorously check for primeness before reporting 
        # divisibility.  It's left as an exercise for the reader to use, 
        # for example, 
        #     https://github.com/wacchoz/APR_CL/blob/master/APR_CL.py
        # to modify the "if s % p == 0:" block above to check rigorously 
        # for primality before reporting a "success".
        s += sum ([ jTerm := jTerm * (j**2) for j in range(p,  pNext) ])
        p = pNext
        count +=1
    
        if count % 10000 == 0:
            print("# About to test p = ", p)

from decimal import *
if __name__ == "__main__":
    print("The Splain() call may take longer than you expect.  The Sreuse() call will return much faster.")
    print("Splain(10,000) = ",  Splain(10_000))
    print("Sreuse(10,000) = ",  Sreuse(10_000))
    SconsecutivePrimes()