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 mid-February, 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 1e-11. I'm a bit worried about precision, as I've worked quite hard to make all other functions have full 1e-16 precision. Internally these other functions use extended precision for intermediate calculations (24 digit fixed-point numbers in [0,1), implemented in l3fp-extended: https://github.com/latex3/latex3/blob/master/l3kernel/l3fp-extended.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 1e-16 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 code-base, 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 (1e-16) of the Gamma function which can be computed using 24-digit fixed-point numbers. For speed it is best to limit as much as possible the number of divisions; besides, I don't remember if full 24-digit precision is available or if it produces only 16 correct digits. I can look if needed. - Implement it using l3fp-extended (I may need to do that; if you are highly motivated you can take a look at implementation of sin in l3fp-trig.dtx). - Implement error checking (ditto, a matter of an hour). - Teach it to l3fp-parse (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 l3fp-parse.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 l3fp-trig.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; }