relibc_openlibm/slatec/d9knus.f

*DECK D9KNUS
      SUBROUTINE D9KNUS (XNU, X, BKNU, BKNU1, ISWTCH)
C***BEGIN PROLOGUE  D9KNUS
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      DOUBLE 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.50000E-01
C                                        with weighted error   2.16E-32
C                                         log weighted error  31.67
C                               significant figures required  30.86
C                                    decimal places required  32.40
C
C Series for ZNU1       on the interval -7.00000E-01 to  0.
C                                        with weighted error   2.45E-33
C                                         log weighted error  32.61
C                               significant figures required  31.85
C                                    decimal places required  33.26
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCSEVL, DGAMMA, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  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   900727  Added EXTERNAL statement.  (WRB)
C   920618  Removed space from variable names.  (RWC, WRB)
C***END PROLOGUE  D9KNUS
      DOUBLE PRECISION XNU, X, BKNU, BKNU1, ALPHA(32), BETA(32), A(32),
     1  C0KCS(29), ZNU1CS(20), ALNZ, ALN2, A0, BKNUD, BKNU0,
     2  B0, C0, EULER, EXPX, P1, P2, P3, QQ, RESULT, SQPI2, SQRTX, V,
     3  VLNZ, XI, XMU, XNUSML, XSML, X2N, X2TOV, Z, ZTOV, ALNSML,
     4  ALNBIG
      REAL ALNEPS
      DOUBLE PRECISION D1MACH, DCSEVL, DGAMMA
      LOGICAL FIRST
      EXTERNAL DGAMMA
      SAVE C0KCS, ZNU1CS, EULER, SQPI2, ALN2, NTC0K,
     1 NTZNU1, XNUSML, XSML, ALNSML, ALNBIG, ALNEPS, FIRST
      DATA C0KCS(  1) / +.6018305724 2626108387 5774451803 29 D-1     /
      DATA C0KCS(  2) / -.1536487143 3017286092 9597559431 24 D+0     /
      DATA C0KCS(  3) / -.1175117600 8210492040 0682292262 13 D-1     /
      DATA C0KCS(  4) / -.8524878889 1979509827 0484015509 87 D-3     /
      DATA C0KCS(  5) / -.6132983876 7496791874 0981769221 11 D-4     /
      DATA C0KCS(  6) / -.4405228124 5510444562 6798895485 05 D-5     /
      DATA C0KCS(  7) / -.3163124672 8384488192 9154458921 99 D-6     /
      DATA C0KCS(  8) / -.2271071938 2899588330 6737717933 96 D-7     /
      DATA C0KCS(  9) / -.1630564460 8077609552 2746205153 60 D-8     /
      DATA C0KCS( 10) / -.1170693929 9414776568 7560440431 30 D-9     /
      DATA C0KCS( 11) / -.8405206378 6464437174 5465934137 92 D-11    /
      DATA C0KCS( 12) / -.6034667011 8979991487 0960507371 98 D-12    /
      DATA C0KCS( 13) / -.4332696033 5681371952 0459973669 03 D-13    /
      DATA C0KCS( 14) / -.3110735803 0203546214 6346977722 37 D-14    /
      DATA C0KCS( 15) / -.2233407822 6736982254 4861334098 40 D-15    /
      DATA C0KCS( 16) / -.1603514671 6864226300 6357915286 10 D-16    /
      DATA C0KCS( 17) / -.1151271736 3666556196 0356977053 05 D-17    /
      DATA C0KCS( 18) / -.8265759174 6836959105 1694790892 58 D-19    /
      DATA C0KCS( 19) / -.5934548080 6383948172 3334366959 84 D-20    /
      DATA C0KCS( 20) / -.4260813819 6467143926 4996130239 76 D-21    /
      DATA C0KCS( 21) / -.3059126686 4812876299 2636983705 42 D-22    /
      DATA C0KCS( 22) / -.2196354142 6734575224 9755018155 16 D-23    /
      DATA C0KCS( 23) / -.1576911326 1495836071 1057506847 60 D-24    /
      DATA C0KCS( 24) / -.1132171393 5950320948 7577310480 56 D-25    /
      DATA C0KCS( 25) / -.8128624883 4598404082 7923497144 33 D-27    /
      DATA C0KCS( 26) / -.5836090089 3453226552 8293493159 49 D-28    /
      DATA C0KCS( 27) / -.4190124162 3610922519 4523377809 05 D-29    /
      DATA C0KCS( 28) / -.3008373796 0206435069 5305042128 62 D-30    /
      DATA C0KCS( 29) / -.2159915206 7808647728 3421680898 32 D-31    /
      DATA ZNU1CS(  1) / +.2033067569 9419172967 4444001216 911 D+0    /
      DATA ZNU1CS(  2) / +.1400779334 1321977106 2943670790 563 D+0    /
      DATA ZNU1CS(  3) / +.7916796961 0016135284 0972241972 320 D-2    /
      DATA ZNU1CS(  4) / +.3398011825 3210404535 2930092205 750 D-3    /
      DATA ZNU1CS(  5) / +.1174197568 8989336666 4507228352 690 D-4    /
      DATA ZNU1CS(  6) / +.3393575706 1226168033 3825865475 121 D-6    /
      DATA ZNU1CS(  7) / +.8425941769 7621991019 4629891264 803 D-8    /
      DATA ZNU1CS(  8) / +.1833366770 2485008918 4748150900 090 D-9    /
      DATA ZNU1CS(  9) / +.3549698447 0441631086 3007064469 557 D-11   /
      DATA ZNU1CS( 10) / +.6190324964 6988733220 5244342078 407 D-13   /
      DATA ZNU1CS( 11) / +.9819645356 8043942496 0346115456 527 D-15   /
      DATA ZNU1CS( 12) / +.1428513143 9649047421 1473563005 985 D-16   /
      DATA ZNU1CS( 13) / +.1918949218 8782529896 6162467488 436 D-18   /
      DATA ZNU1CS( 14) / +.2394309797 3949891416 2313140597 128 D-20   /
      DATA ZNU1CS( 15) / +.2788902468 1534735483 5870465474 995 D-22   /
      DATA ZNU1CS( 16) / +.3046066506 3303344258 2845214092 865 D-24   /
      DATA ZNU1CS( 17) / +.3131732370 4219181577 1564260932 089 D-26   /
      DATA ZNU1CS( 18) / +.3041330989 8785495164 5174908005 034 D-28   /
      DATA ZNU1CS( 19) / +.2798403846 3683308434 3185097659 733 D-30   /
      DATA ZNU1CS( 20) / +.2446371862 7449759648 5238794922 666 D-32   /
      DATA EULER / 0.5772156649 0153286060 6512090082 40D0 /
      DATA SQPI2 / +1.253314137 3155002512 0788264240 55 D0      /
      DATA ALN2 / 0.6931471805 5994530941 7232121458 18D0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  D9KNUS
      IF (FIRST) THEN
         ETA = 0.1D0*D1MACH(3)
         NTC0K = INITDS (C0KCS, 29, ETA)
         NTZNU1 = INITDS (ZNU1CS, 20, ETA)
