Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/cdtr.ssp
There are 2 other files named cdtr.ssp in the archive. Click here to see a list.
C                                                                       CDTR  10
C     ..................................................................CDTR  20
C                                                                       CDTR  30
C        SUBROUTINE CDTR                                                CDTR  40
C                                                                       CDTR  50
C        PURPOSE                                                        CDTR  60
C           COMPUTES P(X) = PROBABILITY THAT THE RANDOM VARIABLE U,     CDTR  70
C           DISTRIBUTED ACCORDING TO THE CHI-SQUARE DISTRIBUTION WITH G CDTR  80
C           DEGREES OF FREEDOM, IS LESS THAN OR EQUAL TO X.  F(G,X), THECDTR  90
C           ORDINATE OF THE CHI-SQUARE DENSITY AT X, IS ALSO COMPUTED.  CDTR 100
C                                                                       CDTR 110
C        USAGE                                                          CDTR 120
C           CALL CDTR(X,G,P,D,IER)                                      CDTR 130
C                                                                       CDTR 140
C        DESCRIPTION OF PARAMETERS                                      CDTR 150
C           X   - INPUT SCALAR FOR WHICH P(X) IS COMPUTED.              CDTR 160
C           G   - NUMBER OF DEGREES OF FREEDOM OF THE CHI-SQUARE        CDTR 170
C                 DISTRIBUTION.  G IS A CONTINUOUS PARAMETER.           CDTR 180
C           P   - OUTPUT PROBABILITY.                                   CDTR 190
C           D   - OUTPUT DENSITY.                                       CDTR 200
C           IER - RESULTANT ERROR CODE WHERE                            CDTR 210
C               IER= 0 --- NO ERROR                                     CDTR 220
C               IER=-1 --- AN INPUT PARAMETER IS INVALID.  X IS LESS    CDTR 230
C                          THAN 0.0, OR G IS LESS THAN 0.5 OR GREATER   CDTR 240
C                          THAN 2*10**(+5).  P AND D ARE SET TO -1.7E38. CDTR 250
C               IER=+1 --- INVALID OUTPUT.  P IS LESS THAN ZERO OR      CDTR 260
C                          GREATER THAN ONE, OR SERIES FOR T1 (SEE      CDTR 270
C                          MATHEMATICAL DESCRIPTION) HAS FAILED TO      CDTR 280
C                          CONVERGE.  P IS SET TO 1.7E38.                CDTR 290
C                                                                       CDTR 300
C        REMARKS                                                        CDTR 310
C           SEE MATHEMATICAL DESCRIPTION.                               CDTR 320
C                                                                       CDTR 330
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  CDTR 340
C           DLGAM                                                       CDTR 350
C           NDTR                                                        CDTR 360
C                                                                       CDTR 370
C        METHOD                                                         CDTR 380
C           REFER TO R.E. BARGMANN AND S.P. GHOSH, STATISTICAL          CDTR 390
C           DISTRIBUTION PROGRAMS FOR A COMPUTER LANGUAGE,              CDTR 400
C           IBM RESEARCH REPORT RC-1094, 1963.                          CDTR 410
C                                                                       CDTR 420
C     ..................................................................CDTR 430
C                                                                       CDTR 440
      SUBROUTINE CDTR(X,G,P,D,IER)                                      CDTR 450
      DOUBLE PRECISION XX,DLXX,X2,DLX2,GG,G2,DLT3,THETA,THP1,           CDTR 460
     1GLG2,DD,T11,SER,CC,XI,FAC,TLOG,TERM,GTH,A2,A,B,C,DT2,DT3,THPI     CDTR 470
C                                                                       CDTR 480
C        TEST FOR VALID INPUT DATA                                      CDTR 490
C                                                                       CDTR 500
      IF(G-(.5-1.E-5)) 590,10,10                                        CDTR 510
   10 IF(G-2.E+5) 20,20,590                                             CDTR 520
   20 IF(X) 590,30,30                                                   CDTR 530
