Home Explore Help
Sign In
DragonOS-Community
/
relibc_openlibm
mirror of https://gitlab.redox-os.org/redox-os/openlibm.git
2
0
Fork 0
Files Issues 0 Wiki
Tree: 50ce8eff11
Branches Tags
relibc_openlibm
/
slatec
/
dpsi.f

dpsi.f 7.5 KB
History Raw

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163 
  1. *DECK DPSI
    2. 

  2.       DOUBLE PRECISION FUNCTION DPSI (X)
    3. 

  3. C***BEGIN PROLOGUE  DPSI
    4. 

  4. C***PURPOSE  Compute the Psi (or Digamma) function.
    5. 

  5. C***LIBRARY   SLATEC (FNLIB)
    6. 

  6. C***CATEGORY  C7C
    7. 

  7. C***TYPE      DOUBLE PRECISION (PSI-S, DPSI-D, CPSI-C)
    8. 

  8. C***KEYWORDS  DIGAMMA FUNCTION, FNLIB, PSI FUNCTION, SPECIAL FUNCTIONS
    9. 

  9. C***AUTHOR  Fullerton, W., (LANL)
    10. 

  10. C***DESCRIPTION
    11. 

  11. C
    12. 

  12. C DPSI calculates the double precision Psi (or Digamma) function for
    13. 

  13. C double precision argument X.  PSI(X) is the logarithmic derivative
    14. 

  14. C of the Gamma function of X.
    15. 

  15. C
    16. 

  16. C Series for PSI        on the interval  0.          to  1.00000E+00
    17. 

  17. C                                        with weighted error   5.79E-32
    18. 

  18. C                                         log weighted error  31.24
    19. 

  19. C                               significant figures required  30.93
    20. 

  20. C                                    decimal places required  32.05
    21. 

  21. C
    22. 

  22. C
    23. 

  23. C Series for APSI       on the interval  0.          to  1.00000E-02
    24. 

  24. C                                        with weighted error   7.75E-33
    25. 

  25. C                                         log weighted error  32.11
    26. 

  26. C                               significant figures required  28.88
    27. 

  27. C                                    decimal places required  32.71
    28. 

  28. C
    29. 

  29. C***REFERENCES  (NONE)
    30. 

  30. C***ROUTINES CALLED  D1MACH, DCOT, DCSEVL, INITDS, XERMSG
    31. 

  31. C***REVISION HISTORY  (YYMMDD)
    32. 

  32. C   770601  DATE WRITTEN
    33. 

  33. C   890531  Changed all specific intrinsics to generic.  (WRB)
    34. 

  34. C   890911  Removed unnecessary intrinsics.  (WRB)
    35. 

  35. C   890911  REVISION DATE from Version 3.2
    36. 

  36. C   891214  Prologue converted to Version 4.0 format.  (BAB)
    37. 

  37. C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    38. 

  38. C   900727  Added EXTERNAL statement.  (WRB)
    39. 

  39. C   920618  Removed space from variable name.  (RWC, WRB)
    40. 

  40. C***END PROLOGUE  DPSI
    41. 

  41.       DOUBLE PRECISION X, PSICS(42), APSICS(16), AUX, DXREL, PI, XBIG,
    42. 

  42.      1  Y, DCOT, DCSEVL, D1MACH
    43. 

  43.       LOGICAL FIRST
    44. 

  44.       EXTERNAL DCOT
    45. 

  45.       SAVE PSICS, APSICS, PI, NTPSI, NTAPSI, XBIG, DXREL, FIRST
    46. 

  46.       DATA PSICS(  1) / -.3805708083 5217921520 4376776670 39 D-1     /
    47. 

  47.       DATA PSICS(  2) / +.4914153930 2938712748 2046996542 77 D+0     /
    48. 

  48.       DATA PSICS(  3) / -.5681574782 1244730242 8920647340 81 D-1     /
    49. 

  49.       DATA PSICS(  4) / +.8357821225 9143131362 7756507478 62 D-2     /
    50. 

  50.       DATA PSICS(  5) / -.1333232857 9943425998 0792741723 93 D-2     /
    51. 

  51.       DATA PSICS(  6) / +.2203132870 6930824892 8723979795 21 D-3     /
    52. 

  52.       DATA PSICS(  7) / -.3704023817 8456883592 8890869492 29 D-4     /
    53. 

  53.       DATA PSICS(  8) / +.6283793654 8549898933 6514187176 90 D-5     /
    54. 

  54.       DATA PSICS(  9) / -.1071263908 5061849855 2835417470 74 D-5     /
    55. 

  55.       DATA PSICS( 10) / +.1831283946 5484165805 7315898103 78 D-6     /
    56. 

  56.       DATA PSICS( 11) / -.3135350936 1808509869 0057797968 85 D-7     /
    57. 

  57.       DATA PSICS( 12) / +.5372808776 2007766260 4719191436 15 D-8     /
    58. 

  58.       DATA PSICS( 13) / -.9211681415 9784275717 8806326247 30 D-9     /
    59. 

  59.       DATA PSICS( 14) / +.1579812652 1481822782 2528840328 23 D-9     /
    60. 

  60.       DATA PSICS( 15) / -.2709864613 2380443065 4405894097 07 D-10    /
    61. 

  61.       DATA PSICS( 16) / +.4648722859 9096834872 9473195295 49 D-11    /
    62. 

  62.       DATA PSICS( 17) / -.7975272563 8303689726 5047977727 37 D-12    /
    63. 

  63.       DATA PSICS( 18) / +.1368272385 7476992249 2510538928 38 D-12    /
    64. 

  64.       DATA PSICS( 19) / -.2347515606 0658972717 3206779807 19 D-13    /
    65. 

  65.       DATA PSICS( 20) / +.4027630715 5603541107 9079250062 81 D-14    /
    66. 

  66.       DATA PSICS( 21) / -.6910251853 1179037846 5474229747 71 D-15    /
    67. 

  67.       DATA PSICS( 22) / +.1185604713 8863349552 9291395257 68 D-15    /
    68. 

  68.       DATA PSICS( 23) / -.2034168961 6261559308 1542104842 23 D-16    /
    69. 

  69.       DATA PSICS( 24) / +.3490074968 6463043850 3742329323 51 D-17    /
    70. 

  70.       DATA PSICS( 25) / -.5988014693 4976711003 0110813934 93 D-18    /
    71. 

  71.       DATA PSICS( 26) / +.1027380162 8080588258 3980057122 13 D-18    /
    72. 

  72.       DATA PSICS( 27) / -.1762704942 4561071368 3592601053 86 D-19    /
    73. 

  73.       DATA PSICS( 28) / +.3024322801 8156920457 4540354901 33 D-20    /
    74. 

  74.       DATA PSICS( 29) / -.5188916830 2092313774 2860888746 66 D-21    /
    75. 

  75.       DATA PSICS( 30) / +.8902773034 5845713905 0058874879 99 D-22    /
    76. 

  76.       DATA PSICS( 31) / -.1527474289 9426728392 8949719040 00 D-22    /
    77. 

  77.       DATA PSICS( 32) / +.2620731479 8962083136 3583180799 99 D-23    /
    78. 

  78.       DATA PSICS( 33) / -.4496464273 8220696772 5983880533 33 D-24    /
    79. 

  79.       DATA PSICS( 34) / +.7714712959 6345107028 9193642666 66 D-25    /
    80. 

  80.       DATA PSICS( 35) / -.1323635476 1887702968 1026389333 33 D-25    /
    81. 

  81.       DATA PSICS( 36) / +.2270999436 2408300091 2773119999 99 D-26    /
    82. 

  82.       DATA PSICS( 37) / -.3896419021 5374115954 4913919999 99 D-27    /
    83. 

  83.       DATA PSICS( 38) / +.6685198138 8855302310 6798933333 33 D-28    /
    84. 

  84.       DATA PSICS( 39) / -.1146998665 4920864872 5299199999 99 D-28    /
    85. 

  85.       DATA PSICS( 40) / +.1967938588 6541405920 5154133333 33 D-29    /
    86. 

  86.       DATA PSICS( 41) / -.3376448818 9750979801 9072000000 00 D-30    /
    87. 

  87.       DATA PSICS( 42) / +.5793070319 3214159246 6773333333 33 D-31    /
    88. 

  88.       DATA APSICS(  1) / -.8327107910 6929076017 4456932269 D-3        /
    89. 

  89.       DATA APSICS(  2) / -.4162518421 9273935282 1627121990 D-3        /
    90. 

  90.       DATA APSICS(  3) / +.1034315609 7874129117 4463193961 D-6        /
    91. 

  91.       DATA APSICS(  4) / -.1214681841 3590415298 7299556365 D-9        /
    92. 

  92.       DATA APSICS(  5) / +.3113694319 9835615552 1240278178 D-12       /
    93. 

  93.       DATA APSICS(  6) / -.1364613371 9317704177 6516100945 D-14       /
    94. 

  94.       DATA APSICS(  7) / +.9020517513 1541656513 0837974000 D-17       /
    95. 

  95.       DATA APSICS(  8) / -.8315429974 2159146482 9933635466 D-19       /
    96. 

  96.       DATA APSICS(  9) / +.1012242570 7390725418 8479482666 D-20       /
    97. 

  97.       DATA APSICS( 10) / -.1562702494 3562250762 0478933333 D-22       /
    98. 

  98.       DATA APSICS( 11) / +.2965427168 0890389613 3226666666 D-24       /
    99. 

  99.       DATA APSICS( 12) / -.6746868867 6570216374 1866666666 D-26       /
    100. 

  100.       DATA APSICS( 13) / +.1803453116 9718990421 3333333333 D-27       /
    101. 

  101.       DATA APSICS( 14) / -.5569016182 4598360746 6666666666 D-29       /
    102. 

  102.       DATA APSICS( 15) / +.1958679226 0773625173 3333333333 D-30       /
    103. 

  103.       DATA APSICS( 16) / -.7751958925 2333568000 0000000000 D-32       /
    104. 

  104.       DATA PI / 3.1415926535 8979323846 2643383279 50 D0 /
    105. 

  105.       DATA FIRST /.TRUE./
    106. 

  106. C***FIRST EXECUTABLE STATEMENT  DPSI
    107. 

  107.       IF (FIRST) THEN
    108. 

  108.          NTPSI = INITDS (PSICS, 42, 0.1*REAL(D1MACH(3)) )
    109. 

  109.          NTAPSI = INITDS (APSICS, 16, 0.1*REAL(D1MACH(3)) )
    110. 

  110. C
    111. 

  111.          XBIG = 1.0D0/SQRT(D1MACH(3))
    112. 

  112.          DXREL = SQRT(D1MACH(4))
    113. 

  113.       ENDIF
    114. 

  114.       FIRST = .FALSE.
    115. 

  115. C
    116. 

  116.       Y = ABS(X)
    117. 

  117. C
    118. 

  118.       IF (Y.GT.10.0D0) GO TO 50
    119. 

  119. C
    120. 

  120. C DPSI(X) FOR ABS(X) .LE. 2
    121. 

  121. C
    122. 

  122.       N = X
    123. 

  123.       IF (X.LT.0.D0) N = N - 1
    124. 

  124.       Y = X - N
    125. 

  125.       N = N - 1
    126. 

  126.       DPSI = DCSEVL (2.D0*Y-1.D0, PSICS, NTPSI)
    127. 

  127.       IF (N.EQ.0) RETURN
    128. 

  128. C
    129. 

  129.       IF (N.GT.0) GO TO 30
    130. 

  130. C
    131. 

  131.       N = -N
    132. 

  132.       IF (X .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DPSI', 'X IS 0', 2, 2)
    133. 

  133.       IF (X .LT. 0.D0 .AND. X+N-2 .EQ. 0.D0) CALL XERMSG ('SLATEC',
    134. 

  134.      +   'DPSI', 'X IS A NEGATIVE INTEGER', 3, 2)
    135. 

  135.       IF (X .LT. (-0.5D0) .AND. ABS((X-AINT(X-0.5D0))/X) .LT. DXREL)
    136. 

  136.      +   CALL XERMSG ('SLATEC', 'DPSI',
    137. 

  137.      +   'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',
    138. 

  138.      +   1, 1)
    139. 

  139. C
    140. 

  140.       DO 20 I=1,N
    141. 

  141.         DPSI = DPSI - 1.D0/(X+I-1)
    142. 

  142.  20   CONTINUE
    143. 

  143.       RETURN
    144. 

  144. C
    145. 

  145. C DPSI(X) FOR X .GE. 2.0 AND X .LE. 10.0
    146. 

  146. C
    147. 

  147.  30   DO 40 I=1,N
    148. 

  148.         DPSI = DPSI + 1.0D0/(Y+I)
    149. 

  149.  40   CONTINUE
    150. 

  150.       RETURN
    151. 

  151. C
    152. 

  152. C DPSI(X) FOR ABS(X) .GT. 10.0
    153. 

  153. C
    154. 

  154.  50   AUX = 0.D0
    155. 

  155.       IF (Y.LT.XBIG) AUX = DCSEVL (2.D0*(10.D0/Y)**2-1.D0, APSICS,
    156. 

  156.      1  NTAPSI)
    157. 

  157. C
    158. 

  158.       IF (X.LT.0.D0) DPSI = LOG(ABS(X)) - 0.5D0/X + AUX
    159. 

  159.      1  - PI*DCOT(PI*X)
    160. 

  160.       IF (X.GT.0.D0) DPSI = LOG(X) - 0.5D0/X + AUX
    161. 

  161.       RETURN
    162. 

  162. C
    163. 

  163.       END
    164.