relibc_openlibm/slatec/r9knus.f

*DECK R9KNUS
      SUBROUTINE R9KNUS (XNU, X, BKNU, BKNU1, ISWTCH)
C***BEGIN PROLOGUE  R9KNUS
C***SUBSIDIARY
C***PURPOSE  Compute Bessel functions EXP(X)*K-SUB-XNU(X) and EXP(X)*
C            K-SUB-XNU+1(X) for 0.0 .LE. XNU .LT. 1.0.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C10B3
C***TYPE      SINGLE PRECISION (R9KNUS-S, D9KNUS-D)
C***KEYWORDS  BESSEL FUNCTION, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Compute Bessel functions EXP(X) * K-sub-XNU (X)  and
C EXP(X) * K-sub-XNU+1 (X) for 0.0 .LE. XNU .LT. 1.0 .
C
C Series for C0K        on the interval  0.          to  2.50000D-01
C                                        with weighted error   1.60E-17
C                                         log weighted error  16.79
C                               significant figures required  15.99
C                                    decimal places required  17.40
C
C Series for ZNU1       on the interval -7.00000D-01 to  0.
C                                        with weighted error   1.43E-17
C                                         log weighted error  16.85
C                               significant figures required  16.08
C                                    decimal places required  17.38
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL, GAMMA, INITS, R1MACH, 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   900720  Routine changed from user-callable to subsidiary.  (WRB)
C   900727  Added EXTERNAL statement.  (WRB)
C   920618  Removed space from variable names.  (RWC, WRB)
C***END PROLOGUE  R9KNUS
      DIMENSION ALPHA(15), BETA(15), A(15), C0KCS(16), ZNU1CS(12)
      LOGICAL FIRST
      EXTERNAL GAMMA
      SAVE C0KCS, ZNU1CS, EULER, SQPI2, ALN2, NTC0K, NTZNU1,
     1 XNUSML, XSML, ALNSML, ALNBIG, ALNEPS, FIRST
      DATA C0KCS( 1) /    .0601830572 42626108E0 /
      DATA C0KCS( 2) /   -.1536487143 3017286E0 /
      DATA C0KCS( 3) /   -.0117511760 08210492E0 /
      DATA C0KCS( 4) /   -.0008524878 88919795E0 /
      DATA C0KCS( 5) /   -.0000613298 38767496E0 /
      DATA C0KCS( 6) /   -.0000044052 28124551E0 /
      DATA C0KCS( 7) /   -.0000003163 12467283E0 /
      DATA C0KCS( 8) /   -.0000000227 10719382E0 /
      DATA C0KCS( 9) /   -.0000000016 30564460E0 /
      DATA C0KCS(10) /   -.0000000001 17069392E0 /
      DATA C0KCS(11) /   -.0000000000 08405206E0 /
      DATA C0KCS(12) /   -.0000000000 00603466E0 /
      DATA C0KCS(13) /   -.0000000000 00043326E0 /
      DATA C0KCS(14) /   -.0000000000 00003110E0 /
      DATA C0KCS(15) /   -.0000000000 00000223E0 /
      DATA C0KCS(16) /   -.0000000000 00000016E0 /
      DATA ZNU1CS( 1) /    .2033067569 9419173E0 /
      DATA ZNU1CS( 2) /    .1400779334 1321977E0 /
      DATA ZNU1CS( 3) /    .0079167969 61001613E0 /
      DATA ZNU1CS( 4) /    .0003398011 82532104E0 /
      DATA ZNU1CS( 5) /    .0000117419 75688989E0 /
      DATA ZNU1CS( 6) /    .0000003393 57570612E0 /
      DATA ZNU1CS( 7) /    .0000000084 25941769E0 /
      DATA ZNU1CS( 8) /    .0000000001 83336677E0 /
      DATA ZNU1CS( 9) /    .0000000000 03549698E0 /
      DATA ZNU1CS(10) /    .0000000000 00061903E0 /
      DATA ZNU1CS(11) /    .0000000000 00000981E0 /
      DATA ZNU1CS(12) /    .0000000000 00000014E0 /
      DATA EULER / 0.5772156649 0153286E0 /
      DATA SQPI2 / 1.253314137 3155003E0 /
      DATA ALN2 / 0.693147180 55994531E0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  R9KNUS
      IF (FIRST) THEN
         NTC0K = INITS (C0KCS, 16, 0.1*R1MACH(3))
         NTZNU1 = INITS (ZNU1CS, 12, 0.1*R1MACH(3))
