# HG changeset patch
# User David Bateman <dbateman@free.fr>
# Date 1205890368 14400
# Node ID 8a939b21786329cab992c0e1136a4c7e90d0128a
# Parent  24abf5a702d93057c7a788db2a98a624c1eb5802
Treat negative values to lgamma and beta correctly

diff -r 24abf5a702d9 -r 8a939b217863 ChangeLog
--- a/ChangeLog	Tue Mar 18 20:27:50 2008 -0400
+++ b/ChangeLog	Tue Mar 18 21:32:48 2008 -0400
@@ -1,3 +1,7 @@
+2008-03-18  David Bateman  <dbateman@free.fr>
+
+	* configure.in (AC_CHECK_FUNCS): Also check lgamma_r.
+
 2008-03-11  David Bateman  <dbateman@free.fr>
 
 	* run-octave.in: Fix typo.
diff -r 24abf5a702d9 -r 8a939b217863 configure.in
--- a/configure.in	Tue Mar 18 20:27:50 2008 -0400
+++ b/configure.in	Tue Mar 18 21:32:48 2008 -0400
@@ -1297,11 +1297,11 @@
 AC_CHECK_FUNCS(atexit basename bcopy bzero canonicalize_file_name chmod \
   dup2 endgrent endpwent execvp fcntl fork getcwd getegid geteuid \
   getgid getgrent getgrgid getgrnam getpgrp getpid getppid getpwent \
-  getpwuid gettimeofday getuid getwd _kbhit kill lgamma link localtime_r \
-  lstat memmove mkdir mkfifo mkstemp on_exit pipe poll putenv raise \
-  readlink realpath rename resolvepath rindex rmdir round select setgrent \
-  setlocale setpwent setvbuf sigaction siglongjmp sigpending sigprocmask \
-  sigsuspend snprintf stat strcasecmp strdup strerror stricmp \
+  getpwuid gettimeofday getuid getwd _kbhit kill lgamma lgamma_r link \
+  localtime_r lstat memmove mkdir mkfifo mkstemp on_exit pipe poll putenv \
+  raise readlink realpath rename resolvepath rindex rmdir round select \
+  setgrent setlocale setpwent setvbuf sigaction siglongjmp sigpending \
+  sigprocmask sigsuspend snprintf stat strcasecmp strdup strerror stricmp \
   strncasecmp strnicmp strptime strsignal symlink tempnam tgamma umask \
   uname unlink usleep utime vfprintf vsprintf vsnprintf waitpid \
   _chmod _snprintf x_utime _utime32)
diff -r 24abf5a702d9 -r 8a939b217863 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Tue Mar 18 20:27:50 2008 -0400
+++ b/liboctave/ChangeLog	Tue Mar 18 21:32:48 2008 -0400
@@ -1,5 +1,9 @@
 2008-03-18  David Bateman  <dbateman@free.fr>
 
+	* lo-specfun.cc (Complex xlgamma (const Complex&)): New function.
+	* lo-specfun.h (Complex xlgamma (const Complex&)): Declare it.
+	* randpoison.c (xlgamma): Use lgamma if HAVE_LGAMMA is defined.
+	
 	* dNDArray.cc (NDArray::min, NDArraymax): chop trailing singletons.
 	* CNDarray.cc (ComplexNDArray::min, CompelxNDArray::max): ditto.
 	* intNDarray.cc (intNDArray<T>::min, intNDArray<T>::max): ditto.
diff -r 24abf5a702d9 -r 8a939b217863 liboctave/lo-specfun.cc
--- a/liboctave/lo-specfun.cc	Tue Mar 18 20:27:50 2008 -0400
+++ b/liboctave/lo-specfun.cc	Tue Mar 18 21:32:48 2008 -0400
@@ -196,6 +196,34 @@
 #endif
 }
 
+Complex
+xlgamma (const Complex& xc)
+{
+  // Can only be called with a real value of x.
+  double x = xc.real ();
+  double result;
+
+#if defined (HAVE_LGAMMA_R)
+  int sgngam;
+  result = lgamma_r (x, &sgngam);    
+#else
+  double sgngam;
+
+  if (xisnan (x))
+    result = x;
+  else if (xisinf (x))
+    result = octave_Inf;
+  else
+    F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));
+
+#endif
+
+  if (sgngam < 0)
+    return result + Complex (0., M_PI);
+  else
+    return result;
+}
+
 static inline Complex
 zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
 
