After reading attentively the code, the 2nd parameter of the l_gamma function is very important and use only for it's kind.
So it is in fact a huge change to remove the doubling.

Here is some test on log(10!):

Original code with doubling            : 15.1044130325
PR code removing doubling            : 15.1044120789
print (math.log(math.factorial(10))) : 15.104412573075516

@cyrilgandon, this is very interesting. It looks like the error here is close to epsilon(0.0), so, perhaps it would be a good idea to shift all real variables to real(dp)? So, the internal accuracy would be higher, and anyways independent of the compiler's default real choice.

@perazz Yes!! I think it's an easy fix and should answer the problem. For most people having the default real(4), it will not change the result.

If I read the OP's comment well, a double precision version (integer(int64) input as far as I understand) may benefit from a quad-precision internal calculation.
However, the real accuracy has not been templated here, so I believe it is good enough to promote it to a fixed 64-bit precision, like this PR is suggesting. Maybe, it would make sense to spread the same precision to all real variables, including the function result. That would not introduce a backward incompatible change, as a function result can be automatically converted (and the new version provides higher precision).

Yes that makes sense, also rename zero_k2 to zero_dp, as the k2 paramter did not exist in this function