Mercurial > hg > octave-avbm
changeset 4004:ca854fb51a88
[project @ 2002-07-25 06:31:33 by jwe]
author | jwe |
---|---|
date | Thu, 25 Jul 2002 06:32:37 +0000 |
parents | 54babc5054eb |
children | fc000ebb19df |
files | libcruft/ChangeLog libcruft/slatec-fn/xgmainc.f liboctave/ChangeLog liboctave/lo-specfun.cc liboctave/lo-specfun.h |
diffstat | 5 files changed, 167 insertions(+), 14 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,3 +1,7 @@ +2002-07-25 John W. Eaton <jwe@bevo.che.wisc.edu> + + * slatec-fn/xgmainc.f: New file. + 2002-07-12 John W. Eaton <jwe@bevo.che.wisc.edu> * dasrt: New subdirectory.
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/xgmainc.f @@ -0,0 +1,96 @@ + subroutine xgammainc (a, x, result) + +c -- jwe, based on DGAMIT. +c +c -- Do a better job than dgami for large values of x. + + double precision a, x, result + intrinsic exp, log, sqrt, sign, aint + + DOUBLE PRECISION AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX, + $ BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, DGAMR, D9GMIT, + $ D9LGIC, D9LGIT, DLNGAM, DGAMI + + LOGICAL FIRST + + SAVE ALNEPS, SQEPS, BOT, FIRST + + DATA FIRST /.TRUE./ + + if (x .eq. 0.0d0) then + + if (a .eq. 0.0d0) then + result = 1.0d0 + else + result = 0.0d0 + endif + + else + + IF (FIRST) THEN + ALNEPS = -LOG (D1MACH(3)) + SQEPS = SQRT(D1MACH(4)) + BOT = LOG (D1MACH(1)) + ENDIF + FIRST = .FALSE. +C + IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'XGMAINC', 'X IS NEGATIVE' + + , 2, 2) +C + IF (X.NE.0.D0) ALX = LOG (X) + SGA = 1.0D0 + IF (A.NE.0.D0) SGA = SIGN (1.0D0, A) + AINTA = AINT (A + 0.5D0*SGA) + AEPS = A - AINTA +C +C IF (X.GT.0.D0) GO TO 20 +C DGAMIT = 0.0D0 +C IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0) +C RETURN +C + 20 IF (X.GT.1.D0) GO TO 30 + IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1, + 1 SGNGAM) +C DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX) + result = exp (a*alx + log (D9GMIT (A, X, ALGAP1, SGNGAM, ALX))) + RETURN +C + 30 IF (A.LT.X) GO TO 40 + T = D9LGIT (A, X, DLNGAM(A+1.0D0)) + IF (T.LT.BOT) CALL XERCLR +C DGAMIT = EXP (T) + result = EXP (a*alx + T) + RETURN +C + 40 ALNG = D9LGIC (A, X, ALX) +C +C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X)) +C + H = 1.0D0 + IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50 +C + CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM) + T = LOG (ABS(A)) + ALNG - ALGAP1 + IF (T.GT.ALNEPS) GO TO 60 +C + IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T) + IF (ABS(H).GT.SQEPS) GO TO 50 +C + CALL XERCLR + CALL XERMSG ('SLATEC', 'XGMAINC', 'RESULT LT HALF PRECISION', 1, + + 1) +C +C 50 T = -A*ALX + LOG(ABS(H)) +C IF (T.LT.BOT) CALL XERCLR +C DGAMIT = SIGN (EXP(T), H) + 50 result = H + RETURN +C +C 60 T = T - A*ALX + 60 IF (T.LT.BOT) CALL XERCLR + result = -SGA * SGNGAM * EXP(T) + RETURN + + endif + return + end
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,9 @@ +2002-07-25 John W. Eaton <jwe@bevo.che.wisc.edu> + + * lo-specfun.cc (gammainc): New arg, err, for scalar version. + Use it in matrix versions to avoid spewing multiple errors. + Call xgammainc instead of dgamit. + 2002-07-22 John W. Eaton <jwe@bevo.che.wisc.edu> * CMatrix.cc (ComplexMatrix::ComplexMatrix (const boolMatrix&)):
--- a/liboctave/lo-specfun.cc +++ b/liboctave/lo-specfun.cc @@ -81,7 +81,7 @@ int F77_FUNC (xdgamma, XDGAMMA) (const double&, double&); - int F77_FUNC (xdgamit, XDGAMIT) (const double&, const double&, double&); + int F77_FUNC (xgammainc, XGAMMAINC) (const double&, const double&, double&); int F77_FUNC (dlgams, DLGAMS) (const double&, double&, double&); } @@ -841,16 +841,21 @@ // XXX FIXME XXX -- there is still room for improvement here... double -gammainc (double x, double a) +gammainc (double x, double a, bool& err) { double retval; - F77_XFCN (xdgamit, XDGAMIT, (a, x, retval)); + err = false; - if (x == 0.0) - retval = 0.0; - else if (x > 0.0) - retval = exp (a * log (x) + log (retval)); + if (a < 0.0 || x < 0.0) + { + (*current_liboctave_error_handler) + ("gammainc: A and X must be non-negative"); + + err = true; + } + else + F77_XFCN (xgammainc, XGAMMAINC, (a, x, retval)); return retval; } @@ -861,11 +866,23 @@ int nr = a.rows (); int nc = a.cols (); - Matrix retval (nr, nc); + Matrix result (nr, nc); + Matrix retval; + + bool err; for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) - retval(i,j) = gammainc (x, a(i,j)); + { + result(i,j) = gammainc (x, a(i,j), err); + + if (err) + goto done; + } + + retval = result; + + done: return retval; } @@ -876,11 +893,23 @@ int nr = x.rows (); int nc = x.cols (); - Matrix retval (nr, nc); + Matrix result (nr, nc); + Matrix retval; + + bool err; for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) - retval(i,j) = gammainc (x(i,j), a); + { + result(i,j) = gammainc (x(i,j), a, err); + + if (err) + goto done; + } + + retval = result; + + done: return retval; } @@ -888,6 +917,7 @@ Matrix gammainc (const Matrix& x, const Matrix& a) { + Matrix result; Matrix retval; int nr = x.rows (); @@ -898,17 +928,28 @@ if (nr == a_nr && nc == a_nc) { - retval.resize (nr, nc); + result.resize (nr, nc); + + bool err; for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) - retval(i,j) = gammainc (x(i,j), a(i,j)); + { + result(i,j) = gammainc (x(i,j), a(i,j), err); + + if (err) + goto done; + } + + retval = result; } else (*current_liboctave_error_handler) ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)", nr, nc, a_nr, a_nc); + done: + return retval; }
--- a/liboctave/lo-specfun.h +++ b/liboctave/lo-specfun.h @@ -188,11 +188,17 @@ extern Matrix betainc (const Matrix& x, const Matrix& a, double b); extern Matrix betainc (const Matrix& x, const Matrix& a, const Matrix& b); -extern double gammainc (double x, double a); +extern double gammainc (double x, double a, bool& err); extern Matrix gammainc (double x, const Matrix& a); extern Matrix gammainc (const Matrix& x, double a); extern Matrix gammainc (const Matrix& x, const Matrix& a); +inline double gammainc (double x, double a) +{ + bool err; + return gammainc (x, a, err); +} + #endif /*