C                                                                       CDTR 540
C        TEST FOR X NEAR 0.0                                            CDTR 550
C                                                                       CDTR 560
   30 IF(X-1.E-8) 40,40,80                                              CDTR 570
   40 P=0.0                                                             CDTR 580
      IF(G-2.) 50,60,70                                                 CDTR 590
   50 D=1.7E38                                                           CDTR 600
      GO TO 610                                                         CDTR 610
   60 D=0.5                                                             CDTR 620
      GO TO 610                                                         CDTR 630
   70 D=0.0                                                             CDTR 640
      GO TO 610                                                         CDTR 650
C                                                                       CDTR 660
C        TEST FOR X GREATER THAN 1.E+6                                  CDTR 670
C                                                                       CDTR 680
   80 IF(X-1.E+6) 100,100,90                                            CDTR 690
   90 D=0.0                                                             CDTR 700
      P=1.0                                                             CDTR 710
      GO TO 610                                                         CDTR 720
C                                                                       CDTR 730
C        SET PROGRAM PARAMETERS                                         CDTR 740
C                                                                       CDTR 750
  100 XX=DBLE(X)                                                        CDTR 760
      DLXX=DLOG(XX)                                                     CDTR 770
      X2=XX/2.D0                                                        CDTR 780
      DLX2=DLOG(X2)                                                     CDTR 790
      GG=DBLE(G)                                                        CDTR 800
      G2=GG/2.D0                                                        CDTR 810
C                                                                       CDTR 820
C        COMPUTE ORDINATE                                               CDTR 830
C                                                                       CDTR 840
      CALL DLGAM(G2,GLG2,IOK)                                           CDTR 850
      DD=(G2-1.D0)*DLXX-X2-G2*.6931471805599453 -GLG2                   CDTR 860
      IF(DD-1.68D02) 110,110,120                                        CDTR 870
  110 IF(DD+1.68D02) 130,130,140                                        CDTR 880
  120 D=1.7E38                                                           CDTR 890
      GO TO 150                                                         CDTR 900
  130 D=0.0                                                             CDTR 910
      GO TO 150                                                         CDTR 920
  140 DD=DEXP(DD)                                                       CDTR 930
      D=SNGL(DD)                                                        CDTR 940
C                                                                       CDTR 950
C        TEST FOR G GREATER THAN 1000.0                                 CDTR 960
C        TEST FOR X GREATER THAN 2000.0                                 CDTR 970
C                                                                       CDTR 980
  150 IF(G-1000.) 160,160,180                                           CDTR 990
  160 IF(X-2000.) 190,190,170                                           CDTR1000
  170 P=1.0                                                             CDTR1010
      GO TO 610                                                         CDTR1020
  180 A=DLOG(XX/GG)/3.D0                                                CDTR1030
      A=DEXP(A)                                                         CDTR1040
      B=2.D0/(9.D0*GG)                                                  CDTR1050
      C=(A-1.D0+B)/DSQRT(B)                                             CDTR1060
      SC=SNGL(C)                                                        CDTR1070
      CALL NDTR(SC,P,DUMMY)                                             CDTR1080
      GO TO 490                                                         CDTR1090
C                                                                       CDTR1100
C        COMPUTE THETA                                                  CDTR1110
C                                                                       CDTR1120
  190 K= IDINT(G2)                                                      CDTR1130
      THETA=G2-DFLOAT(K)                                                CDTR1140
      IF(THETA-1.D-8) 200,200,210                                       CDTR1150
  200 THETA=0.D0                                                        CDTR1160
  210 THP1=THETA+1.D0                                                   CDTR1170
C                                                                       CDTR1180
C        SELECT METHOD OF COMPUTING T1                                  CDTR1190
C                                                                       CDTR1200
      IF(THETA) 230,230,220                                             CDTR1210
  220 IF(XX-10.D0) 260,260,320                                          CDTR1220