C
         XNUSML = SQRT (R1MACH(3)/8.0)
         XSML = 0.1*R1MACH(3)
         ALNSML = LOG (R1MACH(1))
         ALNBIG = LOG (R1MACH(2))
         ALNEPS = LOG (0.1*R1MACH(3))
      ENDIF
      FIRST = .FALSE.
C
      IF (XNU .LT. 0. .OR. XNU .GE. 1.0) CALL XERMSG ('SLATEC',
     +   'R9KNUS', 'XNU MUST BE GE 0 AND LT 1', 1, 2)
      IF (X .LE. 0.) CALL XERMSG ('SLATEC', 'R9KNUS', 'X MUST BE GT 0',
     +   2, 2)
C
      ISWTCH = 0
      IF (X.GT.2.0) GO TO 50
C
C X IS SMALL.  COMPUTE K-SUB-XNU (X) AND THE DERIVATIVE OF K-SUB-XNU (X)
C THEN FIND K-SUB-XNU+1 (X).  XNU IS REDUCED TO THE INTERVAL (-.5,+.5)
C THEN TO (0., .5), BECAUSE K OF NEGATIVE ORDER (-NU) = K OF POSITIVE
C ORDER (+NU).
C
      V = XNU
      IF (XNU.GT.0.5) V = 1.0 - XNU
C
C CAREFULLY FIND (X/2)**XNU AND Z**XNU WHERE Z = X*X/4.
      ALNZ = 2.0 * (LOG(X) - ALN2)
C
      IF (X.GT.XNU) GO TO 20
      IF (-0.5*XNU*ALNZ-ALN2-LOG(XNU) .GT. ALNBIG) CALL XERMSG
     +   ('SLATEC', 'R9KNUS', 'X SO SMALL BESSEL K-SUB-XNU OVERFLOWS',
     +   3, 2)
C
 20   VLNZ = V*ALNZ
      X2TOV = EXP (0.5*VLNZ)
      ZTOV = 0.0
      IF (VLNZ.GT.ALNSML) ZTOV = X2TOV**2
C
      A0 = 0.5*GAMMA(1.0+V)
      B0 = 0.5*GAMMA(1.0-V)
      C0 = -EULER
      IF (ZTOV.GT.0.5 .AND. V.GT.XNUSML) C0 = -0.75 +
     1  CSEVL ((8.0*V)*V-1., C0KCS, NTC0K)
C
      IF (ZTOV.LE.0.5) ALPHA(1) = (A0-ZTOV*B0)/V
      IF (ZTOV.GT.0.5) ALPHA(1) = C0 - ALNZ*(0.75 +
     1  CSEVL (VLNZ/0.35+1.0, ZNU1CS, NTZNU1))*B0
      BETA(1) = -0.5*(A0+ZTOV*B0)
C
      Z = 0.0
      IF (X.GT.XSML) Z = 0.25*X*X
      NTERMS = MAX (2.0, 11.0+(8.*ALNZ-25.19-ALNEPS)/(4.28-ALNZ))
      DO 30 I=2,NTERMS
        XI = I - 1
        A0 = A0/(XI*(XI-V))
        B0 = B0/(XI*(XI+V))
        ALPHA(I) = (ALPHA(I-1)+2.0*XI*A0)/(XI*(XI+V))
        BETA(I) = (XI-0.5*V)*ALPHA(I) - ZTOV*B0
 30   CONTINUE
