The answer to the question is no, here is a slightly smaller solution:
$$
494\;804\;473\ .
$$
sage: factor( factorial(40) + 494804473 )
8912658466556113 * 9232286052422441 * 9915818101022081
It was obtained also in sage, the program was running one day, and this is in this intermediate situation a found solution. (It is still running. It was not written to touch the numbers one by one starting with $40!$ in order, so this may also not be the minimal number doing the job.)
Note: As i also said in the comment, the fact that $40!^{1/3}$ is so close to the number $999\dots 9$ with $16$ digits makes it harder to get such decompositions with three factors of $16,16,16$ digits each. Because the range of the primes is not really between the $16$ digits numbers $1000\dots 0$ and $999\dots 9$. The computer was founding for instance also
k= 33377 40!+k = 1056334373719901 * 22413347482811693 * 34461718591042889
k= 42269 40!+k = 1859373333290887 * 2719883435120179 * 161334845671400153
k= 47189 40!+k = 1822739844251111 * 15851904710438411 * 28238324719193609
k= 55049 40!+k = 1101465899839957 * 3378301824756127 * 219268155834100091
with factors having $\ge 16$ digits,
and also insisting that the smaller prime factor among all three of the three is bigger than $7000\dots 0$ ($16$ digits)...
k=300722963 40!+k = 7073009263847467 * 8126384986908847 * 14195263159044887
k=349245887 40!+k = 7209029739646403 * 8613116168531887 * 13140380656971067
k= 77639477 40!+k = 7252468942583363 * 7992027613367573 * 14076744068584523
k=300731429 40!+k = 7321097013439289 * 7404841114350401 * 15050577417998861
k=339499973 40!+k = 7398940260351053 * 10008305130239719 * 11018310576981439
k= 96999373 40!+k = 7532964349652891 * 9232067350535819 * 11732219831109637
k= 48531323 40!+k = 7564931816003599 * 10225749107742383 * 10547387873068219
k=465611681 40!+k = 7866432514968563 * 7912332084219151 * 13108794049712437
k=465608309 40!+k = 8221566296360779 * 9835813811956651 * 10089745359967421
k=126205139 40!+k = 8278202201940731 * 9790182376796351 * 10067421590328119
k=388112077 40!+k = 8601556182258047 * 8687613919887803 * 10918614312712297
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
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).