Mercurial > hg > octave-lojdl
changeset 7601:8a939b217863
Treat negative values to lgamma and beta correctly
author | David Bateman <dbateman@free.fr> |
---|---|
date | Tue, 18 Mar 2008 21:32:48 -0400 |
parents | 24abf5a702d9 |
children | 7bfaa9611558 |
files | ChangeLog configure.in liboctave/ChangeLog liboctave/lo-specfun.cc liboctave/lo-specfun.h liboctave/randpoisson.c scripts/ChangeLog scripts/specfun/beta.m src/ChangeLog src/mappers.cc src/ov-re-mat.cc src/ov-re-sparse.cc src/ov-scalar.cc |
diffstat | 13 files changed, 81 insertions(+), 13 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -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.
--- a/configure.in +++ b/configure.in @@ -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)
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -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.
--- a/liboctave/lo-specfun.cc +++ b/liboctave/lo-specfun.cc @@ -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);
--- a/liboctave/lo-specfun.h +++ b/liboctave/lo-specfun.h @@ -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);
--- a/liboctave/randpoisson.c +++ b/liboctave/randpoisson.c @@ -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; }
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -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.
--- a/scripts/specfun/beta.m +++ b/scripts/specfun/beta.m @@ -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)
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -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>
--- a/src/mappers.cc +++ b/src/mappers.cc @@ -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") {
--- a/src/ov-re-mat.cc +++ b/src/ov-re-mat.cc @@ -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)
--- a/src/ov-re-sparse.cc +++ b/src/ov-re-sparse.cc @@ -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)
--- a/src/ov-scalar.cc +++ b/src/ov-scalar.cc @@ -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)