From 1ba89195c451e9eacea9a1d9f6104806dbe17887 Mon Sep 17 00:00:00 2001
From: Taylor R Campbell <campbell@mumble.net>
Date: Sun, 30 Jun 2019 23:28:06 +0000
Subject: [PATCH] Add flo:signed-lgamma.

---
 doc/ref-manual/numbers.texi          | 35 ++++++++++++++++++++++++++++
 src/microcode/flonum.c               | 20 ++++++++++++++++
 src/runtime/primitive-arithmetic.scm |  6 ++++-
 src/runtime/runtime.pkg              |  1 +
 tests/runtime/test-arith.scm         | 10 +++++++-
 5 files changed, 70 insertions(+), 2 deletions(-)

diff --git a/doc/ref-manual/numbers.texi b/doc/ref-manual/numbers.texi
index fbd8b9f66..4975258bb 100644
--- a/doc/ref-manual/numbers.texi
+++ b/doc/ref-manual/numbers.texi
@@ -2036,6 +2036,41 @@ This is the flonum version of @code{atan} with two arguments.  When
 compiled, it does not check the types of its arguments.
 @end deffn
 
+@deffn procedure flo:signed-lgamma x
+Returns two values,
+@iftex
+@tex
+$$m = \log \left|\Gamma(x)\right| \quad\hbox{and}\quad s = \mathop{\rm sign} \Gamma(x),$$
+@end tex
+@end iftex
+@ifnottex
+
+@example
+@group
+m = log(|Gamma(@var{x})|)
+
+and
+
+s = sign(Gamma(@var{x})),
+@end group
+@end example
+
+@end ifnottex
+respectively a flonum and an exact integer either @code{-1} or
+@code{1}, so that
+@iftex
+@tex
+$$\Gamma(x) = s \cdot e^m.$$
+@end tex
+@end iftex
+@ifnottex
+
+@example
+Gamma(x) = s * e^m.
+@end example
+
+@end ifnottex
+@end deffn
 
 @deffn procedure flo:min x1 x2
 @deffnx procedure flo:max x1 x2
diff --git a/src/microcode/flonum.c b/src/microcode/flonum.c
index 806c4f697..b8c43fe0e 100644
--- a/src/microcode/flonum.c
+++ b/src/microcode/flonum.c
@@ -334,6 +334,26 @@ DEFINE_PRIMITIVE ("FLONUM-LGAMMA", Prim_flonum_lgamma, 1, 1, 0)
   FLONUM_RESULT (result);
 }
 
+DEFINE_PRIMITIVE ("FLONUM-SIGNED-LGAMMA", Prim_flonum_signed_lgamma, 1, 1, 0)
+{
+  double x;
+  double result;
+  int sign;
+  PRIMITIVE_HEADER (1);
+
+  x = (arg_flonum (1));
+#ifdef HAVE_LGAMMA_R
+  result = (lgamma_r (x, (&sign)));
+#else
+  result = (lgamma (x));
+  sign = signgam;
+#endif
+
+  assert (LONG_TO_FIXNUM_P (sign));
+  PRIMITIVE_RETURN
+    (cons ((double_to_flonum (result)), (LONG_TO_FIXNUM (sign))));
+}
+
 DEFINE_PRIMITIVE ("FLONUM-GAMMA", Prim_flonum_gamma, 1, 1, 0)
      SIMPLE_TRANSCENDENTAL_FUNCTION (tgamma)
 DEFINE_PRIMITIVE ("FLONUM-ERF", Prim_flonum_erf, 1, 1, 0)
diff --git a/src/runtime/primitive-arithmetic.scm b/src/runtime/primitive-arithmetic.scm
index 0c5a7c4fb..b535bf196 100644
--- a/src/runtime/primitive-arithmetic.scm
+++ b/src/runtime/primitive-arithmetic.scm
@@ -274,7 +274,7 @@ USA.
 
 (define (flo:total-order-mag x y)
   (flo:total-order (flo:abs x) (flo:abs y)))
-
+
 (define (flo:total< x y)
   (if (or (flo:nan? x) (flo:nan? y))
       ;; Must handle NaNs first and carefully to avoid exception on
@@ -420,6 +420,10 @@ USA.
   (if (and (flo:finite? x) (not (flo:safe-zero? x)))
       (fix:- (cdr ((ucode-primitive flonum-normalize 1) x)) 1)
       (begin (flo:raise-exceptions! (flo:exception:invalid-operation)) #f)))
+
+(define (flo:signed-lgamma x)
+  (let ((p ((ucode-primitive flonum-signed-lgamma 1) x)))
+    (values (car p) (cdr p))))
 
 ;;;; Exact integers
 
diff --git a/src/runtime/runtime.pkg b/src/runtime/runtime.pkg
index 145ed4b58..acb35ec28 100644
--- a/src/runtime/runtime.pkg
+++ b/src/runtime/runtime.pkg
@@ -359,6 +359,7 @@ USA.
 	  flo:safe>
 	  flo:safe>=
 	  flo:sign-negative?
+	  flo:signed-lgamma
 	  flo:sin
 	  flo:sinh
 	  flo:snan
diff --git a/tests/runtime/test-arith.scm b/tests/runtime/test-arith.scm
index 712bc68fb..5cc1f2370 100644
--- a/tests/runtime/test-arith.scm
+++ b/tests/runtime/test-arith.scm
@@ -1097,4 +1097,12 @@ USA.
         (assert-eqv (exact->inexact x) y)
         (if (not (= x y))
             (assert-except/no-traps (flo:exception:inexact-result)
-                                    (lambda () (exact->inexact x))))))))
\ No newline at end of file
+                                    (lambda () (exact->inexact x))))))))
+
+(define-enumerated-test 'flo:lgamma
+  (list (list -0.123))
+  (lambda (x)
+    (receive (log-gamma sign) (flo:signed-lgamma x)
+      (assert-eqv (flo:lgamma x) log-gamma)
+      (let ((gamma (* sign (exp log-gamma))))
+	(assert-<= (relerr (flo:gamma x) gamma) 1e-15)))))
\ No newline at end of file
-- 
2.25.1