C
         XNUSML = SQRT(D1MACH(3)/8.D0)
         XSML = 0.1D0*D1MACH(3)
         ALNSML = LOG (D1MACH(1))
         ALNBIG = LOG (D1MACH(2))
         ALNEPS = LOG (0.1D0*D1MACH(3))
      ENDIF
      FIRST = .FALSE.
C
      IF (XNU .LT. 0.D0 .OR. XNU .GE. 1.D0) CALL XERMSG ('SLATEC',
     +   'D9KNUS', 'XNU MUST BE GE 0 AND LT 1', 1, 2)
      IF (X .LE. 0.) CALL XERMSG ('SLATEC', 'D9KNUS', 'X MUST BE GT 0',
     +   2, 2)
C
      ISWTCH = 0
      IF (X.GT.2.0D0) 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.5D0) V = 1.0D0 - XNU
C
C CAREFULLY FIND (X/2)**XNU AND Z**XNU WHERE Z = X*X/4.
      ALNZ = 2.D0 * (LOG(X) - ALN2)
C
      IF (X.GT.XNU) GO TO 20
      IF (-0.5D0*XNU*ALNZ-ALN2-LOG(XNU) .GT. ALNBIG) CALL XERMSG
     +   ('SLATEC', 'D9KNUS', 'X SO SMALL BESSEL K-SUB-XNU OVERFLOWS',
     +   3, 2)
