Mercurial > octave-nkf

changeset 4004:ca854fb51a88

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[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	Thu Jul 25 01:19:17 2002 +0000
+++ b/libcruft/ChangeLog	Thu Jul 25 06:32:37 2002 +0000
@@ -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.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/slatec-fn/xgmainc.f	Thu Jul 25 06:32:37 2002 +0000
@@ -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	Thu Jul 25 01:19:17 2002 +0000
+++ b/liboctave/ChangeLog	Thu Jul 25 06:32:37 2002 +0000
@@ -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	Thu Jul 25 01:19:17 2002 +0000
+++ b/liboctave/lo-specfun.cc	Thu Jul 25 06:32:37 2002 +0000
@@ -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	Thu Jul 25 01:19:17 2002 +0000
+++ b/liboctave/lo-specfun.h	Thu Jul 25 06:32:37 2002 +0000
@@ -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
 
 /*