Mercurial > hg > octave-max
changeset 3111:fe6f9bd9d0e6
[project @ 1997-11-26 07:52:06 by jwe]
| author | jwe |
|---|---|
| date | Wed, 26 Nov 1997 07:53:04 +0000 |
| parents | fe2d1ae8926b |
| children | 92394a9e4784 |
| files | ChangeLog Makeconf.in configure.in libcruft/ChangeLog libcruft/slatec-fn/Makefile libcruft/slatec-fn/d9gmit.f libcruft/slatec-fn/d9lgic.f libcruft/slatec-fn/d9lgit.f libcruft/slatec-fn/dbetai.f libcruft/slatec-fn/dgami.f libcruft/slatec-fn/dgamit.f libcruft/slatec-fn/dgamr.f libcruft/slatec-fn/dlbeta.f libcruft/slatec-fn/dlnrel.f src/ChangeLog src/Makefile.in src/version.h |
| diffstat | 17 files changed, 751 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,9 @@ +Wed Nov 26 00:38:31 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * configure.in (SPECIAL_MATH_LIB): If libdxml exists on DU + systems, define SPECIAL_MATH_LIB. + * Makeconf.in (SPECIAL_MATH_LIB): Substitute it. + Wed Nov 19 01:54:11 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * aclocal.m4 (OCTAVE_CXX_NEW_FRIEND_TEMPLATE_DECL): Don't forget
--- a/Makeconf.in +++ b/Makeconf.in @@ -140,6 +140,7 @@ LIBDLFCN = @LIBDLFCN@ LIBPLPLOT = @LIBPLPLOT@ LIBREADLINE = @LIBREADLINE@ +SPECIAL_MATH_LIB = @SPECIAL_MATH_LIB@ # The arguments passed to configure. config_opts = @config_opts@
--- a/configure.in +++ b/configure.in @@ -21,7 +21,7 @@ ### Software Foundation, 59 Temple Place - Suite 330, Boston, MA ### 02111-1307, USA. -AC_REVISION($Revision: 1.288 $) +AC_REVISION($Revision: 1.289 $) AC_PREREQ(2.9) AC_INIT(src/octave.cc) AC_CONFIG_HEADER(config.h) @@ -683,6 +683,13 @@ AC_CHECK_LIB(sun, getpwnam) AC_CHECK_LIB(socket, gethostname) +case "$canonical_host_type" in + alpha-dec-osf*) + AC_CHECK_LIB(dxml, dgemm_, [SPECIAL_MATH_LIB=-ldxml]) + ;; +esac +AC_SUBST(SPECIAL_MATH_LIB) + ### How big are ints and how are they oriented? These could probably ### be eliminated in favor of run-time checks.
--- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,3 +1,10 @@ +Wed Nov 26 01:49:47 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * slatec-fn/d9gmit.f, slatec-fn/d9lgic.f, slatec-fn/d9lgit.f, + slatec-fn/dbetai.f, slatec-fn/dgami.f, slatec-fn/dgamit.f, + slatec-fn/dgamr.f, slatec-fn/dlbeta.f, slatec-fn/dlnrel.f: + New files for incomplete beta and incomplete gamma functions. + Thu Jul 17 13:18:57 1997 Klaus Gebhardt <gebhardt@crunch.ikp.physik.th-darmstadt.de> * blas/xerbla.f (xerbla): Call XSTOPX instead of using STOP.
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/Makefile @@ -0,0 +1,19 @@ +# Generated automatically from Makefile.in by configure. +# +# Makefile for octave's libcruft/slatec-fn directory +# +# John W. Eaton +# jwe@bevo.che.wisc.edu +# University of Wisconsin-Madison +# Department of Chemical Engineering + +TOPDIR = ../.. + +srcdir = . +top_srcdir = ../.. + +EXTERNAL_DISTFILES = $(DISTFILES) + +include $(TOPDIR)/Makeconf + +include ../Makerules
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/d9gmit.f @@ -0,0 +1,91 @@ +*DECK D9GMIT + DOUBLE PRECISION FUNCTION D9GMIT (A, X, ALGAP1, SGNGAM, ALX) +C***BEGIN PROLOGUE D9GMIT +C***SUBSIDIARY +C***PURPOSE Compute Tricomi's incomplete Gamma function for small +C arguments. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C7E +C***TYPE DOUBLE PRECISION (R9GMIT-S, D9GMIT-D) +C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, SMALL X, +C SPECIAL FUNCTIONS, TRICOMI +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C Compute Tricomi's incomplete gamma function for small X. +C +C***REFERENCES (NONE) +C***ROUTINES CALLED D1MACH, DLNGAM, XERMSG +C***REVISION HISTORY (YYMMDD) +C 770701 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890911 Removed unnecessary intrinsics. (WRB) +C 890911 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) +C 900720 Routine changed from user-callable to subsidiary. (WRB) +C***END PROLOGUE D9GMIT + DOUBLE PRECISION A, X, ALGAP1, SGNGAM, ALX, AE, AEPS, ALGS, ALG2, + 1 BOT, EPS, FK, S, SGNG2, T, TE, D1MACH, DLNGAM + LOGICAL FIRST + SAVE EPS, BOT, FIRST + DATA FIRST /.TRUE./ +C***FIRST EXECUTABLE STATEMENT D9GMIT + IF (FIRST) THEN + EPS = 0.5D0*D1MACH(3) + BOT = LOG (D1MACH(1)) + ENDIF + FIRST = .FALSE. +C + IF (X .LE. 0.D0) CALL XERMSG ('SLATEC', 'D9GMIT', + + 'X SHOULD BE GT 0', 1, 2) +C + MA = A + 0.5D0 + IF (A.LT.0.D0) MA = A - 0.5D0 + AEPS = A - MA +C + AE = A + IF (A.LT.(-0.5D0)) AE = AEPS +C + T = 1.D0 + TE = AE + S = T + DO 20 K=1,200 + FK = K + TE = -X*TE/FK + T = TE/(AE+FK) + S = S + T + IF (ABS(T).LT.EPS*ABS(S)) GO TO 30 + 20 CONTINUE + CALL XERMSG ('SLATEC', 'D9GMIT', + + 'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2) +C + 30 IF (A.GE.(-0.5D0)) ALGS = -ALGAP1 + LOG(S) + IF (A.GE.(-0.5D0)) GO TO 60 +C + ALGS = -DLNGAM(1.D0+AEPS) + LOG(S) + S = 1.0D0 + M = -MA - 1 + IF (M.EQ.0) GO TO 50 + T = 1.0D0 + DO 40 K=1,M + T = X*T/(AEPS-(M+1-K)) + S = S + T + IF (ABS(T).LT.EPS*ABS(S)) GO TO 50 + 40 CONTINUE +C + 50 D9GMIT = 0.0D0 + ALGS = -MA*LOG(X) + ALGS + IF (S.EQ.0.D0 .OR. AEPS.EQ.0.D0) GO TO 60 +C + SGNG2 = SGNGAM * SIGN (1.0D0, S) + ALG2 = -X - ALGAP1 + LOG(ABS(S)) +C + IF (ALG2.GT.BOT) D9GMIT = SGNG2 * EXP(ALG2) + IF (ALGS.GT.BOT) D9GMIT = D9GMIT + EXP(ALGS) + RETURN +C + 60 D9GMIT = EXP (ALGS) + RETURN +C + END
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/d9lgic.f @@ -0,0 +1,54 @@ +*DECK D9LGIC + DOUBLE PRECISION FUNCTION D9LGIC (A, X, ALX) +C***BEGIN PROLOGUE D9LGIC +C***SUBSIDIARY +C***PURPOSE Compute the log complementary incomplete Gamma function +C for large X and for A .LE. X. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C7E +C***TYPE DOUBLE PRECISION (R9LGIC-S, D9LGIC-D) +C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, LARGE X, +C LOGARITHM, SPECIAL FUNCTIONS +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C Compute the log complementary incomplete gamma function for large X +C and for A .LE. X. +C +C***REFERENCES (NONE) +C***ROUTINES CALLED D1MACH, XERMSG +C***REVISION HISTORY (YYMMDD) +C 770701 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890531 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) +C 900720 Routine changed from user-callable to subsidiary. (WRB) +C***END PROLOGUE D9LGIC + DOUBLE PRECISION A, X, ALX, EPS, FK, P, R, S, T, XMA, XPA, D1MACH + SAVE EPS + DATA EPS / 0.D0 / +C***FIRST EXECUTABLE STATEMENT D9LGIC + IF (EPS.EQ.0.D0) EPS = 0.5D0*D1MACH(3) +C + XPA = X + 1.0D0 - A + XMA = X - 1.D0 - A +C + R = 0.D0 + P = 1.D0 + S = P + DO 10 K=1,300 + FK = K + T = FK*(A-FK)*(1.D0+R) + R = -T/((XMA+2.D0*FK)*(XPA+2.D0*FK)+T) + P = R*P + S = S + P + IF (ABS(P).LT.EPS*S) GO TO 20 + 10 CONTINUE + CALL XERMSG ('SLATEC', 'D9LGIC', + + 'NO CONVERGENCE IN 300 TERMS OF CONTINUED FRACTION', 1, 2) +C + 20 D9LGIC = A*ALX - X + LOG(S/XPA) +C + RETURN + END
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/d9lgit.f @@ -0,0 +1,67 @@ +*DECK D9LGIT + DOUBLE PRECISION FUNCTION D9LGIT (A, X, ALGAP1) +C***BEGIN PROLOGUE D9LGIT +C***SUBSIDIARY +C***PURPOSE Compute the logarithm of Tricomi's incomplete Gamma +C function with Perron's continued fraction for large X and +C A .GE. X. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C7E +C***TYPE DOUBLE PRECISION (R9LGIT-S, D9LGIT-D) +C***KEYWORDS FNLIB, INCOMPLETE GAMMA FUNCTION, LOGARITHM, +C PERRON'S CONTINUED FRACTION, SPECIAL FUNCTIONS, TRICOMI +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C Compute the log of Tricomi's incomplete gamma function with Perron's +C continued fraction for large X and for A .GE. X. +C +C***REFERENCES (NONE) +C***ROUTINES CALLED D1MACH, XERMSG +C***REVISION HISTORY (YYMMDD) +C 770701 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890531 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) +C 900720 Routine changed from user-callable to subsidiary. (WRB) +C***END PROLOGUE D9LGIT + DOUBLE PRECISION A, X, ALGAP1, AX, A1X, EPS, FK, HSTAR, P, R, S, + 1 SQEPS, T, D1MACH + LOGICAL FIRST + SAVE EPS, SQEPS, FIRST + DATA FIRST /.TRUE./ +C***FIRST EXECUTABLE STATEMENT D9LGIT + IF (FIRST) THEN + EPS = 0.5D0*D1MACH(3) + SQEPS = SQRT(D1MACH(4)) + ENDIF + FIRST = .FALSE. +C + IF (X .LE. 0.D0 .OR. A .LT. X) CALL XERMSG ('SLATEC', 'D9LGIT', + + 'X SHOULD BE GT 0.0 AND LE A', 2, 2) +C + AX = A + X + A1X = AX + 1.0D0 + R = 0.D0 + P = 1.D0 + S = P + DO 20 K=1,200 + FK = K + T = (A+FK)*X*(1.D0+R) + R = T/((AX+FK)*(A1X+FK)-T) + P = R*P + S = S + P + IF (ABS(P).LT.EPS*S) GO TO 30 + 20 CONTINUE + CALL XERMSG ('SLATEC', 'D9LGIT', + + 'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2) +C + 30 HSTAR = 1.0D0 - X*S/A1X + IF (HSTAR .LT. SQEPS) CALL XERMSG ('SLATEC', 'D9LGIT', + + 'RESULT LESS THAN HALF PRECISION', 1, 1) +C + D9LGIT = -X - ALGAP1 - LOG(HSTAR) + RETURN +C + END
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/dbetai.f @@ -0,0 +1,121 @@ + +*DECK DBETAI + DOUBLE PRECISION FUNCTION DBETAI (X, PIN, QIN) +C***BEGIN PROLOGUE DBETAI +C***PURPOSE Calculate the incomplete Beta function. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C7F +C***TYPE DOUBLE PRECISION (BETAI-S, DBETAI-D) +C***KEYWORDS FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C DBETAI calculates the DOUBLE PRECISION incomplete beta function. +C +C The incomplete beta function ratio is the probability that a +C random variable from a beta distribution having parameters PIN and +C QIN will be less than or equal to X. +C +C -- Input Arguments -- All arguments are DOUBLE PRECISION. +C X upper limit of integration. X must be in (0,1) inclusive. +C PIN first beta distribution parameter. PIN must be .GT. 0.0. +C QIN second beta distribution parameter. QIN must be .GT. 0.0. +C +C***REFERENCES Nancy E. Bosten and E. L. Battiste, Remark on Algorithm +C 179, Communications of the ACM 17, 3 (March 1974), +C pp. 156. +C***ROUTINES CALLED D1MACH, DLBETA, XERMSG +C***REVISION HISTORY (YYMMDD) +C 770701 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890911 Removed unnecessary intrinsics. (WRB) +C 890911 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) +C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) +C***END PROLOGUE DBETAI + DOUBLE PRECISION X, PIN, QIN, ALNEPS, ALNSML, C, EPS, FINSUM, P, + 1 PS, Q, SML, TERM, XB, XI, Y, D1MACH, DLBETA, P1 + LOGICAL FIRST + SAVE EPS, ALNEPS, SML, ALNSML, FIRST + DATA FIRST /.TRUE./ +C***FIRST EXECUTABLE STATEMENT DBETAI + IF (FIRST) THEN + EPS = D1MACH(3) + ALNEPS = LOG (EPS) + SML = D1MACH(1) + ALNSML = LOG (SML) + ENDIF + FIRST = .FALSE. +C + IF (X .LT. 0.D0 .OR. X .GT. 1.D0) CALL XERMSG ('SLATEC', 'DBETAI', + + 'X IS NOT IN THE RANGE (0,1)', 1, 2) + IF (PIN .LE. 0.D0 .OR. QIN .LE. 0.D0) CALL XERMSG ('SLATEC', + + 'DBETAI', 'P AND/OR Q IS LE ZERO', 2, 2) +C + Y = X + P = PIN + Q = QIN + IF (Q.LE.P .AND. X.LT.0.8D0) GO TO 20 + IF (X.LT.0.2D0) GO TO 20 + Y = 1.0D0 - Y + P = QIN + Q = PIN +C + 20 IF ((P+Q)*Y/(P+1.D0).LT.EPS) GO TO 80 +C +C EVALUATE THE INFINITE SUM FIRST. TERM WILL EQUAL +C Y**P/BETA(PS,P) * (1.-PS)-SUB-I * Y**I / FAC(I) . +C + PS = Q - AINT(Q) + IF (PS.EQ.0.D0) PS = 1.0D0 + XB = P*LOG(Y) - DLBETA(PS,P) - LOG(P) + DBETAI = 0.0D0 + IF (XB.LT.ALNSML) GO TO 40 +C + DBETAI = EXP (XB) + TERM = DBETAI*P + IF (PS.EQ.1.0D0) GO TO 40 + N = MAX (ALNEPS/LOG(Y), 4.0D0) + DO 30 I=1,N + XI = I + TERM = TERM * (XI-PS)*Y/XI + DBETAI = DBETAI + TERM/(P+XI) + 30 CONTINUE +C +C NOW EVALUATE THE FINITE SUM, MAYBE. +C + 40 IF (Q.LE.1.0D0) GO TO 70 +C + XB = P*LOG(Y) + Q*LOG(1.0D0-Y) - DLBETA(P,Q) - LOG(Q) + IB = MAX (XB/ALNSML, 0.0D0) + TERM = EXP(XB - IB*ALNSML) + C = 1.0D0/(1.D0-Y) + P1 = Q*C/(P+Q-1.D0) +C + FINSUM = 0.0D0 + N = Q + IF (Q.EQ.DBLE(N)) N = N - 1 + DO 50 I=1,N + IF (P1.LE.1.0D0 .AND. TERM/EPS.LE.FINSUM) GO TO 60 + XI = I + TERM = (Q-XI+1.0D0)*C*TERM/(P+Q-XI) +C + IF (TERM.GT.1.0D0) IB = IB - 1 + IF (TERM.GT.1.0D0) TERM = TERM*SML +C + IF (IB.EQ.0) FINSUM = FINSUM + TERM + 50 CONTINUE +C + 60 DBETAI = DBETAI + FINSUM + 70 IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI + DBETAI = MAX (MIN (DBETAI, 1.0D0), 0.0D0) + RETURN +C + 80 DBETAI = 0.0D0 + XB = P*LOG(MAX(Y,SML)) - LOG(P) - DLBETA(P,Q) + IF (XB.GT.ALNSML .AND. Y.NE.0.0D0) DBETAI = EXP(XB) + IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI +C + RETURN + END
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/dgami.f @@ -0,0 +1,47 @@ + +*DECK DGAMI + DOUBLE PRECISION FUNCTION DGAMI (A, X) +C***BEGIN PROLOGUE DGAMI +C***PURPOSE Evaluate the incomplete Gamma function. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C7E +C***TYPE DOUBLE PRECISION (GAMI-S, DGAMI-D) +C***KEYWORDS FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C Evaluate the incomplete gamma function defined by +C +C DGAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) . +C +C DGAMI is evaluated for positive values of A and non-negative values +C of X. A slight deterioration of 2 or 3 digits accuracy will occur +C when DGAMI is very large or very small, because logarithmic variables +C are used. The function and both arguments are double precision. +C +C***REFERENCES (NONE) +C***ROUTINES CALLED DGAMIT, DLNGAM, XERMSG +C***REVISION HISTORY (YYMMDD) +C 770701 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890531 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) +C***END PROLOGUE DGAMI + DOUBLE PRECISION A, X, FACTOR, DLNGAM, DGAMIT +C***FIRST EXECUTABLE STATEMENT DGAMI + IF (A .LE. 0.D0) CALL XERMSG ('SLATEC', 'DGAMI', + + 'A MUST BE GT ZERO', 1, 2) + IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMI', + + 'X MUST BE GE ZERO', 2, 2) +C + DGAMI = 0.D0 + IF (X.EQ.0.0D0) RETURN +C +C THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW. + FACTOR = EXP (DLNGAM(A) + A*LOG(X)) +C + DGAMI = FACTOR * DGAMIT (A, X) +C + RETURN + END
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/dgamit.f @@ -0,0 +1,119 @@ +*DECK DGAMIT + DOUBLE PRECISION FUNCTION DGAMIT (A, X) +C***BEGIN PROLOGUE DGAMIT +C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C7E +C***TYPE DOUBLE PRECISION (GAMIT-S, DGAMIT-D) +C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, +C SPECIAL FUNCTIONS, TRICOMI +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C Evaluate Tricomi's incomplete Gamma function defined by +C +C DGAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) * +C T**(A-1.) +C +C for A .GT. 0.0 and by analytic continuation for A .LE. 0.0. +C GAMMA(X) is the complete gamma function of X. +C +C DGAMIT is evaluated for arbitrary real values of A and for non- +C negative values of X (even though DGAMIT is defined for X .LT. +C 0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite, +C which is a fatal error. +C +C The function and both arguments are DOUBLE PRECISION. +C +C A slight deterioration of 2 or 3 digits accuracy will occur when +C DGAMIT is very large or very small in absolute value, because log- +C arithmic variables are used. Also, if the parameter A is very +C close to a negative integer (but not a negative integer), there is +C a loss of accuracy, which is reported if the result is less than +C half machine precision. +C +C***REFERENCES W. Gautschi, A computational procedure for incomplete +C gamma functions, ACM Transactions on Mathematical +C Software 5, 4 (December 1979), pp. 466-481. +C W. Gautschi, Incomplete gamma functions, Algorithm 542, +C ACM Transactions on Mathematical Software 5, 4 +C (December 1979), pp. 482-489. +C***ROUTINES CALLED D1MACH, D9GMIT, D9LGIC, D9LGIT, DGAMR, DLGAMS, +C DLNGAM, XERCLR, XERMSG +C***REVISION HISTORY (YYMMDD) +C 770701 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890531 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) +C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) +C***END PROLOGUE DGAMIT + DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX, + 1 BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, DGAMR, D9GMIT, D9LGIT, + 2 DLNGAM, D9LGIC + LOGICAL FIRST + SAVE ALNEPS, SQEPS, BOT, FIRST + DATA FIRST /.TRUE./ +C***FIRST EXECUTABLE STATEMENT DGAMIT + 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', 'DGAMIT', '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 + IF (X.GT.0.D0) GO TO 20 + DGAMIT = 0.0D0 + IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0) + 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) + DGAMIT = 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 + DGAMIT = EXP (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', 'DGAMIT', 'RESULT LT HALF PRECISION', 1, + + 1) +C + 50 T = -A*ALX + LOG(ABS(H)) + IF (T.LT.BOT) CALL XERCLR + DGAMIT = SIGN (EXP(T), H) + RETURN +C + 60 T = T - A*ALX + IF (T.LT.BOT) CALL XERCLR + DGAMIT = -SGA * SGNGAM * EXP(T) + RETURN +C + END
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/dgamr.f @@ -0,0 +1,44 @@ +*DECK DGAMR + DOUBLE PRECISION FUNCTION DGAMR (X) +C***BEGIN PROLOGUE DGAMR +C***PURPOSE Compute the reciprocal of the Gamma function. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C7A +C***TYPE DOUBLE PRECISION (GAMR-S, DGAMR-D, CGAMR-C) +C***KEYWORDS FNLIB, RECIPROCAL GAMMA FUNCTION, SPECIAL FUNCTIONS +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C DGAMR(X) calculates the double precision reciprocal of the +C complete Gamma function for double precision argument X. +C +C***REFERENCES (NONE) +C***ROUTINES CALLED DGAMMA, DLGAMS, XERCLR, XGETF, XSETF +C***REVISION HISTORY (YYMMDD) +C 770701 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890531 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900727 Added EXTERNAL statement. (WRB) +C***END PROLOGUE DGAMR + DOUBLE PRECISION X, ALNGX, SGNGX, DGAMMA + EXTERNAL DGAMMA +C***FIRST EXECUTABLE STATEMENT DGAMR + DGAMR = 0.0D0 + IF (X.LE.0.0D0 .AND. AINT(X).EQ.X) RETURN +C + CALL XGETF (IROLD) + CALL XSETF (1) + IF (ABS(X).GT.10.0D0) GO TO 10 + DGAMR = 1.0D0/DGAMMA(X) + CALL XERCLR + CALL XSETF (IROLD) + RETURN +C + 10 CALL DLGAMS (X, ALNGX, SGNGX) + CALL XERCLR + CALL XSETF (IROLD) + DGAMR = SGNGX * EXP(-ALNGX) + RETURN +C + END
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/dlbeta.f @@ -0,0 +1,62 @@ +*DECK DLBETA + DOUBLE PRECISION FUNCTION DLBETA (A, B) +C***BEGIN PROLOGUE DLBETA +C***PURPOSE Compute the natural logarithm of the complete Beta +C function. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C7B +C***TYPE DOUBLE PRECISION (ALBETA-S, DLBETA-D, CLBETA-C) +C***KEYWORDS FNLIB, LOGARITHM OF THE COMPLETE BETA FUNCTION, +C SPECIAL FUNCTIONS +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C DLBETA(A,B) calculates the double precision natural logarithm of +C the complete beta function for double precision arguments +C A and B. +C +C***REFERENCES (NONE) +C***ROUTINES CALLED D9LGMC, DGAMMA, DLNGAM, DLNREL, XERMSG +C***REVISION HISTORY (YYMMDD) +C 770701 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890531 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) +C 900727 Added EXTERNAL statement. (WRB) +C***END PROLOGUE DLBETA + DOUBLE PRECISION A, B, P, Q, CORR, SQ2PIL, D9LGMC, DGAMMA, DLNGAM, + 1 DLNREL + EXTERNAL DGAMMA + SAVE SQ2PIL + DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 / +C***FIRST EXECUTABLE STATEMENT DLBETA + P = MIN (A, B) + Q = MAX (A, B) +C + IF (P .LE. 0.D0) CALL XERMSG ('SLATEC', 'DLBETA', + + 'BOTH ARGUMENTS MUST BE GT ZERO', 1, 2) +C + IF (P.GE.10.D0) GO TO 30 + IF (Q.GE.10.D0) GO TO 20 +C +C P AND Q ARE SMALL. +C + DLBETA = LOG (DGAMMA(P) * (DGAMMA(Q)/DGAMMA(P+Q)) ) + RETURN +C +C P IS SMALL, BUT Q IS BIG. +C + 20 CORR = D9LGMC(Q) - D9LGMC(P+Q) + DLBETA = DLNGAM(P) + CORR + P - P*LOG(P+Q) + 1 + (Q-0.5D0)*DLNREL(-P/(P+Q)) + RETURN +C +C P AND Q ARE BIG. +C + 30 CORR = D9LGMC(P) + D9LGMC(Q) - D9LGMC(P+Q) + DLBETA = -0.5D0*LOG(Q) + SQ2PIL + CORR + (P-0.5D0)*LOG(P/(P+Q)) + 1 + Q*DLNREL(-P/(P+Q)) + RETURN +C + END
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/dlnrel.f @@ -0,0 +1,98 @@ +*DECK DLNREL + DOUBLE PRECISION FUNCTION DLNREL (X) +C***BEGIN PROLOGUE DLNREL +C***PURPOSE Evaluate ln(1+X) accurate in the sense of relative error. +C***LIBRARY SLATEC (FNLIB) +C***CATEGORY C4B +C***TYPE DOUBLE PRECISION (ALNREL-S, DLNREL-D, CLNREL-C) +C***KEYWORDS ELEMENTARY FUNCTIONS, FNLIB, LOGARITHM +C***AUTHOR Fullerton, W., (LANL) +C***DESCRIPTION +C +C DLNREL(X) calculates the double precision natural logarithm of +C (1.0+X) for double precision argument X. This routine should +C be used when X is small and accurate to calculate the logarithm +C accurately (in the relative error sense) in the neighborhood +C of 1.0. +C +C Series for ALNR on the interval -3.75000E-01 to 3.75000E-01 +C with weighted error 6.35E-32 +C log weighted error 31.20 +C significant figures required 30.93 +C decimal places required 32.01 +C +C***REFERENCES (NONE) +C***ROUTINES CALLED D1MACH, DCSEVL, INITDS, XERMSG +C***REVISION HISTORY (YYMMDD) +C 770601 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890531 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) +C***END PROLOGUE DLNREL + DOUBLE PRECISION ALNRCS(43), X, XMIN, DCSEVL, D1MACH + LOGICAL FIRST + SAVE ALNRCS, NLNREL, XMIN, FIRST + DATA ALNRCS( 1) / +.1037869356 2743769800 6862677190 98 D+1 / + DATA ALNRCS( 2) / -.1336430150 4908918098 7660415531 33 D+0 / + DATA ALNRCS( 3) / +.1940824913 5520563357 9261993747 50 D-1 / + DATA ALNRCS( 4) / -.3010755112 7535777690 3765377765 92 D-2 / + DATA ALNRCS( 5) / +.4869461479 7154850090 4563665091 37 D-3 / + DATA ALNRCS( 6) / -.8105488189 3175356066 8099430086 22 D-4 / + DATA ALNRCS( 7) / +.1377884779 9559524782 9382514960 59 D-4 / + DATA ALNRCS( 8) / -.2380221089 4358970251 3699929149 35 D-5 / + DATA ALNRCS( 9) / +.4164041621 3865183476 3918599019 89 D-6 / + DATA ALNRCS( 10) / -.7359582837 8075994984 2668370319 98 D-7 / + DATA ALNRCS( 11) / +.1311761187 6241674949 1522943450 11 D-7 / + DATA ALNRCS( 12) / -.2354670931 7742425136 6960923301 75 D-8 / + DATA ALNRCS( 13) / +.4252277327 6034997775 6380529625 67 D-9 / + DATA ALNRCS( 14) / -.7719089413 4840796826 1081074933 00 D-10 / + DATA ALNRCS( 15) / +.1407574648 1359069909 2153564721 91 D-10 / + DATA ALNRCS( 16) / -.2576907205 8024680627 5370786275 84 D-11 / + DATA ALNRCS( 17) / +.4734240666 6294421849 1543950059 38 D-12 / + DATA ALNRCS( 18) / -.8724901267 4742641745 3012632926 75 D-13 / + DATA ALNRCS( 19) / +.1612461490 2740551465 7398331191 15 D-13 / + DATA ALNRCS( 20) / -.2987565201 5665773006 7107924168 15 D-14 / + DATA ALNRCS( 21) / +.5548070120 9082887983 0413216972 79 D-15 / + DATA ALNRCS( 22) / -.1032461915 8271569595 1413339619 32 D-15 / + DATA ALNRCS( 23) / +.1925023920 3049851177 8785032448 68 D-16 / + DATA ALNRCS( 24) / -.3595507346 5265150011 1897078442 66 D-17 / + DATA ALNRCS( 25) / +.6726454253 7876857892 1945742267 73 D-18 / + DATA ALNRCS( 26) / -.1260262416 8735219252 0824256375 46 D-18 / + DATA ALNRCS( 27) / +.2364488440 8606210044 9161589555 19 D-19 / + DATA ALNRCS( 28) / -.4441937705 0807936898 8783891797 33 D-20 / + DATA ALNRCS( 29) / +.8354659446 4034259016 2412939946 66 D-21 / + DATA ALNRCS( 30) / -.1573155941 6479562574 8992535210 66 D-21 / + DATA ALNRCS( 31) / +.2965312874 0247422686 1543697066 66 D-22 / + DATA ALNRCS( 32) / -.5594958348 1815947292 1560132266 66 D-23 / + DATA ALNRCS( 33) / +.1056635426 8835681048 1872841386 66 D-23 / + DATA ALNRCS( 34) / -.1997248368 0670204548 3149994666 66 D-24 / + DATA ALNRCS( 35) / +.3778297781 8839361421 0498559999 99 D-25 / + DATA ALNRCS( 36) / -.7153158688 9081740345 0381653333 33 D-26 / + DATA ALNRCS( 37) / +.1355248846 3674213646 5020245333 33 D-26 / + DATA ALNRCS( 38) / -.2569467304 8487567430 0798293333 33 D-27 / + DATA ALNRCS( 39) / +.4874775606 6216949076 4595199999 99 D-28 / + DATA ALNRCS( 40) / -.9254211253 0849715321 1323733333 33 D-29 / + DATA ALNRCS( 41) / +.1757859784 1760239233 2697600000 00 D-29 / + DATA ALNRCS( 42) / -.3341002667 7731010351 3770666666 66 D-30 / + DATA ALNRCS( 43) / +.6353393618 0236187354 1802666666 66 D-31 / + DATA FIRST /.TRUE./ +C***FIRST EXECUTABLE STATEMENT DLNREL + IF (FIRST) THEN + NLNREL = INITDS (ALNRCS, 43, 0.1*REAL(D1MACH(3))) + XMIN = -1.0D0 + SQRT(D1MACH(4)) + ENDIF + FIRST = .FALSE. +C + IF (X .LE. (-1.D0)) CALL XERMSG ('SLATEC', 'DLNREL', 'X IS LE -1' + + , 2, 2) + IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'DLNREL', + + 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR -1', 1, 1) +C + IF (ABS(X).LE.0.375D0) DLNREL = X*(1.D0 - + 1 X*DCSEVL (X/.375D0, ALNRCS, NLNREL)) +C + IF (ABS(X).GT.0.375D0) DLNREL = LOG (1.0D0+X) +C + RETURN + END
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +Wed Nov 26 00:39:34 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * Makefile.in (OCTAVE_LIBS): Include $(SPECIAL_MATH_LIB) just + ahead of -lcruft. + Thu Nov 20 15:16:22 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * octave.cc (maximum_braindamage): Bind implicit_num_to_str_ok to 1.
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -171,7 +171,7 @@ OCTAVE_LFLAGS = -L../liboctave -L../libcruft -L../glob \ -L../dlfcn -L. $(RLD_FLAG) -OCTAVE_LIBS = -loctinterp -loctave -lcruft $(LIBPLPLOT) \ +OCTAVE_LIBS = -loctinterp -loctave $(SPECIAL_MATH_LIB) -lcruft $(LIBPLPLOT) \ $(LIBREADLINE) ../kpathsea/libkpathsea.$(LIBEXT) -lglob $(LIBDLFCN) LIBS = @LIBS@