MINICERN/mathlib/gen/c/rpsipg64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:01  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 *FCA :          Fri Mar 26 17:27:50 CET 1999 by  Federico Carminati
  10 *               restore the entry for DDIGAM, even if obsolete
  11 #include "gen/pilot.h"
  12 #if defined(CERNLIB_DOUBLE)
  13       FUNCTION DDIGAM(X)
  14 #include "gen/imp64.inc"
  15       DDIGAM = DPSIPG(X,0)
  16       END
  17 C
  18       FUNCTION DPSIPG(X,K)
  19 C
  20 #include "gen/imp64.inc"
  21 C
  22       CHARACTER*(*) NAME
  23       PARAMETER(NAME='RPSIPG/DPSIPG')
  24 #endif
  25 #if !defined(CERNLIB_DOUBLE)
  26       FUNCTION RPSIPG(X,K)
  27 C
  28       CHARACTER*(*) NAME
  29       PARAMETER(NAME='RPSIPG')
  30 #endif
  31 C
  32       DIMENSION B(0:20,6),C(7,6),NB(6),P1(0:7),Q1(0:7),P2(0:4),Q2(0:4)
  33       DIMENSION SGN(6),SGF(0:6),SGH(6)
  34
  35       PARAMETER (DELTA = 1D-13)
  36       PARAMETER (Z1 = 1, HF = Z1/2)
  37       PARAMETER (PI = 3.14159 26535 89793 24D0)
  38       PARAMETER (C1 = -PI**2, C2 = 2*PI**3, C3 = 2*PI**4)
  39       PARAMETER (C4 = -8*PI**5, C5 = -8*PI**6, C6 = 16*PI**7)
  40
  41       CHARACTER*80 ERRTXT
  42
  43       DATA NB /16,17,17,18,19,20/
  44       DATA SGN /-1,1,-1,1,-1,1/, SGF /1,-1,2,-6,24,-120,720/
  45       DATA SGH /-0.5D0,1,-3,12,-60,360/
  46       DATA X0 /1.46163 21449 68362 34D0/
  47
  48       DATA (P1(J),Q1(J),J=0,7)
  49      A/ 1.35249 99667 72634 64D+4, 6.93891 11753 76344 44D-7,
  50      1  4.52856 01699 54728 97D+4, 1.97685 74263 04673 64D+4,
  51      2  4.51351 68469 73666 26D+4, 4.12551 60835 35383 23D+4,
  52      3  1.85290 11818 58261 02D+4, 2.93902 87119 93268 19D+4,
  53      4  3.32915 25149 40693 55D+3, 9.08196 66074 85517 03D+3,
  54      5  2.40680 32474 35720 18D+2, 1.24474 77785 67085 60D+3,
  55      6  5.15778 92000 13908 47D+0, 6.74291 29516 37859 38D+1,
  56      7  6.22835 06918 98474 58D-3, 1/
  57
  58       DATA (P2(J),Q2(J),J=0,4)
  59      A/-2.72817 57513 15296 78D-15,7.77788 54852 29616 04D+0,
  60      1 -6.48157 12376 61965 10D-1, 5.46117 73810 32150 70D+1,
  61      2 -4.48616 54391 80193 58D+0, 8.92920 70048 18613 70D+1,
  62      3 -7.01677 22776 67586 64D+0, 3.22703 49379 11433 61D+1,
  63      7 -2.12940 44513 10105 17D+0, 1/
  64
  65       DATA B( 0,1) / 0.33483 86979 10949 386D0/
  66       DATA B( 1,1) /-0.05518 74820 48730 095D0/
  67       DATA B( 2,1) / 0.00451 01907 36011 502D0/
  68       DATA B( 3,1) /-0.00036 57058 88303 721D0/
  69       DATA B( 4,1) / 0.00002 94346 27468 223D0/
  70       DATA B( 5,1) /-0.00000 23527 76815 151D0/
  71       DATA B( 6,1) / 0.00000 01868 53176 633D0/
  72       DATA B( 7,1) /-0.00000 00147 50720 184D0/
  73       DATA B( 8,1) / 0.00000 00011 57993 337D0/
  74       DATA B( 9,1) /-0.00000 00000 90439 179D0/
  75       DATA B(10,1) / 0.00000 00000 07029 627D0/
  76       DATA B(11,1) /-0.00000 00000 00543 989D0/
  77       DATA B(12,1) / 0.00000 00000 00041 925D0/
  78       DATA B(13,1) /-0.00000 00000 00003 219D0/
  79       DATA B(14,1) / 0.00000 00000 00000 246D0/
  80       DATA B(15,1) /-0.00000 00000 00000 019D0/
  81       DATA B(16,1) / 0.00000 00000 00000 001D0/
  82
  83       DATA B( 0,2) /-0.11259 29353 45473 830D0/
  84       DATA B( 1,2) / 0.03655 70017 42820 941D0/
  85       DATA B( 2,2) /-0.00443 59424 96027 282D0/
  86       DATA B( 3,2) / 0.00047 54758 54728 926D0/
  87       DATA B( 4,2) /-0.00004 74718 36382 632D0/
  88       DATA B( 5,2) / 0.00000 45218 15237 353D0/
  89       DATA B( 6,2) /-0.00000 04163 00079 620D0/
  90       DATA B( 7,2) / 0.00000 00373 38998 165D0/
  91       DATA B( 8,2) /-0.00000 00032 79914 474D0/
  92       DATA B( 9,2) / 0.00000 00002 83211 377D0/
  93       DATA B(10,2) /-0.00000 00000 24104 028D0/
  94       DATA B(11,2) / 0.00000 00000 02026 297D0/
  95       DATA B(12,2) /-0.00000 00000 00168 524D0/
  96       DATA B(13,2) / 0.00000 00000 00013 885D0/
  97       DATA B(14,2) /-0.00000 00000 00001 135D0/
  98       DATA B(15,2) / 0.00000 00000 00000 092D0/
  99       DATA B(16,2) /-0.00000 00000 00000 007D0/
 100       DATA B(17,2) / 0.00000 00000 00000 001D0/
 101
 102       DATA B( 0,3) / 0.07601 26046 55110 384D0/
 103       DATA B( 1,3) /-0.03625 71864 81828 739D0/
 104       DATA B( 2,3) / 0.00579 72023 38937 002D0/
 105       DATA B( 3,3) /-0.00076 96465 13610 481D0/
 106       DATA B( 4,3) / 0.00009 14920 82189 884D0/
 107       DATA B( 5,3) /-0.00001 00971 31488 364D0/
 108       DATA B( 6,3) / 0.00000 10557 77442 831D0/
 109       DATA B( 7,3) /-0.00000 01059 29577 481D0/
 110       DATA B( 8,3) / 0.00000 00102 85494 201D0/
 111       DATA B( 9,3) /-0.00000 00009 72314 310D0/
 112       DATA B(10,3) / 0.00000 00000 89884 635D0/
 113       DATA B(11,3) /-0.00000 00000 08153 171D0/
 114       DATA B(12,3) / 0.00000 00000 00727 572D0/
 115       DATA B(13,3) /-0.00000 00000 00064 010D0/
 116       DATA B(14,3) / 0.00000 00000 00005 562D0/
 117       DATA B(15,3) /-0.00000 00000 00000 478D0/
 118       DATA B(16,3) / 0.00000 00000 00000 041D0/
 119       DATA B(17,3) /-0.00000 00000 00000 003D0/
 120
 121       DATA B( 0,4) /-0.07723 47240 56994 793D0/
 122       DATA B( 1,4) / 0.04786 71634 51599 467D0/
 123       DATA B( 2,4) /-0.00944 07021 86674 632D0/
 124       DATA B( 3,4) / 0.00148 95447 40103 448D0/
 125       DATA B( 4,4) /-0.00020 49440 23348 860D0/
 126       DATA B( 5,4) / 0.00002 56714 25065 297D0/
 127       DATA B( 6,4) /-0.00000 30013 93581 584D0/
 128       DATA B( 7,4) / 0.00000 03327 66437 356D0/
 129       DATA B( 8,4) /-0.00000 00353 65412 111D0/
 130       DATA B( 9,4) / 0.00000 00036 30622 927D0/
 131       DATA B(10,4) /-0.00000 00003 62096 951D0/
 132       DATA B(11,4) / 0.00000 00000 35237 509D0/
 133       DATA B(12,4) /-0.00000 00000 03357 440D0/
 134       DATA B(13,4) / 0.00000 00000 00314 068D0/
 135       DATA B(14,4) /-0.00000 00000 00028 908D0/
 136       DATA B(15,4) / 0.00000 00000 00002 623D0/
 137       DATA B(16,4) /-0.00000 00000 00000 235D0/
 138       DATA B(17,4) / 0.00000 00000 00000 021D0/
 139       DATA B(18,4) /-0.00000 00000 00000 002D0/
 140
 141       DATA B( 0,5) / 0.10493 30344 59278 632D0/
 142       DATA B( 1,5) /-0.07887 79016 52793 557D0/
 143       DATA B( 2,5) / 0.01839 74151 12159 397D0/
 144       DATA B( 3,5) /-0.00335 22841 59396 504D0/
 145       DATA B( 4,5) / 0.00052 28782 30918 016D0/
 146       DATA B( 5,5) /-0.00007 31797 85814 740D0/
 147       DATA B( 6,5) / 0.00000 94497 29612 085D0/
 148       DATA B( 7,5) /-0.00000 11463 39856 723D0/
 149       DATA B( 8,5) / 0.00000 01322 69366 108D0/
 150       DATA B( 9,5) /-0.00000 00146 46669 180D0/
 151       DATA B(10,5) / 0.00000 00015 66940 742D0/
 152       DATA B(11,5) /-0.00000 00001 62791 157D0/
 153       DATA B(12,5) / 0.00000 00000 16490 345D0/
 154       DATA B(13,5) /-0.00000 00000 01634 028D0/
 155       DATA B(14,5) / 0.00000 00000 00158 807D0/
 156       DATA B(15,5) /-0.00000 00000 00015 171D0/
 157       DATA B(16,5) / 0.00000 00000 00001 427D0/
 158       DATA B(17,5) /-0.00000 00000 00000 132D0/
 159       DATA B(18,5) / 0.00000 00000 00000 012D0/
 160       DATA B(19,5) /-0.00000 00000 00000 001D0/
 161
 162       DATA B( 0,6) /-0.17861 76221 42502 753D0/
 163       DATA B( 1,6) / 0.15577 64622 00520 579D0/
 164       DATA B( 2,6) /-0.04172 36376 73831 277D0/
 165       DATA B( 3,6) / 0.00859 71413 03245 400D0/
 166       DATA B( 4,6) /-0.00149 62277 61073 229D0/
 167       DATA B( 5,6) / 0.00023 10896 08557 137D0/
 168       DATA B( 6,6) /-0.00003 26320 44778 436D0/
 169       DATA B( 7,6) / 0.00000 42960 97867 090D0/
 170       DATA B( 8,6) /-0.00000 05345 28790 204D0/
 171       DATA B( 9,6) / 0.00000 00634 78151 644D0/
 172       DATA B(10,6) /-0.00000 00072 48699 714D0/
 173       DATA B(11,6) / 0.00000 00008 00521 979D0/
 174       DATA B(12,6) /-0.00000 00000 85888 793D0/
 175       DATA B(13,6) / 0.00000 00000 08985 442D0/
 176       DATA B(14,6) /-0.00000 00000 00919 356D0/
 177       DATA B(15,6) / 0.00000 00000 00092 225D0/
 178       DATA B(16,6) /-0.00000 00000 00009 090D0/
 179       DATA B(17,6) / 0.00000 00000 00000 882D0/
 180       DATA B(18,6) /-0.00000 00000 00000 084D0/
 181       DATA B(19,6) / 0.00000 00000 00000 008D0/
 182       DATA B(20,6) /-0.00000 00000 00000 001D0/
 183
 184       DATA C(1,1) / 1.66666 66666 66666 67D-1/
 185       DATA C(2,1) /-3.33333 33333 33333 33D-2/
 186       DATA C(3,1) / 2.38095 23809 52380 95D-2/
 187       DATA C(4,1) /-3.33333 33333 33333 33D-2/
 188       DATA C(5,1) / 7.57575 75757 57575 76D-2/
 189       DATA C(6,1) /-2.53113 55311 35531 14D-1/
 190       DATA C(7,1) / 1.16666 66666 66666 67D 0/
 191
 192       DATA C(1,2) / 5.00000 00000 00000 00D-1/
 193       DATA C(2,2) /-1.66666 66666 66666 67D-1/
 194       DATA C(3,2) / 1.66666 66666 66666 67D-1/
 195       DATA C(4,2) /-3.00000 00000 00000 00D-1/
 196       DATA C(5,2) / 8.33333 33333 33333 33D-1/
 197       DATA C(6,2) /-3.29047 61904 76190 48D 0/
 198       DATA C(7,2) / 1.75000 00000 00000 00D 1/
 199
 200       DATA C(1,3) / 2.00000 00000 00000 00D 0/
 201       DATA C(2,3) /-1.00000 00000 00000 00D 0/
 202       DATA C(3,3) / 1.33333 33333 33333 33D 0/
 203       DATA C(4,3) /-3.00000 00000 00000 00D 0/
 204       DATA C(5,3) / 1.00000 00000 00000 00D+1/
 205       DATA C(6,3) /-4.60666 66666 66666 67D+1/
 206       DATA C(7,3) / 2.80000 00000 00000 00D+2/
 207
 208       DATA (C(J,4),J=1,7) /10,-7,12,-33,130,-691,4760/
 209       DATA (C(J,5),J=1,7) /60,-56,120,-396,1820,-11056,85680/
 210       DATA (C(J,6),J=1,7) /420,-504,1320,-5148,27300,-187952,1627920/
 211
 212       A=ABS(X)
 213       V=A
 214       IX=X-DELTA
 215       IF(K .LT. 0 .OR. K .GT. 6) THEN
 216        H=0
 217        WRITE(ERRTXT,101) K
 218        CALL MTLPRT(NAME,'C316.1',ERRTXT)
 219       ELSEIF(ABS(IX-X) .LE. DELTA) THEN
 220        H=0
 221        WRITE(ERRTXT,102) X
 222        CALL MTLPRT(NAME,'C316.2',ERRTXT)
 223       ELSEIF(K .EQ. 0) THEN
 224        IF(A .LE. 3) THEN
 225         S=0
 226         IF(A .LT. HF) THEN
 227          S=1/V
 228          V=V+1
 229         ENDIF
 230         AP=P1(7)
 231         AQ=Q1(7)
 232         DO 11 I = 6,0,-1
 233         AP=P1(I)+V*AP
 234    11   AQ=Q1(I)+V*AQ
 235         H=(V-X0)*AP/AQ-S
 236        ELSE
 237         R=1/V**2
 238         AP=P2(4)
 239         AQ=Q2(4)
 240         DO 12 I = 3,0,-1
 241         AP=P2(I)+R*AP
 242    12   AQ=Q2(I)+R*AQ
 243         H=LOG(V)-HF/V+AP/AQ
 244        ENDIF
 245        IF(X .LT. 0) H=H+1/A+PI/TAN(PI*A)
 246       ELSE
 247        K1=K+1
 248        IF(A .LE. 10) THEN
 249         IF(A .LT. 3) THEN
 250          S=-1/V**K1
 251          DO 1 J = 1,2-INT(A)
 252          V=V+1
 253     1    S=S-1/V**K1
 254          V=V+1
 255         ELSEIF(A .LE. 4) THEN
 256          S=0
 257         ELSE
 258          V=V-1
 259          S=1/V**K1
 260          DO 5 J = 1,INT(A)-4
 261          V=V-1
 262     5    S=S+1/V**K1
 263         ENDIF
 264         H=2*V-7
 265         ALFA=H+H
 266         B1=0
 267         B2=0
 268         DO 2 J = NB(K),0,-1
 269         B0=B(J,K)+ALFA*B1-B2
 270         B2=B1
 271     2   B1=B0
 272         H=B0-H*B2+SGF(K)*S
 273        ELSE
 274         S=0
 275         IF(A .LT. 15) THEN
 276          S=1/V**K1
 277          DO 3 J = 1,14-INT(A)
 278          V=V+1
 279     3    S=S+1/V**K1
 280          V=V+1
 281         ENDIF
 282         R=1/V**2
 283         P=R*C(7,K)
 284         DO 4 J = 6,1,-1
 285     4   P=R*(C(J,K)+P)
 286         H=((SGF(K-1)-SGN(K)*P)*V-SGH(K))/V**K1-SGF(K)*S
 287        ENDIF
 288        IF(X .LT. 0) THEN
 289         P=PI*A
 290         IF(K .EQ. 1) THEN
 291          V=C1/SIN(P)**2
 292         ELSEIF(K .EQ. 2) THEN
 293          V=C2*COS(P)/SIN(P)**3
 294         ELSEIF(K .EQ. 3) THEN
 295          S=SIN(P)**2
 296          V=C3*(2*S-3)/S**2
 297         ELSEIF(K .EQ. 4) THEN
 298          S=SIN(P)
 299          V=C4*COS(P)*(S**2-3)/S**5
 300         ELSEIF(K .EQ. 5) THEN
 301          S=SIN(P)**2
 302          V=C5*(15-15*S+2*S**2)/S**3
 303         ELSEIF(K .EQ. 6) THEN
 304          S=SIN(P)
 305          V=C6*COS(P)*(45-30*S**2+2*S**4)/S**7
 306         ENDIF
 307         H=SGN(K)*(H+V+SGF(K)/A**K1)
 308        ENDIF
 309       ENDIF
 310 #if defined(CERNLIB_DOUBLE)
 311       DPSIPG=H
 312 #endif
 313 #if !defined(CERNLIB_DOUBLE)
 314       RPSIPG=H
 315 #endif
 316       RETURN
 317   101 FORMAT('K = ',I5,'  (< 0  OR  > 6)')
 318   102 FORMAT('ARGUMENT EQUALS NON-POSITIVE INTEGER =',1P,E15.6)
 319       END