\documentclass{article}
\usepackage{expl3}

\ExplSyntaxOn

\cs_new:Npn \fp_factorial:n #1
 % computes #1!
 {
  % TODO: signal error for #1 < 0
  \fp_compare:nTF { (#1) == 0 }
   { 1 }
   { \fp_eval:n { (#1) * \fp_factorial:n { \fp_eval:n { (#1) - 1 } } } }
 }

\cs_new:Npn \fp_gamma:n #1
 % computer gamma(#1)
 {
  \fp_compare:nTF { (#1) < 0.5 }
   { \__fp_gamma_reflection:n { \fp_eval:n { #1 } } }
   { \__fp_gamma_lanczos:n { \fp_eval:n { #1 - 1 } } }
 }

\cs_new:Npn \__fp_gamma_reflection:n #1
 % #1 = x
 {
  \fp_eval:n
   {
    pi / ( sin( pi * #1 ) * \fp_gamma:n { \fp_eval:n { 1 - #1 } } )
   }
 }

% https://en.wikipedia.org/wiki/Lanczos_approximation
\cs_new:Npn \__fp_gamma_lanczos:n #1
 % #1 = x-1
 {
  \fp_eval:n
   {
    sqrt(2*pi) * (#1 + 7 + 0.5)^(#1+0.5) *
    exp(-(#1 + 7 + 0.5)) *
    ( 0.99999999999980993
    +  676.5203681218851    / (#1 + 1)
    - 1259.1392167224028    / (#1 + 2)
    + 771.32342877765313    / (#1 + 3)
    - 176.61502916214059    / (#1 + 4)
    + 12.507343278686905    / (#1 + 5)
    - 0.13857109526572012   / (#1 + 6)
    + 9.9843695780195716e-6 / (#1 + 7)
    + 1.5056327351493116e-7 / (#1 + 8) )
   }
 }

% Macros for tests
\msg_new:nnn { fp } { test-failed } { Test ~ failed! ~ #1 }
\cs_new:Npn \test:nn #1#2
 {
  \edef\temp{#1}
  \fp_compare:nF { (\temp < #2 + 1e-11) && (\temp > #2 - 1e-11) }
   {
    \exp_args:Nnnx
    \msg_error:nnn { fp } { test-failed }
     { \tl_to_str:n { #1 } ~ = ~ \temp\space != ~ #2 }
   }
 }

% Tests for factorial
\test:nn { \fp_factorial:n { 0 } } { 1 }
\test:nn { \fp_factorial:n { 8 } } { 40320 }

% Tests for gamma
\test:nn { \fp_gamma:n { -0.5 } } {   -3.5449077018110295 }
\test:nn { \fp_gamma:n {  1   } } {    1                  }
\test:nn { \fp_gamma:n {  2.5 } } {    1.3293403881791384 }
\test:nn { \fp_gamma:n {  8   } } { 5040.000000000011     }

\ExplSyntaxOff

\begin{document}
\end{document}