C
 20   VLNZ = V*ALNZ
      X2TOV = EXP (0.5D0*VLNZ)
      ZTOV = 0.0D0
      IF (VLNZ.GT.ALNSML) ZTOV = X2TOV**2
C
      A0 = 0.5D0*DGAMMA(1.0D0+V)
      B0 = 0.5D0*DGAMMA(1.0D0-V)
      C0 = -EULER
      IF (ZTOV.GT.0.5D0 .AND. V.GT.XNUSML) C0 = -0.75D0 +
     1  DCSEVL ((8.0D0*V)*V-1.0D0, C0KCS, NTC0K)
C
      IF (ZTOV.LE.0.5D0) ALPHA(1) = (A0-ZTOV*B0)/V
      IF (ZTOV.GT.0.5D0) ALPHA(1) = C0 - ALNZ*(0.75D0 +
     1  DCSEVL (VLNZ/0.35D0+1.0D0, ZNU1CS, NTZNU1))*B0
      BETA(1) = -0.5D0*(A0+ZTOV*B0)
C
      Z = 0.0D0
      IF (X.GT.XSML) Z = 0.25D0*X*X
      NTERMS = MAX (2.0, 11.0+(8.*REAL(ALNZ)-25.19-ALNEPS)
     1  /(4.28-REAL(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.0D0*XI*A0)/(XI*(XI+V))
        BETA(I) = (XI-0.5D0*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.5D0*(XNU+1.D0)*ALNZ-2.0D0*ALN2.GT.ALNBIG) ISWTCH = 1
      IF (ISWTCH.EQ.1) RETURN
      BKNUD = EXPX*BKNUD*2.0D0/(X2TOV*X)
C
      IF (XNU.LE.0.5D0) BKNU1 = V*BKNU/X - BKNUD
      IF (XNU.LE.0.5D0) RETURN
C
      BKNU0 = BKNU
      BKNU = -V*BKNU/X - BKNUD
      BKNU1 = 2.0D0*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.0D0/XSML) GO TO 90
      AN = -0.60 - 1.02/REAL(X)
      BN = -0.27 - 0.53/REAL(X)
      NTERMS = MIN (32, MAX1 (3.0, AN+BN*ALNEPS))
C
      DO 80 INU=1,2
        XMU = 0.D0
        IF (INU.EQ.1 .AND. XNU.GT.XNUSML) XMU = (4.0D0*XNU)*XNU
        IF (INU.EQ.2) XMU = 4.0D0*(ABS(XNU)+1.D0)**2
C
        A(1) = 1.0D0 - XMU
        A(2) = 9.0D0 - XMU
        A(3) = 25.0D0 - XMU
        IF (A(2).EQ.0.D0) RESULT = SQPI2*(16.D0*X+XMU+7.D0) /
     1    (16.D0*X*SQRTX)
        IF (A(2).EQ.0.D0) GO TO 70
C
        ALPHA(1) = 1.0D0
        ALPHA(2) = (16.D0*X+A(2))/A(2)
        ALPHA(3) = ((768.D0*X+48.D0*A(3))*X + A(2)*A(3))/(A(2)*A(3))
C
        BETA(1) = 1.0D0
        BETA(2) = (16.D0*X+(XMU+7.D0))/A(2)
        BETA(3) = ((768.D0*X+48.D0*(XMU+23.D0))*X +
     1    ((XMU+62.D0)*XMU+129.D0))/(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.D0)**2 - XMU
          QQ = 16.D0*X2N/A(I)
          P1 = -X2N*((12*N*N-20*N)-A(1))/((X2N-2.D0)*A(I))
     1      - QQ*X
          P2 = ((12*N*N-28*N+8)-A(1))/A(I) - QQ*X
          P3 = -X2N*A(I-3)/((X2N-2.D0)*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