Sage program for 40!+k project

edit

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.