Mercurial > hg > octave-nkf

changeset 3117:f735f3ea1ee7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1997-11-29 19:43:23 by jwe]
author jwe
date Sat, 29 Nov 1997 19:43:24 +0000
parents 8f31cbbcf555
children 74cc8e2fe2c0
files libcruft/specfun/ribesl.f libcruft/specfun/rjbesl.f libcruft/specfun/rkbesl.f libcruft/specfun/rybesl.f
diffstat 4 files changed, 79 insertions(+), 32 deletions(-) [+]
line wrap: on
line diff
--- a/libcruft/specfun/ribesl.f
+++ b/libcruft/specfun/ribesl.f
@@ -167,8 +167,9 @@
 C               Argonne, IL  60439
 C
 C-------------------------------------------------------------------
+      LOGICAL FIRST
       INTEGER IZE,K,L,MAGX,N,NB,NBMX,NCALC,NEND,NSIG,NSTART
-      DOUBLE PRECISION DGAMMA,             
+      DOUBLE PRECISION DGAMMA, D1MACH
      1 ALPHA,B,CONST,CONV,EM,EMPAL,EMP2AL,EN,ENMTEN,ENSIG,
      2 ENTEN,EXPARG,FUNC,HALF,HALFX,ONE,P,PLAST,POLD,PSAVE,PSAVEL,
      3 RTNSIG,SUM,TEMPA,TEMPB,TEMPC,TEST,TOVER,TWO,X,XLARGE,ZERO
@@ -176,13 +177,25 @@
 C-------------------------------------------------------------------
 C  Mathematical constants
 C-------------------------------------------------------------------
-      DATA ONE,TWO,ZERO,HALF,CONST/1.0D0,2.0D0,0.0D0,0.5D0,1.585D0/
+      PARAMETER (ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0)
+      PARAMETER (HALF = 0.5D0, CONST = 1.585D0)
 C-------------------------------------------------------------------
 C  Machine-dependent parameters
 C-------------------------------------------------------------------
       DATA NSIG,XLARGE,EXPARG /16,1.0D4,709.0D0/
       DATA ENTEN,ENSIG,RTNSIG/1.0D308,1.0D16,1.0D-4/
       DATA ENMTEN/8.9D-308/
