Here is the code i tried:

```
S = 553276188 # 553.276.187
P = prod( list( primes(20) ) )
B = floor(S/P) + 1 # 58
F = factorial(40)
LIST = [ a for a in xrange(P) if gcd(a, P) == 1 ]
P1 = prod( list( primes(20, 40) ) )
P2 = prod( list( primes(40, 80) ) )
P3 = prod( list( primes(80, 160) ) )
counter = 0
for a in LIST:
if a > counter:
print "a=%s ... still searching" % a
counter += 100
for b in range(B):
# 58 is related to S/P
k = b*P + a
n = F + k
if gcd(n, P1) > 1 or gcd(n, P2) > 1 or gcd(n, P3) > 1:
continue
factorization = n.factor()
if len(factorization) != 3:
continue
(f1, mul1), (f2, mul2), (f3, mul3) = factorization
if len(str(f1)) < 16:
continue
print "a=%6s b=%2s k=%9s 40!+k = %16s * %16s * %s" % ( a, b, k, f1, f2, f3 )
if ( mul1 == 1 and
mul2 == 1 and
mul3 == 1 and
f1 > 10^15 and
f3 < 10^16):
print 100*'='
print "SOLUTION !!! "
print "k = %s" % k
print "40! + k = %s" % factorization
print "1.st factor is %s" % f1
print "2.nd factor is %s" % f2
print "3.rd factor is %s" % f3
print 100*'='
print
```

It it not optimal, uses only some obvious ways to eliminate $k$ values in
$$
40!+k
$$
that have no chance to work. For instance, the $k$ divisible by one of the primes till $40$. This is done as follows.
Consider the product `P`

of all primes up to $20$. Then eliminate all numbers with a gcd bigger one w.r.t. `P`

up to `P`

. This delivers a list with $\phi(9699690)$ elements,

```
sage: euler_phi(P)
1658880
sage: ( euler_phi(P) / P ) . n()
0.171024022417211
```

so we consider less then one fifth of the numbers $a$ in the interval `[ 1..P ]`

. For each such $a$ consider all $k$--values up to the known solution $S=553276187$, which satisfy $k\equiv a$ modulo $P$. (This new enumeration of the relevant $k$ values in the interval was trying to do "something else" then taking them one by one, since the message in the OP of MSE was there is no solution in the first two mil. values.)

I had to stop the run after some(two) days. (Life is too short.) The last computed lines were:

```
a=125401 ... still searching
a=125483 b=44 k=426911843 40!+k = 3990918283471681 * 11226713503201489 * 18210404503813427
a=125501 ... still searching
a=125603 ... still searching
a=125687 b=20 k=194119487 40!+k = 2150616161684849 * 2660769607977791 * 142585345306056593
a=125701 ... still searching
a=125803 ... still searching
a=125863 b=54 k=523909123 40!+k = 3542386248615503 * 6801209511973739 * 33865929607036519
and the only solution found so far was extracted from the lines
a=120203 ... still searching
a=120283 b=51 k=494804473 40!+k = 8912658466556113 * 9232286052422441 * 9915818101022081
====================================================================================================
SOLUTION !!!
k = %s
40! + k = 8912658466556113 * 9232286052422441 * 9915818101022081
1.st factor is 8912658466556113
2.nd factor is 9232286052422441
3.rd factor is 9915818101022081
====================================================================================================
a=120301 ... still searching
```

As seen above, there was an error in the running code, that `k = %s`

should have given the solution. Fortunately i had it in the previous line.

The "obvious elimination" did after was to first impose also the condition that $40!+k$ has no small primes inside. I did this in three tranches, imposing the conditions that the gcd of $40!+k$ and `P1`

is one. Well, before factorizing, i decided to also use `P2`

, `P3`

in a similar way, since statistically (for the eye and my feeling) i saw too many "small primes". And only after this pass to the full factorization. Where we reject all factorization with (strictly) more or less than three factors.

For a factorization with exactly three factors, the smaller one, `f1`

, should have $16$ digits. I am printing each such constellation.

Note: It is hard to find a solution, since $40!^{1/3}$ is very close to $\underbrace{999...9}_{16\text{ digits}}$, it is:

```
sage: print '%19.2f' % factorial(40)^(1./3.)
9344334059477234.00
```

so at least one prime factor among the three, `f3`

, the bigger one is between this number and `999...999`

. But then, for such a fixed prime `f1`

, the other two prime numbers, `f1`

, `f2`

, also do not have to a too big range to live in, so that the product is "very close to" $40!$.

This is all i have so far.

Homework ?

@Emmanuel Charpentier no it is a project in mathexchange.

Also posted at https://stackoverflow.com/questions/5...