Simple brute force search in Maple reveals that smallest positive integer $n$ is $25536$, the probable prime being
$$\lfloor 19543 \cdot \pi^{25536}\rfloor = \underbrace{3236982484 \dots 5580309289}_{12700 \text{ digits}}$$
Perhaps someone can verify primality (at least probable).
Although this doesn't answer your question, it might still be helpful:
Using Mathematica (and an exhaustive amount of computing power), I have checked every number $(9n)!+n!+1$ for $n\le 2000$ with no prime found.
Finding a prime now seems very hard. For example, if we were to model a number $n$ „being prime“ as a Bernoulli random variable with parameter $\frac{1}{\ln(n)}$ [albeit motivated by the prime number theorem, this is a very rough model, for instance one could do much better already by distinguishing even and odd numbers], and if we assume the Bernoullis to be independent for $2001\le n\le 3000$ [once again very rough], then the probability of at least one success in our model is $$1-\prod_{n=2001}^{3000}\left(1-\frac{1}{\ln\left(\left(9n\right)!+n!+1\right)}\right)\approx0.00499232.$$
If you nonetheless want to continue the search for primes, here is my Mathematica code (just replace STARTHERE and STOPHERE by the lower and upper bounds of $n$ to check):
SetSharedVariable[primes, checked]; primes = {}; checked = {};
Monitor[
ParallelDo[
If[! PrimeQ[n + 1],
If[PrimeQ[(9 n)! + n! + 1], AppendTo[primes, n]]
];
AppendTo[checked, n],
{n, STARTHERE, STOPHERE}, Method -> "FinestGrained"
],
{Sort[checked], primes}
]
EDIT: I have updated the source code because we can skip numbers $n$ for which $n+1$ is a prime number as pointed out in the comments by Sil.
EDIT 2: Here is Python source code (sadly, the prime checking function of SymPy seems to be about ten times slower than that of Mathematica)
from multiprocessing import Pool
from os import cpu_count
from sympy.ntheory import primetest
from math import factorial
import time
START = 200
STOP = 200
def check(n):
num = factorial(9*n)+factorial(n)+1
if primetest.isprime(num):
print("Found prime for", n)
return n
if __name__ == '__main__':
start_time = time.time()
with Pool(cpu_count()) as p:
primes = p.map(check, list(range(START, STOP+1)))
primes = [prime for prime in primes if prime]
print("--- {} seconds ---".format(time.time() - start_time))
Best Answer
Comment: probably this idea helps:
Wilson's theorem says:
$(p-1)+1\equiv 0 \bmod p$
Case 1-If (n+1)=q is prime we can write:
$[(n+1)-1]!+1\equiv 0\bmod (n+1)$
We rewrite $n!-(n+1)^k$ as:
$P=n!-(n+1)^K\equiv -1\bmod (n+1)+(n+1)^k\equiv -1 \bmod (n+1)$$\space\space\space\space (1)$
Case 2- If (n-1) is prime we have:
$[(n-1)-1]!+1+1\equiv 1\bmod (n-1)$
$P=1\mod (n-1)-(n+1)^k\equiv [1-(n+1)^k]\bmod (n-1)=m(n-1)-(n+1)^k+1$
considering the fact that:
$n!>(n+1)^{\lfloor \frac{n-1}2 \rfloor}$; if n is odd.
$n!>(n+1)^{\lceil \frac{n-1}2 \rceil}$; if n is even.
Substituting n by (n-1) we get:
$(n-1)!>n^{\lceil \frac{n-2}2 \rceil}$; if n is odd.
$(n-1)!>n^{\lfloor \frac{n-2}2 \rfloor}$; if n is even.
Suppose n is even and $n+1$ must be prime as OP suggest , then problem reduces to:
$$\begin {cases}P=m(n+1)-1; m,intger \\(n-1)!>n^{\lfloor \frac{n-2}2 \rfloor}\end{cases}$$
Update: The application of this theorem to relation (1) may also help:
Lioville theorem:
For prime number $p>5$ and natural number $m$ the equality $(p-1)!+1=p^m$ is not possible.
I saw this theorem and it's proof in this book:
"250 problemes de theorie elementaire des number"
By W. Sierpinsky.