lgamma->gammaln

[Paul Kienzle]

Paul Kienzle <address@hidden>

        * miscellaneous/bincoeff.m: Replace lgamma with gammaln.
        * specfun/beta.m: Likewise.
        * special-matrix/invhilb.m: Likewise (but it is only in a comment).
        * statistics/distributions/gamma_pdf.m: Likewise.
        * statistics/distributions/poisson_pdf.m: Likewise.

Index: miscellaneous/bincoeff.m
===================================================================
RCS file: /cvs/octave/scripts/miscellaneous/bincoeff.m,v
retrieving revision 1.5
diff -c -r1.5 bincoeff.m
*** miscellaneous/bincoeff.m    2002/05/01 06:48:35     1.5
--- miscellaneous/bincoeff.m    2002/08/09 17:29:49
***************
*** 83,96 ****
  
    ind = find ((k > 0) & ((n == real (round (n))) & (n < 0)));
    if any (ind)
!     b(ind) = (-1) .^ k(ind) .* exp (lgamma (abs (n(ind)) + k(ind)) ...
!         - lgamma (k(ind) + 1) - lgamma (abs (n(ind))));
    endif
  
    ind = find ((k > 0) & ((n != real (round (n))) | (n >= k)));
    if (length (ind) > 0)
!     b(ind) = exp (lgamma (n(ind) + 1) - lgamma (k(ind) + 1) ...
!         - lgamma (n(ind) - k(ind) + 1));
    endif
  
    ## clean up rounding errors
--- 83,96 ----
  
    ind = find ((k > 0) & ((n == real (round (n))) & (n < 0)));
    if any (ind)
!     b(ind) = (-1) .^ k(ind) .* exp (gammaln (abs (n(ind)) + k(ind)) ...
!         - gammaln (k(ind) + 1) - gammaln (abs (n(ind))));
    endif
  
    ind = find ((k > 0) & ((n != real (round (n))) | (n >= k)));
    if (length (ind) > 0)
!     b(ind) = exp (gammaln (n(ind) + 1) - gammaln (k(ind) + 1) ...
!         - gammaln (n(ind) - k(ind) + 1));
    endif
  
    ## clean up rounding errors
Index: specfun/beta.m
===================================================================
RCS file: /cvs/octave/scripts/specfun/beta.m,v
retrieving revision 1.15
diff -c -r1.15 beta.m
*** specfun/beta.m      2000/01/13 08:40:27     1.15
--- specfun/beta.m      2002/08/09 17:29:49
***************
*** 45,50 ****
      usage ("beta (a, b)");
    endif
  
!   retval = exp (lgamma (a) + lgamma (b) - lgamma (a+b));
  
  endfunction
--- 45,50 ----
      usage ("beta (a, b)");
    endif
  
!   retval = exp (gammaln (a) + gammaln (b) - gammaln (a+b));
  
  endfunction
Index: special-matrix/invhilb.m
===================================================================
RCS file: /cvs/octave/scripts/special-matrix/invhilb.m,v
retrieving revision 1.16
diff -c -r1.16 invhilb.m
*** special-matrix/invhilb.m    2002/04/04 23:03:15     1.16
--- special-matrix/invhilb.m    2002/08/09 17:29:49
***************
*** 85,91 ****
      ## machine number, the result is also exact.  Otherwise we calculate
      ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)).
      ##
!     ## The Octave bincoeff routine uses transcendental functions (lgamma
      ## and exp) rather than multiplications, for the sake of speed.  
      ## However, it rounds the answer to the nearest integer, which 
      ## justifies the claim about exactness made above.
--- 85,91 ----
      ## machine number, the result is also exact.  Otherwise we calculate
      ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)).
      ##
!     ## The Octave bincoeff routine uses transcendental functions (gammaln
      ## and exp) rather than multiplications, for the sake of speed.  
      ## However, it rounds the answer to the nearest integer, which 
      ## justifies the claim about exactness made above.
Index: statistics/distributions/gamma_pdf.m
===================================================================
RCS file: /cvs/octave/scripts/statistics/distributions/gamma_pdf.m,v
retrieving revision 1.5
diff -c -r1.5 gamma_pdf.m
*** statistics/distributions/gamma_pdf.m        2002/06/27 14:14:08     1.5
--- statistics/distributions/gamma_pdf.m        2002/08/09 17:29:49
***************
*** 59,65 ****
    k = find ((x > 0) & (a > 1) & (b > 0));
    if (any (k))
      pdf(k) = exp (a(k) .* log (b(k)) + (a(k)-1) .* log (x(k))
!                 - b(k) .* x(k) - lgamma (a(k)));
    endif
  
    pdf = reshape (pdf, r, c);
--- 59,65 ----
    k = find ((x > 0) & (a > 1) & (b > 0));
    if (any (k))
      pdf(k) = exp (a(k) .* log (b(k)) + (a(k)-1) .* log (x(k))
!                 - b(k) .* x(k) - gammaln (a(k)));
    endif
  
    pdf = reshape (pdf, r, c);
Index: statistics/distributions/poisson_pdf.m
===================================================================
RCS file: /cvs/octave/scripts/statistics/distributions/poisson_pdf.m,v
retrieving revision 1.4
diff -c -r1.4 poisson_pdf.m
*** statistics/distributions/poisson_pdf.m      2002/05/01 06:48:36     1.4
--- statistics/distributions/poisson_pdf.m      2002/08/09 17:29:49
***************
*** 50,56 ****
  
    k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (l > 0));
    if (any (k))
!     pdf(k) = exp (x(k) .* log (l(k)) - l(k) - lgamma (x(k) + 1));
    endif
  
    pdf = reshape (pdf, r, c);
--- 50,56 ----
  
    k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (l > 0));
    if (any (k))
!     pdf(k) = exp (x(k) .* log (l(k)) - l(k) - gammaln (x(k) + 1));
    endif
  
    pdf = reshape (pdf, r, c);