fix(factorial): result should be of the default kind, not double default kind

[cyrilgandon]

Fixes #875

The parameter 1.0D0 is promoted to a real kind 16 when use in conjunction with the flag of gfortran default-real-8, and there is not such an overload for the l_gamma function.
Using a double precision number seems a bug here, nowhere else in the stdlib we have this usage.

Especially when the result is defined as a default real:

real :: res

I don't see a reason to pass a double precision parameter to the l_gamma function. I propose to drop the double precision here.

Edit: here is a comment by OP

#625 (comment)

By using higher precision during calculation, the desired precision digits are guaranteed to be correct. So single precision uses double precision for calculation, double precision uses quadruple precision for calculation. Extended precision calculation may need higher precision a little over quadruple, while quadruple precision calculation need much higher precision not available right now.

So I'm not sure anymore that I should drop that double precision.
Maybe @Jim-215-Fisher can help here. What should happen if default real size is 64bit, should we still do calculation with 128bits real? I think not since most of architecture are not 128bits..

I'd say that double precision logic is included is the gamma functions, no need to have that at the level of the l_factorial.

[perazz]

LGTM, this looks like the right thing to do here. Thank you @cyrilgandon

[cyrilgandon]

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

[perazz]

@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.

[cyrilgandon]

@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.

[perazz]

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).

[cyrilgandon]

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

[perazz]

LGTM. @jvdp1 what do you think?

[jvdp1]

LGTM. Thank you @cyrilgandon

[jvdp1]

Thank you @cyrilgandon . I will merge this PR.