+      DATA FIRST /.TRUE./
+      IF (FIRST) THEN
+        NSIG = NINT (-LOG (D1MACH (1)))
+        ENTEN = 1.0D1 ** (INT (LOG10 (D1MACH (2)))
+        ENSIG = 1.0D1 ** NSIG
+        RTNSIG = 1.0D1 ** (-NINT (NSIG / 4.0))
+        ENMTEN = 4.0D0 * D1MACH (1)
+        EXPARG = LOG (D1MACH (2))
+        XLARGE = 1.0D4
+        FIRST = .FALSE.
+      ENDIF
 C-------------------------------------------------------------------
 C  Statement functions for conversion
 C-------------------------------------------------------------------
--- a/libcruft/specfun/rjbesl.f
+++ b/libcruft/specfun/rjbesl.f
@@ -137,8 +137,9 @@
 C          Argonne, IL  60439
 C
 C---------------------------------------------------------------------
-      INTEGER I,J,K,L,M,MAGX,N,NB,NBMX,NCALC,NEND,NSTART
-      DOUBLE PRECISION  DGAMMA,
+      LOGICAL FIRST
+      INTEGER I,J,K,L,M,MAGX,N,NB,NBMX,NCALC,NEND,NSIG,NSTART
+      DOUBLE PRECISION  DGAMMA, D1MACH
      1 ALPHA,ALPEM,ALP2EM,B,CAPP,CAPQ,CONV,EIGHTH,EM,EN,ENMTEN,ENSIG,
      2 ENTEN,FACT,FOUR,FUNC,GNU,HALF,HALFX,ONE,ONE30,P,PI2,PLAST,
      3 POLD,PSAVE,PSAVEL,RTNSIG,S,SUM,T,T1,TEMPA,TEMPB,TEMPC,TEST,
@@ -153,16 +154,13 @@
 C   TWOPI2 - (2*PI - TWOPI) to working precision, i.e.,
 C            TWOPI1 + TWOPI2 = 2 * PI to extra precision.
 C---------------------------------------------------------------------
-      DATA PI2, TWOPI1, TWOPI2 /0.636619772367581343075535D0,6.28125D0,
-     1 1.935307179586476925286767D-3/
-      DATA ZERO, EIGHTH, HALF, ONE /0.0D0,0.125D0,0.5D0,1.0D0/
-      DATA TWO, THREE, FOUR, TWOFIV /2.0D0,3.0D0,4.0D0,25.0D0/
-      DATA ONE30, THREE5 /130.0D0,35.0D0/
-C---------------------------------------------------------------------
-C  Machine-dependent parameters
-C---------------------------------------------------------------------
-      DATA ENTEN, ENSIG, RTNSIG /1.0D38,1.0D17,1.0D-4/
-      DATA ENMTEN, XLARGE /1.2D-37,1.0D4/
+      PARAMETER (PI2 = 0.636619772367581343075535D0)
+      PARAMETER (TWOPI1 = 6.28125D0,
+      PARAMETER (TWOPI2 = 1.935307179586476925286767D-3)
+      PARAMETER (ZERO = 0.0D0, EIGHTH = 0.125D0, HALF = 0.5D0)
+      PARAMETER (ONE = 1.0D0, TWO = 2.0D0, THREE = 3.0D0)
+      PARAMETER (FOUR = 4.0D0, TWOFIV = 2.5D1, ONE30 = 1.3D2)
+      PARAMETER (THREE5 = 3.5D1)
 C---------------------------------------------------------------------
 C     Factorial(N)
 C---------------------------------------------------------------------
@@ -173,6 +171,19 @@
      4 5.109094217170944D19,1.12400072777760768D21,
      5 2.585201673888497664D22,6.2044840173323943936D23/
 C---------------------------------------------------------------------
+C  Machine-dependent parameters
+C---------------------------------------------------------------------
+      DATA FIRST /.TRUE./
+      IF (FIRST) THEN
+        NSIG = NINT (-LOG (D1MACH (1)))
+        ENTEN = 1.0D1 ** (INT (LOG10 (D1MACH (2)))
+        ENSIG = 1.0D1 ** NSIG
+        RTNSIG = 1.0D1 ** (-NINT (NSIG / 4.0))
+        ENMTEN = 4.0D0 * D1MACH (1)
+        XLARGE = 1.0D4
+        FIRST = .FALSE.
+      ENDIF
+C---------------------------------------------------------------------
 C Statement functions for conversion and the gamma function.
 C---------------------------------------------------------------------
       CONV(I) = DBLE(I)
--- a/libcruft/specfun/rkbesl.f
+++ b/libcruft/specfun/rkbesl.f
@@ -141,7 +141,9 @@
 C               Argonne, IL  60439
 C
 C-------------------------------------------------------------------
+      LOGICAL FIRST
       INTEGER I,IEND,ITEMP,IZE,J,K,M,MPLUS1,NB,NCALC
+      DOUBLE PRECISION D1MACH
       DOUBLE PRECISION  
      1    A,ALPHA,BLPHA,BK,BK1,BK2,C,D,DM,D1,D2,D3,ENU,EPS,ESTF,ESTM,
      2    EX,FOUR,F0,F1,F2,HALF,ONE,P,P0,Q,Q0,R,RATIO,S,SQXMIN,T,TINYX,
@@ -152,14 +154,10 @@
 C    A = LOG(2.D0) - Euler's constant
 C    D = SQRT(2.D0/PI)
 C---------------------------------------------------------------------
-      DATA HALF,ONE,TWO,ZERO/0.5D0,1.0D0,2.0D0,0.0D0/
-      DATA FOUR,TINYX/4.0D0,1.0D-10/
-      DATA A/ 0.11593151565841244881D0/,D/0.797884560802865364D0/
-C---------------------------------------------------------------------
-C  Machine dependent parameters
-C---------------------------------------------------------------------
-      DATA EPS/2.22D-16/,SQXMIN/1.49D-154/,XINF/1.79D+308/
-      DATA XMIN/2.23D-308/,XMAX/705.342D0/
+      PARAMETER (HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0)
+      PARAMETER (FOUR = 4.0D0, TINYX = 1.0D-10)
+      PARAMETER (A = 0.11593151565841244881D0)
+      PARAMETER (D = 0.797884560802865364D0)
 C---------------------------------------------------------------------
 C  P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA
 C                                         + Euler's constant
@@ -189,6 +187,18 @@
       DATA ESTF/4.18341D1, 7.1075D0, 6.4306D0, 4.25110D1, 1.35633D0,
      1          8.45096D1, 2.0D1/
 C---------------------------------------------------------------------
+C  Machine dependent parameters
+C---------------------------------------------------------------------
+      DATA FIRST /.TRUE./
+      IF (FIRST) THEN
+        EPS = D1MACH (4)
+        XINF = D1MACH (2)
+        XMIN = D1MACH (1)
+        SQXMIN = SQRT (XMIN)
+        XMAX = 0.99D0 * LOG (XINF)
+        FIRST = .FALSE.
+      ENDIF
+C---------------------------------------------------------------------
       EX = X
       ENU = ALPHA
       NCALC = MIN(NB,0)-2
--- a/libcruft/specfun/rybesl.f
+++ b/libcruft/specfun/rybesl.f
@@ -134,7 +134,9 @@
 C               Argonne, IL  60439
 C
 C----------------------------------------------------------------------
+      LOGICAL FIRST
       INTEGER I,K,NA,NB,NCALC
+      DOUBLE PRECISION D1MACH
       DOUBLE PRECISION
      1  ALFA,ALPHA,AYE,B,BY,C,CH,COSMU,D,DEL,DEN,DDIV,DIV,DMU,D1,D2,
      2  E,EIGHT,EN,ENU,EN1,EPS,EVEN,EX,F,FIVPI,G,GAMMA,H,HALF,ODD,
@@ -150,16 +152,14 @@
 C    PIBY2 = PI/2
 C    SQ2BPI = SQUARE ROOT OF 2/PI
 C----------------------------------------------------------------------
-      DATA ZERO,HALF,ONE,TWO,THREE/0.0D0,0.5D0,1.0D0,2.0D0,3.0D0/
-      DATA EIGHT,ONE5,TEN9/8.0D0,15.0D0,1.9D1/
-      DATA FIVPI,PIBY2/1.5707963267948966192D1,1.5707963267948966192D0/
-      DATA PI,SQ2BPI/3.1415926535897932385D0,7.9788456080286535588D-1/
-      DATA PIM5,ONBPI/7.0796326794896619231D-1,3.1830988618379067154D-1/
-C----------------------------------------------------------------------
-C  Machine-dependent constants
-C----------------------------------------------------------------------
-      DATA DEL,XMIN,XINF,EPS/1.0D-8,4.46D-308,1.79D308,1.11D-16/
-      DATA THRESH,XLARGE/16.0D0,1.0D8/
+      PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0)
+      PARAMETER (THREE = 3.0D0, EIGHT = 8.0D0, ONE5 = 1.5D1)
+      PARAMETER (TEN9 = 1.9D1, FIVPI = 1.5707963267948966192D1)
+      PARAMETER (PIBY2 = 1.5707963267948966192D0)
+      PARAMETER (PI = 3.1415926535897932385D0)
+      PARAMETER (SQ2BPI = 7.9788456080286535588D-1)
+      PARAMETER (PIM5 = /7.0796326794896619231D-1)
+      PARAMETER (ONBPI = 3.1830988618379067154D-1)
 C----------------------------------------------------------------------
 C  Coefficients for Chebyshev polynomial expansion of 
 C         1/gamma(1-x), abs(x) .le. .5
@@ -176,6 +176,19 @@
      9        -0.76852840844786673690D-01,-0.28387654227602353814D+00,
      A         0.92187029365045265648D+00/
 C----------------------------------------------------------------------
+C  Machine-dependent constants
+C----------------------------------------------------------------------
+      DATA FIRST /.TRUE./
+      IF (FIRST) THEN
+        EPS = D1MACH (4)
+        XINF = D1MACH (2)
+        XMIN = D1MACH (1)
+        DEL = SQRT (EPS)
+        XLARGE = ONE / DEL
+        THRESH = DINT (-LOG10 (EPS / TWO)) + ONE
+        FIRST = .FALSE.
+      ENDIF
+C----------------------------------------------------------------------
       EX = X
       ENU = ALPHA
       IF ((NB .GT. 0) .AND. (X .GE. XMIN) .AND. (EX .LT. XLARGE)