C                                                                       CDTR1230
C        COMPUTE T1 FOR THETA EQUALS 0.0                                CDTR1240
C                                                                       CDTR1250
  230 IF(X2-1.68D02) 250,240,240                                        CDTR1260
  240 T1=1.0                                                            CDTR1270
      GO TO 400                                                         CDTR1280
  250 T11=1.D0-DEXP(-X2)                                                CDTR1290
      T1=SNGL(T11)                                                      CDTR1300
      GO TO 400                                                         CDTR1310
C                                                                       CDTR1320
C        COMPUTE T1 FOR THETA GREATER THAN 0.0 AND                      CDTR1330
C        X LESS THAN OR EQUAL TO 10.0                                   CDTR1340
C                                                                       CDTR1350
  260 SER=X2*(1.D0/THP1 -X2/(THP1+1.D0))                                CDTR1360
      J=+1                                                              CDTR1370
      CC=DFLOAT(J)                                                      CDTR1380
      DO 270 IT1=3,30                                                   CDTR1390
      XI=DFLOAT(IT1)                                                    CDTR1400
      CALL DLGAM(XI,FAC,IOK)                                            CDTR1410
      TLOG= XI*DLX2-FAC-DLOG(XI+THETA)                                  CDTR1420
      TERM=DEXP(TLOG)                                                   CDTR1430
      TERM=DSIGN(TERM,CC)                                               CDTR1440
      SER=SER+TERM                                                      CDTR1450
      CC=-CC                                                            CDTR1460
      IF(DABS(TERM)-1.D-9) 280,270,270                                  CDTR1470
  270 CONTINUE                                                          CDTR1480
      GO TO 600                                                         CDTR1490
  280 IF(SER) 600,600,290                                               CDTR1500
  290 CALL DLGAM(THP1,GTH,IOK)                                          CDTR1510
      TLOG=THETA*DLX2+DLOG(SER)-GTH                                     CDTR1520
      IF(TLOG+1.68D02) 300,300,310                                      CDTR1530
  300 T1=0.0                                                            CDTR1540
      GO TO 400                                                         CDTR1550
  310 T11=DEXP(TLOG)                                                    CDTR1560
      T1=SNGL(T11)                                                      CDTR1570
      GO TO 400                                                         CDTR1580
C                                                                       CDTR1590
C        COMPUTE T1 FOR THETA GREATER THAN 0.0 AND                      CDTR1600
C        X GREATER THAN 10.0 AND LESS THAN 2000.0                       CDTR1610
C                                                                       CDTR1620
  320 A2=0.D0                                                           CDTR1630
      DO 340 I=1,25                                                     CDTR1640
      XI=DFLOAT(I)                                                      CDTR1650
      CALL DLGAM(THP1,GTH,IOK)                                          CDTR1660
      T11=-(13.D0*XX)/XI +THP1*DLOG(13.D0*XX/XI) -GTH-DLOG(XI)          CDTR1670
      IF(T11+1.68D02) 340,340,330                                       CDTR1680
  330 T11=DEXP(T11)                                                     CDTR1690
      A2=A2+T11                                                         CDTR1700
  340 CONTINUE                                                          CDTR1710
      A=1.01282051+THETA/156.D0-XX/312.D0                               CDTR1720
      B=DABS(A)                                                         CDTR1730
      C= -X2+THP1*DLX2+DLOG(B)-GTH-3.951243718581427                    CDTR1740
      IF(C+1.68D02) 370,370,350                                         CDTR1750
  350 IF (A) 360,370,380                                                CDTR1760
  360 C=-DEXP(C)                                                        CDTR1770
      GO TO 390                                                         CDTR1780
  370 C=0.D0                                                            CDTR1790
      GO TO 390                                                         CDTR1800
  380 C=DEXP(C)                                                         CDTR1810
  390 C=A2+C                                                            CDTR1820
      T11=1.D0-C                                                        CDTR1830
      T1=SNGL(T11)                                                      CDTR1840