diff -r 24abf5a702d9 -r 8a939b217863 liboctave/lo-specfun.h
--- a/liboctave/lo-specfun.h	Tue Mar 18 20:27:50 2008 -0400
+++ b/liboctave/lo-specfun.h	Tue Mar 18 21:32:48 2008 -0400
@@ -59,6 +59,7 @@
 
 extern OCTAVE_API double xgamma (double x);
 extern OCTAVE_API double xlgamma (double x);
+extern OCTAVE_API Complex xlgamma (const Complex& x);
 
 extern OCTAVE_API Complex
 besselj (double alpha, const Complex& x, bool scaled, octave_idx_type& ierr);
diff -r 24abf5a702d9 -r 8a939b217863 liboctave/randpoisson.c
--- a/liboctave/randpoisson.c	Tue Mar 18 20:27:50 2008 -0400
+++ b/liboctave/randpoisson.c	Tue Mar 18 21:32:48 2008 -0400
@@ -59,6 +59,9 @@
 xlgamma (double x)
 {
   double result;
+#ifdef HAVE_LGAMMA
+  result = lgamma (x);
+#else
   double sgngam;
 
   if (lo_ieee_isnan (x))
@@ -67,7 +70,7 @@
     result = octave_Inf;
   else
     F77_XFCN (dlgams, DLGAMS, (&x, &result, &sgngam));
-
+#endif
   return result;
 }
 
diff -r 24abf5a702d9 -r 8a939b217863 scripts/ChangeLog
--- a/scripts/ChangeLog	Tue Mar 18 20:27:50 2008 -0400
+++ b/scripts/ChangeLog	Tue Mar 18 21:32:48 2008 -0400
@@ -1,3 +1,7 @@
+2008-03-18  Ben Abbott <bpabbott@mac.com>
+
+	* specfun/beta.m: Fix for negative inputs.
+
 2008-03-18  Michael D. Godfrey  <godfrey@isl.stanford.edu>
 
 	* plot/__go_draw_axes__.m: Use correct symbol codes.
diff -r 24abf5a702d9 -r 8a939b217863 scripts/specfun/beta.m
--- a/scripts/specfun/beta.m	Tue Mar 18 20:27:50 2008 -0400
+++ b/scripts/specfun/beta.m	Tue Mar 18 21:32:48 2008 -0400
@@ -19,7 +19,7 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Mapping Function} {} beta (@var{a}, @var{b})
-## Return the Beta function,
+## For real inputs, return the Beta function,
 ## @iftex
 ## @tex
 ## $$
@@ -45,7 +45,15 @@
     print_usage ();
   endif
 
-  retval = exp (gammaln (a) + gammaln (b) - gammaln (a+b));
+  if (any (size (a) != size (b)) && numel (a) != 1 && numel (b) != 1)
+    error ("beta: inputs have inconsistent sizes.")
+  endif
+
+  if (! isreal (a) || ! isreal (b))
+    error ("beta: inputs must be real.")
+  endif
+
+  retval = real (exp (gammaln (a) + gammaln (b) - gammaln (a+b)));
 
 endfunction
 
@@ -61,3 +69,16 @@
 
 %!error beta(1);
 
+%!assert (1, beta (1, 1))
+
+%!test
+%! a = 2:10;
+%! tol = 10 * max (a) * eps;
+%! assert (-a, beta (-1./a, 1), tol)
+%! assert (-a, beta (1, -1./a), tol)
+
+%!test
+%! a = 0.25 + (0:5) * 0.5;
+%! tol = 10 * max (a) * eps;
+%! assert (zeros (size (a)), beta (a, -a), tol)
+%! assert (zeros (size (a)), beta (-a, a), tol)
diff -r 24abf5a702d9 -r 8a939b217863 src/ChangeLog
--- a/src/ChangeLog	Tue Mar 18 20:27:50 2008 -0400
+++ b/src/ChangeLog	Tue Mar 18 21:32:48 2008 -0400
@@ -1,5 +1,9 @@
 2008-03-18  David Bateman  <dbateman@free.fr>
 