C
      BKNU = ALPHA(NTERMS)
      BKNUD = BETA(NTERMS)
      DO 40 II=2,NTERMS
        I = NTERMS + 1 - II
        BKNU = ALPHA(I) + BKNU*Z
        BKNUD = BETA(I) + BKNUD*Z
 40   CONTINUE
C
      EXPX = EXP(X)
      BKNU = EXPX*BKNU/X2TOV
C
      IF (-0.5*(XNU+1.)*ALNZ-2.0*ALN2.GT.ALNBIG) ISWTCH = 1
      IF (ISWTCH.EQ.1) RETURN
      BKNUD = EXPX*BKNUD*2.0/(X2TOV*X)
C
      IF (XNU.LE.0.5) BKNU1 = V*BKNU/X - BKNUD
      IF (XNU.LE.0.5) RETURN
C
      BKNU0 = BKNU
      BKNU = -V*BKNU/X - BKNUD
      BKNU1 = 2.0*XNU*BKNU/X + BKNU0
      RETURN
C
C X IS LARGE.  FIND K-SUB-XNU (X) AND K-SUB-XNU+1 (X) WITH Y. L. LUKE-S
C RATIONAL EXPANSION.
C
 50   SQRTX = SQRT(X)
      IF (X.GT.1.0/XSML) GO TO 90
      AN = -1.56 + 4.0/X
      BN = -0.29 - 0.22/X
      NTERMS = MIN (15, MAX1 (3.0, AN+BN*ALNEPS))
C
      DO 80 INU=1,2
        XMU = 0.
        IF (INU.EQ.1 .AND. XNU.GT.XNUSML) XMU = (4.0*XNU)*XNU
        IF (INU.EQ.2) XMU = 4.0*(ABS(XNU)+1.)**2
C
        A(1) = 1.0 - XMU
        A(2) = 9.0 - XMU
        A(3) = 25.0 - XMU
        IF (A(2).EQ.0.) RESULT = SQPI2*(16.*X+XMU+7.)/(16.*X*SQRTX)
        IF (A(2).EQ.0.) GO TO 70
C
        ALPHA(1) = 1.0
        ALPHA(2) = (16.*X+A(2))/A(2)
        ALPHA(3) = ((768.*X+48.*A(3))*X + A(2)*A(3))/(A(2)*A(3))
C
        BETA(1) = 1.0
        BETA(2) = (16.*X+(XMU+7.))/A(2)
        BETA(3) = ((768.*X+48.*(XMU+23.))*X + ((XMU+62.)*XMU+129.))
     1    / (A(2)*A(3))
C
        IF (NTERMS.LT.4) GO TO 65
        DO 60 I=4,NTERMS
          N = I - 1
          X2N = 2*N - 1
C
          A(I) = (X2N+2.)**2 - XMU
          QQ = 16.*X2N/A(I)
          P1 = -X2N*(12*N*N-20*N-A(1))/((X2N-2.)*A(I)) - QQ*X
          P2 = (12*N*N-28*N+8-A(1))/A(I) - QQ*X
          P3 = -X2N*A(I-3)/((X2N-2.)*A(I))
C
          ALPHA(I) = -P1*ALPHA(I-1) - P2*ALPHA(I-2) - P3*ALPHA(I-3)
          BETA(I) = -P1*BETA(I-1) - P2*BETA(I-2) - P3*BETA(I-3)
 60     CONTINUE
C
 65     RESULT = SQPI2*BETA(NTERMS)/(SQRTX*ALPHA(NTERMS))
C
 70     IF (INU.EQ.1) BKNU = RESULT
        IF (INU.EQ.2) BKNU1 = RESULT
 80   CONTINUE
      RETURN
C
 90   BKNU = SQPI2/SQRTX
      BKNU1 = BKNU
      RETURN
C
      END