C                                                                       CDTR1850
C        SELECT PROPER EXPRESSION FOR P                                 CDTR1860
C                                                                       CDTR1870
  400 IF(G-2.) 420,410,410                                              CDTR1880
  410 IF(G-4.) 450,460,460                                              CDTR1890
C                                                                       CDTR1900
C        COMPUTE P FOR G GREATER THAN ZERO AND LESS THAN 2.0            CDTR1910
C                                                                       CDTR1920
  420 CALL DLGAM(THP1,GTH,IOK)                                          CDTR1930
      DT2=THETA*DLXX-X2-THP1*.6931471805599453 -GTH                     CDTR1940
      IF(DT2+1.68D02) 430,430,440                                       CDTR1950
  430 P=T1                                                              CDTR1960
      GO TO 490                                                         CDTR1970
  440 DT2=DEXP(DT2)                                                     CDTR1980
      T2=SNGL(DT2)                                                      CDTR1990
      P=T1+T2+T2                                                        CDTR2000
      GO TO 490                                                         CDTR2010
C                                                                       CDTR2020
C        COMPUTE P FOR G GREATER THAN OR EQUAL TO 2.0                   CDTR2030
C        AND LESS THAN 4.0                                              CDTR2040
C                                                                       CDTR2050
  450 P=T1                                                              CDTR2060
      GO TO 490                                                         CDTR2070
C                                                                       CDTR2080
C        COMPUTE P FOR G GREATER THAN OR EQUAL TO 4.0                   CDTR2090
C        AND LESS THAN OR EQUAL TO 1000.0                               CDTR2100
C                                                                       CDTR2110
  460 DT3=0.D0                                                          CDTR2120
      DO 480 I3=2,K                                                     CDTR2130
      THPI=DFLOAT(I3)+THETA                                             CDTR2140
      CALL DLGAM(THPI,GTH,IOK)                                          CDTR2150
      DLT3=THPI*DLX2-DLXX-X2-GTH                                        CDTR2160
      IF(DLT3+1.68D02) 480,480,470                                      CDTR2170
  470 DT3=DT3+DEXP(DLT3)                                                CDTR2180
  480 CONTINUE                                                          CDTR2190
      T3=SNGL(DT3)                                                      CDTR2200
      P=T1-T3-T3                                                        CDTR2210
C                                                                       CDTR2220
C        SET ERROR INDICATOR                                            CDTR2230
C                                                                       CDTR2240
  490 IF(P) 500,520,520                                                 CDTR2250
  500 IF(ABS(P)-1.E-7) 510,510,600                                      CDTR2260
  510 P=0.0                                                             CDTR2270
      GO TO 610                                                         CDTR2280
  520 IF(1.-P) 530,550,550                                              CDTR2290
  530 IF(ABS(1.-P)-1.E-7) 540,540,600                                   CDTR2300
  540 P=1.0                                                             CDTR2310
      GO TO 610                                                         CDTR2320
  550 IF(P-1.E-8) 560,560,570                                           CDTR2330
  560 P=0.0                                                             CDTR2340
      GO TO 610                                                         CDTR2350
  570 IF((1.0-P)-1.E-8) 580,580,610                                     CDTR2360
  580 P=1.0                                                             CDTR2370
      GO TO 610                                                         CDTR2380
  590 IER=-1                                                            CDTR2390
      D=-1.7E38                                                          CDTR2400
      P=-1.7E38                                                          CDTR2410
      GO TO 620                                                         CDTR2420
  600 IER=+1                                                            CDTR2430
      P= 1.7E38                                                          CDTR2440
      GO TO 620                                                         CDTR2450
  610 IER=0                                                             CDTR2460
  620 RETURN                                                            CDTR2470
      END                                                               CDTR2480