+	* ov-re-mat.cc (lgamma): Convert to a allow negative arguments.
+	* ov-re-sparse.cc (lgamma): ditto.
+	* ov-scalar.cc (lgamma): ditto.
+
 	* DLD-FUNCTIONS/minmax.cc: 64-bit indexing fix.
 
 2008-03-13  John W. Eaton  <jwe@octave.org>
diff -r 24abf5a702d9 -r 8a939b217863 src/mappers.cc
--- a/src/mappers.cc	Tue Mar 18 20:27:50 2008 -0400
+++ b/src/mappers.cc	Tue Mar 18 21:32:48 2008 -0400
@@ -725,8 +725,7 @@
     "-*- texinfo -*-\n\
 @deftypefn {Mapping Function} {} lgamma (@var{x})\n\
 @deftypefnx {Mapping Function} {} gammaln (@var{x})\n\
-Return the natural logarithm of the absolute value of the gamma\n\
-function of @var{x}.\n\
+Return the natural logarithm of the gamma function of @var{x}.\n\
 @seealso{gamma, gammai}\n\
 @end deftypefn")
 {
diff -r 24abf5a702d9 -r 8a939b217863 src/ov-re-mat.cc
--- a/src/ov-re-mat.cc	Tue Mar 18 20:27:50 2008 -0400
+++ b/src/ov-re-mat.cc	Tue Mar 18 21:32:48 2008 -0400
@@ -715,7 +715,7 @@
 ARRAY_MAPPER (erf, NDArray::dmapper, ::erf)
 ARRAY_MAPPER (erfc, NDArray::dmapper, ::erfc)
 ARRAY_MAPPER (gamma, NDArray::dmapper, xgamma)
-ARRAY_MAPPER (lgamma, NDArray::dmapper, xlgamma)
+CD_ARRAY_MAPPER (lgamma, xlgamma, xlgamma, 0.0, octave_Inf)
 ARRAY_MAPPER (abs, NDArray::dmapper, ::fabs)
 ARRAY_MAPPER (acos, NDArray::dmapper, ::acos)
 CD_ARRAY_MAPPER (acosh, ::acosh, ::acosh, 1.0, octave_Inf)
diff -r 24abf5a702d9 -r 8a939b217863 src/ov-re-sparse.cc
--- a/src/ov-re-sparse.cc	Tue Mar 18 20:27:50 2008 -0400
+++ b/src/ov-re-sparse.cc	Tue Mar 18 21:32:48 2008 -0400
@@ -905,7 +905,7 @@
 SPARSE_MAPPER (erf, SparseMatrix::dmapper, ::erf)
 SPARSE_MAPPER (erfc, SparseMatrix::dmapper, ::erfc)
 SPARSE_MAPPER (gamma, SparseMatrix::dmapper, xgamma)
-SPARSE_MAPPER (lgamma, SparseMatrix::dmapper, xlgamma)
+CD_SPARSE_MAPPER (lgamma, xlgamma, xlgamma, 0.0, octave_Inf)
 SPARSE_MAPPER (abs, SparseMatrix::dmapper, ::fabs)
 SPARSE_MAPPER (acos, SparseMatrix::dmapper, ::acos)
 CD_SPARSE_MAPPER (acosh, ::acosh, ::acosh, 1.0, octave_Inf)
diff -r 24abf5a702d9 -r 8a939b217863 src/ov-scalar.cc
--- a/src/ov-scalar.cc	Tue Mar 18 20:27:50 2008 -0400
+++ b/src/ov-scalar.cc	Tue Mar 18 21:32:48 2008 -0400
@@ -310,7 +310,7 @@
 SCALAR_MAPPER (erf, ::erf)
 SCALAR_MAPPER (erfc, ::erfc)
 SCALAR_MAPPER (gamma, xgamma)
-SCALAR_MAPPER (lgamma, xlgamma)
+CD_SCALAR_MAPPER (lgamma, xlgamma, xlgamma, 0.0, octave_Inf)
 SCALAR_MAPPER (abs, ::fabs)
 SCALAR_MAPPER (acos, ::acos)
 CD_SCALAR_MAPPER (acosh, ::acosh, ::acosh, 1.0, octave_Inf)