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
Not a full answer, just a thought on the prime randomness:
The fibonacci numbers are about
$$F_n\sim \frac{\phi^{n}}{\sqrt{5}}$$
where $\phi=\frac{1+\sqrt{5}}{2}$. Then $f(n)$ grows as
$$f(n)=\prod_{i=1}^n F_i\sim \prod_{i=1}^n \frac{\phi^{n}}{\sqrt{5}}=\sqrt{\frac{\phi^{n(n+1)}}{5^{n}}}$$
Then the probability that there is a prime $1$ above or below this number is about
$$P=2\cdot \frac{2}{\log\left(\frac{\phi^{n(n+1)}}{5^{n}}\right)}=\frac{4}{n(n+1)\log(\phi)-n\log(5)}$$
which is $O(n^{-2})$. Since
$$\sum_{n=1}^\infty \frac{1}{n^2}<\infty$$
we would expect that the number of such primes to be finite.