Dear Henri,
Thank you for your contribution. I've logged an issue at
https://github.com/latex3/latex3/issues/297 , without including your
name since I didn't know if you wanted it public. Since it's about
l3fp I'll be the one to work on that issue. Unfortunately I'll be
overworked in the next month or so, so I cannot get back to that until
at least midFebruary, or more probably early March.
> I have implemented a preliminary version of
> \fp_factorial:n and \fp_gamma:n
> The current working precision of \fp_gamma:n is limited to 1e11.
I'm a bit worried about precision, as I've worked quite hard to make
all other functions have full 1e16 precision. Internally these other
functions use extended precision for intermediate calculations (24
digit fixedpoint numbers in [0,1), implemented in l3fpextended:
https://github.com/latex3/latex3/blob/master/l3kernel/l3fpextended.dtx
). It would be interesting if you could find out whether this
precision would be enough to reach 16 digits of precision in all cases
using a suitably extended Lanczos approximation. I'm particularly
interested in a proof that this is the case. It might be better to
use Spouge's approximation
https://en.wikipedia.org/wiki/Spouge's_approximation as coefficients
appear to be easier to compute.
Another can of worms is to be certain that the function is exactly
monotonous where it should be. E.g. I recall having problems with log
in the regard, where it would jump by 1e16 at junctions between
different intervals used in the calculation. I need to look again at
log.
> So far, error checking for the argument of factorial is missing and I
> don't know how to integrate this into l3fp, especially making these
> available in \fp_eval:n as factorial(x) and gamma(x).
>
> If anyone would be so kind as to provide some guidance, I would be
> very grateful.
That would be easy for me as I'm used to the codebase, but a little
bit tricky to explain. See below. Here are the steps that need to be
done on this issue:
 Find a good approximation (1e16) of the Gamma function which can be
computed using 24digit fixedpoint numbers. For speed it is best to
limit as much as possible the number of divisions; besides, I don't
remember if full 24digit precision is available or if it produces
only 16 correct digits. I can look if needed.
 Implement it using l3fpextended (I may need to do that; if you are
highly motivated you can take a look at implementation of sin in
l3fptrig.dtx).
 Implement error checking (ditto, a matter of an hour).
 Teach it to l3fpparse (that's a matter of minutes).
Best regards,
Bruno

For functions currently in l3fp, there are two parts: one is to teach
them to the parsing step (in l3fpparse.dtx): e.g. for asecd that is
equivalent to doing
\cs_new_nopar:Npn \@@_parse_word_asecd:N
{ \@@_parse_unary_function:nNN {asec} \use_ii:nn }
The next is to define the macro \@@_asec_o:w (defined in
l3fptrig.dtx), which ends up being called followed by \use_ii:nn
(because of the def above) and an internal floating point and @.
Treating all special cases gets messy... The following code says that
if #2 (the "type" of floating point) is 0 (floating point is +0 or 0)
then we have an invalid operation since 0 is outside the domain of
arcsecant; if #2 is 1 (normal floating point) then we call an internal
function implementing acsc (and reverse arguments); if #2 is 2
(positive or negative infinity) then we call an internal function
shared with atan; if #3 is 3 (NaN), just leave the NaN be.
\cs_new:Npn \@@_asec_o:w #1 \s_@@ \@@_chk:w #2#3; @
{
\if_case:w #2 \exp_stop_f:
\@@_case_use:nw
{ \@@_invalid_operation_o:fw { #1 { asec } { asecd } } }
\or:
\@@_case_use:nw
{
\@@_acsc_normal_o:NfwNnw #1 { #1 { asec } { asecd } }
\@@_reverse_args:Nww
}
\or: \@@_case_use:nw { \@@_atan_inf_o:NNNw #1 0 \c_four }
\else: \@@_case_return_same_o:w
\fi:
\s_@@ \@@_chk:w #2 #3;
}
