MIME-Version: 1.0
References: <CAD7iri+qhPKJPu5w4=Ku=C3+8G8f3hY5uceHXMPfBfxrD_aOoA@mail.gmail.com>
Content-Type: text/plain; charset=UTF-8
Message-ID:  <CANQYN6xKKZ9QTSxW3XLoGi_faPDCiW6mCWCqXofkZgQwXtcc+g@mail.gmail.com>
Date:         Fri, 15 Jan 2016 09:31:04 -0500
Reply-To: Mailing list for the LaTeX3 project
              <LATEX-L@LISTSERV.UNI-HEIDELBERG.DE>
Sender: Mailing list for the LaTeX3 project <LATEX-L@LISTSERV.UNI-HEIDELBERG.DE>
From: Bruno Le Floch <blflatex@GMAIL.COM>
Subject: Re: Gamma function and factorial
To: LATEX-L@LISTSERV.UNI-HEIDELBERG.DE
In-Reply-To:  <CAD7iri+qhPKJPu5w4=Ku=C3+8G8f3hY5uceHXMPfBfxrD_aOoA@mail.gmail.com>
Precedence: list
Envelope-To: <rainer.schoepf@GMX.NET>
Status: R

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;
  }