C             ***** BENCHMARK  SP1111 *****
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER NN=5000
      COMMON/MISC/NGANGB
      COMMON/GEOM/RAB,RABSQ
      COMMON/PGEOM/GP(36),EP(36),DP00P(36),DP01P(36),DP10P(36),APP(36),
     $ BPP(36)
      COMMON/GINF/GCD
      COMMON/EABECD/EAB,ECD
      COMMON/QGEOM/AQZ,QPERP,QPERP2
      COMMON/CONST/CONP(36)
      COMMON/MAXC/ISMLP(36),ISMLQ,ISML
      COMMON/H/H(160)
      COMMON/PICON/PITO52,PIDIV4
      COMMON/IO/IN,IOUT,IPUNCH
C
      DATA TWO/2.0D0/,TWOPT5/2.5D0/,FOUR/4.0D0/,TEN/10.0D0/,F326/3.26077
     $D0/,F409/4.091326D0/,ONE/1.0D0/,PT2/0.2D0/,PT1/0.1D0/
      DATA ZERO/0.0D0/
      DATA F160/160.0D0/
      DATA CHEK/4.590475501D0/,THR/1.0D-08/
C
 2001 FORMAT(' FOR I = ',I6, '     SUM = ',D20.13)
 2002 FORMAT(' THERE WERE A TOTAL OF ',I6,' FAILURES.')
C
      IN=5
      IOUT=6
      IPUNCH=7
C
C
      CALL HEADER('BENCHMARK SP1111',16)
      CALL TIMRB
C
C
      PI=FOUR*DATAN(ONE)
      PITO52=TWO*(PI**TWOPT5)
      PIDIV4=PI/FOUR
C
      NGANGB=36
C
      RABSQ=PI/TEN
      RAB=DSQRT(RABSQ)
C
      EXXA=F326
      EXXB=F409
      GCD=EXXA+EXXB
C
      ECD=ONE/GCD
C
      DO 10 I=1,36
      GP(I)=TWO*F326
      EP(I)=ONE/GP(I)
      DP00P(I)=F326/TEN
      DP01P(I)=F326/TEN
      DP10P(I)=F326/TEN
      APP(I)=GP(I)
      BPP(I)=EP(I)
      CONP(I)=PITO52*GP(I)
   10 ISMLP(I)=0
C
      AQZ=PT1*PI
      QPERP2=PT2*PI
      QPERP=DSQRT(QPERP2)
C
      ISMLQ=0
C
      CALL GAMGEN(0)
C
      IFAIL=0
      DO 40 I=1,NN
      DO 20 J=1,160
   20 H(J)=ZERO
      CALL SP1111(I)
      SUM=ZERO
      DO 30 J=1,160
   30 SUM=SUM+DABS(H(J))
      SUM=SUM/F160
      IF(DABS(SUM-CHEK).LE.THR)GO TO 40
      WRITE(IOUT,2001)I,SUM
      IFAIL=IFAIL+1
   40 CONTINUE
C
      CALL TIMRE
C
C
      WRITE(IOUT,2002)IFAIL
C
      STOP
      END
      SUBROUTINE SP1111(IQXX)                                               0010
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                    0020
C                                                                           0030
C     --------------------------                                            0040
C     GAUSSIAN 74 (QCPE VERSION)                                            0050
C     DECEMBER 1974                                                         0060
C     UNIVAC 1108/CDC 7600                                                  0070
C     --------------------------                                            0080
C                                                                           0090
C                                                                           0100
      COMMON/MISC/NGANGB
      COMMON/GEOM/RAB,RABSQ
      COMMON/PGEOM/GP(36),EP(36),DP00P(36),DP01P(36),DP10P(36)
     *,APP(36),BPP(36)                                                      0190
      COMMON/PQGEOM/AP,BP
      COMMON/GINF/GCD
      COMMON/EABECD/EAB,ECD                                                 0230
      COMMON/DPQ/DP00,DP01,DP10
      COMMON/QGEOM/AQZ,QPERP,QPERP2
      COMMON/CONST/CONP(36)
      COMMON/MAXC/ISMLP(36),ISMLQ,ISML
      COMMON/ASTORE/QQ,THETA,N                                              0290
      COMMON/TABLE/                                                         0300
     *AA(400),BA(400),CA(400),                                              0310
     *AB(400),BB(400),CB(400),                                              0320
     *AC(400),BC(400),CC(400),                                              0330
     *AD(400),BD(400),CD(400),                                              0340
     *AE(400),BE(400),CE(400)                                               0350
      COMMON/AUXVAR/AUXVAR,VAR1,VAR2                                        0360
      COMMON/H/                                                             0370
     *H0000,H0001,H0002,H0003,H0011,H0012,H0013,H0022,H0023,H0033,          0380
     *H0100,H0101,H0102,H0103,H0111,H0112,H0113,H0122,H0123,H0133,          0390
     *H0200,H0201,H0202,H0203,H0211,H0212,H0213,H0222,H0223,H0233,          0400
     *H0300,H0301,H0302,H0303,H0311,H0312,H0313,H0322,H0323,H0333,          0410
     *H1000,H1001,H1002,H1003,H1011,H1012,H1013,H1022,H1023,H1033,          0420
     *H1100,H1101,H1102,H1103,H1111,H1112,H1113,H1122,H1123,H1133,          0430
     *H1200,H1201,H1202,H1203,H1211,H1212,H1213,H1222,H1223,H1233,          0440
     *H1300,H1301,H1302,H1303,H1311,H1312,H1313,H1322,H1323,H1333,          0450
     *H2000,H2001,H2002,H2003,H2011,H2012,H2013,H2022,H2023,H2033,          0460
     *H2100,H2101,H2102,H2103,H2111,H2112,H2113,H2122,H2123,H2133,          0470
     *H2200,H2201,H2202,H2203,H2211,H2212,H2213,H2222,H2223,H2233,          0480
     *H2300,H2301,H2302,H2303,H2311,H2312,H2313,H2322,H2323,H2333,          0490
     *H3000,H3001,H3002,H3003,H3011,H3012,H3013,H3022,H3023,H3033,          0500
     *H3100,H3101,H3102,H3103,H3111,H3112,H3113,H3122,H3123,H3133,          0510
     *H3200,H3201,H3202,H3203,H3211,H3212,H3213,H3222,H3223,H3233,          0520
     *H3300,H3301,H3302,H3303,H3311,H3312,H3313,H3322,H3323,H3333           0530
      COMMON/PICON/PITO52,PIDIV4
C                                                                           0540
      DATA ZERO/0.0D0/,PT25/0.25D0/,PT5/0.5D0/,ONE/1.0D0/,ONEPT5/1.5D0/,    0550
     *TWO/2.0D0/,TWOPT5/2.5D0/,THREE/3.0D0/,THRPT5/3.5D0/,TWENTY/20.0D0/    0560
C                                                                           0570
      X1=ZERO                                                               0580
      X2=ZERO                                                               0590
      X3=ZERO                                                               0600
      X4=ZERO                                                               0610
      X5=ZERO                                                               0620
      X6=ZERO                                                               0630
      Y1=ZERO                                                               0640
      Y2=ZERO                                                               0650
      Y3=ZERO                                                               0660
      Y4=ZERO                                                               0670
      Y5=ZERO                                                               0680
      Y6=ZERO                                                               0690
      Z1=ZERO                                                               0700
      Z2=ZERO                                                               0710
      Z3=ZERO                                                               0720
      Z4=ZERO                                                               0730
      Z5=ZERO                                                               0740
      Z6=ZERO                                                               0750
      Z7=ZERO                                                               0760
      Z8=ZERO                                                               0770
      Z9=ZERO                                                               0780
      V1=ZERO                                                               0790
      V2=ZERO                                                               0800
      V3=ZERO                                                               0810
      V4=ZERO                                                               0820
      V5=ZERO                                                               0830
      V6=ZERO                                                               0840
      W1=ZERO                                                               0850
      W2=ZERO                                                               0860
      W3=ZERO                                                               0870
      W4=ZERO                                                               0880
      W5=ZERO                                                               0890
      W6=ZERO                                                               0900
      W7=ZERO                                                               0910
      W8=ZERO                                                               0920
      W9=ZERO                                                               0930
      S1=ZERO                                                               0940
      S2=ZERO                                                               0950
      S3=ZERO                                                               0960
      S4=ZERO                                                               0970
      S5=ZERO                                                               0980
      S6=ZERO                                                               0990
      S7=ZERO                                                               1000
      S8=ZERO                                                               1010
      S9=ZERO                                                               1020
      S10=ZERO                                                              1030
      S11=ZERO                                                              1040
      S12=ZERO                                                              1050
      S13=ZERO                                                              1060
      S14=ZERO                                                              1070
      T1=ZERO                                                               1080
      T2=ZERO                                                               1090
      T3=ZERO                                                               1100
      T4=ZERO                                                               1110
      T5=ZERO                                                               1120
      T6=ZERO                                                               1130
      T7=ZERO                                                               1140
      T8=ZERO                                                               1150
      T9=ZERO                                                               1160
      T10=ZERO                                                              1170
      T11=ZERO                                                              1180
      T12=ZERO                                                              1190
      T13=ZERO                                                              1200
      T14=ZERO                                                              1210
      C1=ZERO                                                               1220
      C2=ZERO                                                               1230
      C3=ZERO                                                               1240
      C4=ZERO                                                               1250
      C5=ZERO                                                               1260
      C6=ZERO                                                               1270
      DO 90  IND=1,NGANGB                                                   1280
      ISML=ISMLQ+ISMLP(IND)
      IF(ISML-2)10,90,90                                                    1300
   10 IF(ISML-1)20,30,30
   20 AUXVAR=VAR1                                                           1320
      GO TO 40                                                              1330
   30 AUXVAR=VAR2                                                           1340
   40 EAB=EP(IND)                                                           1350
      DP00=DP00P(IND)                                                       1360
      DP01=DP01P(IND)                                                       1370
      DP10=DP10P(IND)                                                       1380
      AP=APP(IND)                                                           1390
      BP=BPP(IND)                                                           1400
      PQAB=AQZ-AP                                                           1410
      PQAB2=PQAB*PQAB                                                       1420
      G=ONE/(EAB+ECD)                                                       1430
      X=G*(QPERP2+PQAB2)                                                    1440
      IF(X-AUXVAR)60,60,50                                                  1450
   50 F0=CONP(IND)*DSQRT(PIDIV4/(X*(GP(IND)+GCD)))
      GTX=G/X                                                               1470
      F1=PT5*F0*GTX                                                         1480
      F2=ONEPT5*F1*GTX                                                      1490
      F3=TWOPT5*F2*GTX                                                      1500
      F4=THRPT5*F3*GTX                                                      1510
      GO TO 70                                                              1520
   60 Y=CONP(IND)/DSQRT(GP(IND)+GCD)                                        1530
      GY=G*Y                                                                1540
      GGY=G*GY                                                              1550
      GGGY=G*GGY                                                            1560
      QQ=X*TWENTY                                                           1570
      THETA=QQ-DINT(QQ)                                                     1580
      N=QQ-THETA                                                            1590
      THETA2=THETA*(THETA-ONE)                                              1600
      THETA3=THETA2*(THETA-TWO)                                             1610
      THETA4=THETA2*(THETA+ONE)                                             1620
      F0=(AA(N+1)+THETA*BA(N+1)-THETA3*CA(N+1)+THETA4*CA(N+2))*Y            1630
      F1=(AB(N+1)+THETA*BB(N+1)-THETA3*CB(N+1)+THETA4*CB(N+2))*GY           1640
      F2=(AC(N+1)+THETA*BC(N+1)-THETA3*CC(N+1)+THETA4*CC(N+2))*GGY          1650
      F3=(AD(N+1)+THETA*BD(N+1)-THETA3*CD(N+1)+THETA4*CD(N+2))*GGGY         1660
      F4=(AE(N+1)+THETA*BE(N+1)-THETA3*CE(N+1)+THETA4*CE(N+2))*GGGY*G       1670
   70 APBP=AP*BP                                                            1680
      EAB2=EAB*EAB                                                          1690
      BPDP01=BP*DP01                                                        1700
      APDP10=AP*DP10                                                        1710
      EDP01=EAB*DP01                                                        1720
      EDP10=EAB*DP10                                                        1730
      F1PQAB=F1*PQAB                                                        1740
      F2PQAB=F2*PQAB                                                        1750
      F3PQAB=F3*PQAB                                                        1760
      F4PQAB=F4*PQAB                                                        1770
      F1PQA2=F1*PQAB2                                                       1780
      F2PQA2=F2*PQAB2                                                       1790
      F3PQA2=F3*PQAB2                                                       1800
      F4PQA2=F4*PQAB2                                                       1810
      F2PQA3=F2PQA2*PQAB                                                    1820
      F3PQA3=F3PQA2*PQAB                                                    1830
      F4PQA3=F4PQA2*PQAB                                                    1840
      X1=X1+F0    *DP00                                                     1850
      X2=X2+F1    *DP00                                                     1860
      X3=X3+F2    *DP00                                                     1870
      X4=X4+F1PQAB*DP00                                                     1880
      X5=X5+F2PQAB*DP00                                                     1890
      X6=X6+F2PQA2*DP00                                                     1900
      Z1=Z1+F1          *EDP01                                              1910
      Z2=Z2+F2          *EDP01                                              1920
      Z3=Z3+F3          *EDP01                                              1930
      Z4=Z4+F1PQAB      *EDP01
      Z5=Z5+F2PQAB      *EDP01                                              1950
      Z6=Z6+F3PQAB      *EDP01                                              1960
      Z7=Z7+F2PQA2      *EDP01                                              1970
      Z8=Z8+F3PQA2      *EDP01                                              1980
      Z9=Z9+F3PQA3      *EDP01                                              1990
      W1=W1+F1          *EDP10                                              2000
      W2=W2+F2          *EDP10                                              2010
      W3=W3+F3          *EDP10                                              2020
      W4=W4+F1PQAB      *EDP10                                              2030
      W5=W5+F2PQAB      *EDP10                                              2040
      W6=W6+F3PQAB      *EDP10                                              2050
      W7=W7+F2PQA2      *EDP10                                              2060
      W8=W8+F3PQA2      *EDP10                                              2070
      W9=W9+F3PQA3      *EDP10                                              2080
      S1=S1+F0    *EAB                                                      2090
      S2=S2+F1    *EAB                                                      2100
      S3=S3+F2    *EAB                                                      2110
      S4=S4+F3    *EAB                                                      2120
      S6=S6+F1PQAB*EAB                                                      2130
      S7=S7+F2PQAB*EAB                                                      2140
      S8=S8+F3PQAB*EAB                                                      2150
      S9=S9+F1PQA2*EAB                                                      2160
      S10=S10+F2PQA2*EAB                                                    2170
      S11=S11+F3PQA2*EAB                                                    2180
      S12=S12+F2PQA3*EAB                                                    2190
      S13=S13+F3PQA3*EAB                                                    2200
      S14=S14+F3PQA3*PQAB*EAB                                               2210
      T1=T1+F0    *EAB2                                                     2220
      T2=T2+F1    *EAB2                                                     2230
      T3=T3+F2    *EAB2                                                     2240
      T4=T4+F3    *EAB2                                                     2250
      T5=T5+F4    *EAB2                                                     2260
      T6=T6+F2PQAB*EAB2                                                     2270
      T7=T7+F3PQAB*EAB2                                                     2280
      T8=T8+F4PQAB*EAB2                                                     2290
      T9=T9+F2PQA2*EAB2                                                     2300
      T10=T10+F3PQA2*EAB2                                                   2310
      T11=T11+F4PQA2*EAB2                                                   2320
      T12=T12+F3PQA3*EAB2                                                   2330
      T13=T13+F4PQA3*EAB2                                                   2340
      T14=T14+F4PQA3*PQAB*EAB2                                              2350
      IF (RABSQ) 80,90,80                                                   2360
   80 CONTINUE                                                              2370
      Y1=Y1+F0    *BPDP01                                                   2380
      Y2=Y2+F1    *BPDP01                                                   2390
      Y3=Y3+F2    *BPDP01                                                   2400
      Y4=Y4+F1PQAB*BPDP01                                                   2410
      Y5=Y5+F2PQAB*BPDP01                                                   2420
      Y6=Y6+F2PQA2*BPDP01                                                   2430
      V1=V1+F0    *APDP10                                                   2440
      V2=V2+F1    *APDP10                                                   2450
      V3=V3+F2    *APDP10                                                   2460
      V4=V4+F1PQAB*APDP10                                                   2470
      V5=V5+F2PQAB*APDP10                                                   2480
      V6=V6+F2PQA2*APDP10                                                   2490
      C1=C1+F0    *APBP                                                     2500
      C2=C2+F1    *APBP                                                     2510
      C3=C3+F2    *APBP                                                     2520
      C4=C4+F1PQAB*APBP                                                     2530
      C5=C5+F2PQAB*APBP                                                     2540
      C6=C6+F2PQA2*APBP                                                     2550
   90 CONTINUE                                                              2560
      A1=AQZ*S2-S6                                                          2570
      A2=AQZ*S3-S7                                                          2580
      A3=AQZ*S4-S8                                                          2590
      A4=AQZ*S6-S9                                                          2600
      A5=AQZ*S7-S10                                                         2610
      A6=AQZ*S8-S11                                                         2620
      A8=AQZ*S10-S12                                                        2630
      A9=AQZ*S11-S13                                                        2640
      A10=AQZ*S13-S14                                                       2650
      BQZ=AQZ-RAB                                                           2660
      B1=BQZ*S2-S6                                                          2670
      B2=BQZ*S3-S7                                                          2680
      B3=BQZ*S4-S8                                                          2690
      B4=BQZ*S6-S9                                                          2700
      B5=BQZ*S7-S10                                                         2710
      B6=BQZ*S8-S11                                                         2720
      B8=BQZ*S10-S12                                                        2730
      B9=BQZ*S11-S13                                                        2740
      B10=BQZ*S13-S14                                                       2750
      HECD=PT5*ECD                                                          2760
      ECD2=ECD*ECD                                                          2770
      HECD2=PT5*ECD2                                                        2780
      QECD=QPERP*ECD                                                        2790
      HQECD=PT5*QECD                                                        2800
      QECD2=QPERP*ECD2                                                      2810
      HQECD2=PT5*QECD2                                                      2820
      Q2ECD=QPERP2*ECD                                                      2830
      Q3ECD=QPERP*Q2ECD                                                     2840
      Q2ECD2=QPERP2*ECD2                                                    2850
      Q3ECD2=Q2ECD2*QPERP                                                   2860
      H0000=X1                                                              2870
      H0001=QECD*X2                                                         2880
      H0003=-ECD*X4                                                         2890
      H0022=HECD*(X1-ECD*X2)                                                2900
      H0011=H0022+Q2ECD2*X3                                                 2910
      H0013=-QECD2*X5                                                       2920
      H0033=H0022+ECD2*X6                                                   2930
      H0100=-QPERP*Z1                                                       2940
      H0300=Z4+Y1                                                           2950
      H0202=HECD*Z1                                                         2960
      H0101=H0202-Q2ECD*Z2                                                  2970
      H0103=QECD*Z5                                                         2980
      H0301=H0103+QECD*Y2                                                   2990
      H0303=H0202-ECD*Z7-ECD*Y4                                             3000
      H0212=HQECD2*Z2                                                       3010
      H0223=-HECD2*Z5                                                       3020
      H0122=H0212-QPERP*H0202                                               3030
      H0322=H0223+HECD*(H0300-ECD*Y2)                                       3040
      H0113=H0223+Q2ECD2*Z6                                                 3050
      H0313=H0212-QECD2*(Z8+Y5)                                             3060
      H0111=H0122+H0212+H0212-Q3ECD2*Z3                                     3070
      H0133=H0122-QECD2*Z8                                                  3080
      H0311=H0322+Q2ECD2*(Z6+Y3)                                            3090
      H0333=H0322+H0223+H0223+ECD2*(Z9+Y6)                                  3100
      H1000=-QPERP*W1                                                       3110
      H3000=W4+V1                                                           3120
      H2002=HECD*W1                                                         3130
      H1001=H2002-Q2ECD*W2                                                  3140
      H1003=QECD*W5                                                         3150
      H3001=H1003+QECD*V2                                                   3160
      H3003=H2002-ECD*W7-ECD*V4                                             3170
      H2012=HQECD2*W2                                                       3180
      H2023=-HECD2*W5                                                       3190
      H1022=H2012-QPERP*H2002                                               3200
      H3022=H2023+HECD*(H3000-ECD*V2)                                       3210
      H1013=H2023+Q2ECD2*W6                                                 3220
      H3013=H2012-QECD2*(W8+V5)                                             3230
      H1011=H1022+H2012+H2012-Q3ECD2*W3                                     3240
      H1033=H1022-QECD2*W8                                                  3250
      H3011=H3022+Q2ECD2*(W6+V3)                                            3260
      H3033=H3022+H2023+H2023+ECD2*(W9+V6)                                  3270
      H2200=PT5*(S1-T2)                                                     3280
      H1100=H2200+QPERP2*T3                                                 3290
      H1300=-QPERP*(T6+B1)                                                  3300
      H3100=-QPERP*(T6+A1)                                                  3310
      H3300=H2200+T9+A4+B4+C1                                               3320
      H2201=HQECD*(S2-T3)                                                   3330
      H1101=H2201-QECD*T3+Q3ECD*T4                                          3340
      TEMP =HECD*T6-Q2ECD*T7                                                3350
      H1301=TEMP+HECD*B1-Q2ECD*B2                                           3360
      H3101=TEMP+HECD*A1-Q2ECD*A2                                           3370
      H3301=H2201+QECD*(T10+A5+B5+C2)                                       3380
      H1202=-HQECD*T3                                                       3390
      H2102=H1202                                                           3400
      H2302=HECD*(T6+B1)                                                    3410
      H3202=HECD*(T6+A1)                                                    3420
      H2203=HECD*(T6-S6)                                                    3430
      H1103=H2203-Q2ECD*T7                                                  3440
      TEMP =-HQECD*T3+QECD*T10                                              3450
      H1303=TEMP+QECD*B5                                                    3460
      H3103=TEMP+QECD*A5                                                    3470
      H3303=H2203+ECD*(T6-T12-A8-B8-C4)+HECD*(A1+B1)                        3480
      H1212=PT25*ECD2*T3-PT5*Q2ECD2*T4                                      3490
      H2112=H1212                                                           3500
      H2312=HQECD2*(T7+B2)                                                  3510
      H3212=HQECD2*(T7+A2)                                                  3520
      H1223=HQECD2*T7                                                       3530
      H2123=H1223                                                           3540
      H2323=HECD2*(PT5*T3-T10-B5)                                           3550
      H3223=HECD2*(PT5*T3-T10-A5)                                           3560
      HXXYY=PT25*(ECD*(S1-T2)-ECD2*(S2-T3))                                 3570
      H2222=HXXYY+HECD2*T3                                                  3580
      H1122=HXXYY+PT5*(Q2ECD*T3-Q2ECD2*T4)                                  3590
      TEMP =HQECD*(ECD*T7-T6)                                               3600
      H1322=TEMP+HQECD*(ECD*B2-B1)                                          3610
      H3122=TEMP+HQECD*(ECD*A2-A1)                                          3620
      H3322=HXXYY+HECD*(T9+A4+B4+C1)-HECD2*(T10+A5+B5+C2)                   3630
      H2211=HXXYY+PT5*Q2ECD2*(S3-T4)                                        3640
      H1111=HXXYY+(HECD2+PT5*Q2ECD)*T3+Q2ECD2*(-THREE*T4+PT5*S3+QPERP2*T    3650
     15)                                                                    3660
      H1311=ONEPT5*QECD2*(T7+B2)-HQECD*(T6+B1)-Q3ECD2*(B3+T8)               3670
      H3111=ONEPT5*QECD2*(T7+A2)-HQECD*(T6+A1)-Q3ECD2*(A3+T8)               3680
      H3311=HXXYY-HECD2*(QPERP2*T4+T10+A5+B5)+HECD*(T9+A4+B4+C1-ECD*C2)+    3690
     1Q2ECD2*(T11+PT5*S3+A6+B6+C3)                                          3700
      H2213=HQECD2*(T7-S7)                                                  3710
      H1113=ONEPT5*QECD2*T7-HQECD2*S7-Q3ECD2*T8                             3720
      TEMP =HECD2*(PT5*T3-T10)+Q2ECD2*(T11-PT5*T4)                          3730
      H1313=TEMP-HECD2*B5+Q2ECD2*B6                                         3740
      H3113=TEMP-HECD2*A5+Q2ECD2*A6                                         3750
      H3313=QECD2*(ONEPT5*T7-T13-A9-B9-C5)-HQECD2*(S7-A2-B2)                3760
      H2233=HXXYY+HECD2*(S10-T10)                                           3770
      H1133=HXXYY-HECD2*(QPERP2*T4+T10-S10)+PT5*Q2ECD*T3+Q2ECD2*T11         3780
      H1333=QECD2*(ONEPT5*T7-T13-B9)-HQECD*(T6+B1)+HQECD2*B2                3790
      H3133=QECD2*(ONEPT5*T7-T13-A9)-HQECD*(T6+A1)+HQECD2*A2                3800
      H3333=HXXYY+HECD2*(-THREE*(A5+B5)+T3+S10-C2)+ECD2*(-THREE*T10+T14+    3810
     1A10+  B10+C6)+HECD*(T9+A4+B4+C1)                                      3820
      RETURN                                                                3830
      END                                                                   3840
      SUBROUTINE GAMGEN(KOP)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     --------------------------
C     GAUSSIAN 74 (QCPE VERSION)
C     DECEMBER 1974
C     UNIVAC 1108/CDC 7600
C     --------------------------
C
      COMMON/TABLE/C(1200,5)
      DIMENSION Y(410),F(9)
      DATA PT15/0.15D0/,PT05/0.05D0/
      DATA PT184/0.184D0/
      DATA SIX/6.0D0/,FOUR/4.0D0/
      DATA IOGAM/3/
C
C     INITIALIZE THE F(M,T) ROUTINE.
      CALL FMTSET(0,0,0)
C
C     GENERATE THE DESIRED F(M,T) FOR THE COMPLETE RANGE.
C     WE WILL COMPUTE F(0,T) THROUGH F(4,T) IN THIS SECTION.
      T=-PT15
      DO 10 I=1,404
      T=T+PT05
      CALL FMTGEN(F,T,5,ICK)
C     COPY THE RETURNED VALUES INTO PLACES WHERE THEY CAN BE
C     REACHED LATER.
      C(I,2)=F(1)
      C(I,3)=F(2)
      C(I,4)=F(3)
      C(I,5)=F(4)
   10 Y(I  )=F(5)
C
C     COMPUTE THE INTERPOLATION TABLE WITH THE VALUES AVAILABLE.
C     THIS IS SOMEWHAT COMPLICATED BY THE FACT THAT SOME OF THE
C     INTEGRALS ARE IN Y RATHER THAN IN C.
      DO 40 K=1,5
      DO 40 I=1,400
      J=I+2
C     K INDEXES THE VALUE OF M IN F(M,T).
C     I INDEXES THE INTERPOLATION TABLE (C).
C     J INDEXES THE STORED VALUES OF F(M,T).
C     IF K .EQ. 5, WE MUST USE ALTERNATE CODE BECAUSE THE INTEGRALS
C     ARE STORED IN Y.
      IF(K-5)20,30,20
C     INTEGRALS IN C, PROCEED AS NORMAL.
   20 TEMP1=C(J+1,K+1)+C(J-1,K+1)-(C(J,K+1)+C(J,K+1))
      TEMP2=SIX*C(J,K+1)-FOUR*(C(J+1,K+1)+C(J-1,K+1))+C(J-2,K+1)+C(J+2,K
     $+1)
      C(I    ,K)=C(J,K+1)
      C(I+400,K)=C(J+1,K+1)-C(J,K+1)
      C(I+800,K)=(TEMP1-PT184*TEMP2)/SIX
      GO TO 40
C     ALTERNATE CODE USING Y.
   30 TEMP1=Y(J+1)+Y(J-1)-(Y(J)+Y(J))
      TEMP2=SIX*Y(J)-FOUR*(Y(J+1)+Y(J-1))+Y(J-2)+Y(J+2)
      C(I    ,K)=Y(J)
      C(I+400,K)=Y(J+1)-Y(J)
      C(I+800,K)=(TEMP1-PT184*TEMP2)/SIX
   40 CONTINUE
      RETURN
      END
      SUBROUTINE FMTSET(KOP1,KOP2,KOP3)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     ROUTINE TO PRE-SET CUTOFFS FOR FMTGEN.
C
      COMMON/FM/GA(15),RPITWO,FMZERO(15),TOL,CUT0S,CUTSM,CUTML
      COMMON/IO/IN,IOUT,IPUNCH
      COMMON/FMCONS/FOUR,ONE,HALF,TWO,ZERO,TEN,TENM9,F20,F42,F500
      DIMENSION F(15)
      DATA FOUR/4.0D0/,ONE/1.0D0/,HALF/0.5D0/,TWO/2.0D0/,ZERO/0.0D0/
      DATA TEN/10.0D0/,TENM9/1.0D-9/,F20/20.0D0/,F42/42.0D0/,F500/500.0D
     $0/
C
 2001 FORMAT(21H FROM FMTSET, KOPS = ,3I2,12H AND CUTS = ,3E10.3)
C
C     COMPUTE PI RELATED CONSTANTS.
      PI=FOUR*DATAN(ONE)
      GA(1)=DSQRT(PI)
      RPITWO=HALF*GA(1)
C
C     FILL GAMMA FUNCTION ARRAY.
      TOL=HALF
      DO 10 I=2,15
      GA(I)=GA(I-1)*TOL
   10 TOL=TOL+ONE
C
C     FILL FMZERO (FOR ARGUMENT OF ZERO).
      TOL=ONE
      FMZERO(1)=ONE
      DO 20 I=2,15
      TOL=TOL+TWO
   20 FMZERO(I)=ONE/TOL
C
C     OBTAIN CUTOFFS.
C
C        CUT0S = CUTOFF, ZERO TO SMALL.
C        CUTSM = CUTOFF, SMALL TO MEDIUM.
C        CUTML = CUTOFF, MEDIUM TO LARGE.
C
C     CUTOS IS OBTAINED FROM KOP1.
      CUT0S=ZERO
      IF(KOP1)30,40,30
   30 CUT0S=TEN**(-2*KOP1)
C
C     CUTSM IS OBTAINED FROM KOP2.
C     KOP2 SPECIFIES THE DESIRED ACCURACY OF THE ASSYMPTOTIC
C     EXPANSION.  ROUTINE FMTGEN IS CALLED TO CHECK THE VIABILITY
C     OF THE ASSYMPTOTIC EXPANSION FOR VARIOUS VALUES OF CUTSM.
   40 TOL=TENM9
      CUTSM=TEN
      IF(KOP2)50,90,50
   50 TOL=TEN**(-6-KOP2)
      T=F20
   60 CALL FMTGEN(F,T,1,ICK)
      IF(ICK)80,70,80
   70 T=T-ONE
      IF(T-TEN)90,60,60
   80 CUT=T+ONE
C
C     CUTML IS TAKEN FROM KOP3.
   90 CUTML=F42
      IF(KOP3)100,130,100
  100 IF(KOP3-7)120,110,120
  110 CUTML=F500
      GO TO 130
  120 CUTML=DFLOAT(KOP3)*FOUR
  130 IF((KOP1+KOP2+KOP3).NE.0)WRITE(IOUT,2001)KOP1,KOP2,KOP3,CUT0S,
     $ CUTSM,CUTML
      RETURN
      END
      SUBROUTINE FMTGEN(F,T,M,ICK)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     --------------------------
C     GAUSSIAN 74 (QCPE VERSION)
C     DECEMBER 1974
C     UNIVAC 1108/CDC 7600
C     --------------------------
C
      DIMENSION F(M)
      COMMON/IO/IN,IOUT,IPUNCH
      COMMON/FMCONS/FOUR,ONE,HALF,TWO,ZERO,TEN,TENM9,F20,F42,F500
      COMMON/FM/GA(15),RPITWO,FMZERO(15),TOL,CUT0S,CUTSM,CUTML
      EQUIVALENCE(APPROX,OLDSUM)
      DATA ZERO/0.0D0/,TEN/10.0D0/,HALF/0.5D0/,ONE/1.0D0/,F42/42.0D0/
      DATA TWO/2.0D0/
 2001 FORMAT(41H0FAILURE IN FMGEN FOR SMALL T:  IX ' 50, /
     $ 6H IX = ,I3,7H,  T = ,E20.14)
 2002 FORMAT(37H0FAILURE IN FMGEN FOR INTERMEDIATE T:/
     $  6H  T = ,E20.14)
C
C     ICK IS AN ERROR INDICATOR.
C     ON RETURN, ICK=0 IMPLIES THAT ALL IS WELL.
C     IF ON RETURN, ICK IS NON-ZERO, THE ASSYMPTOTIC EXPANSION
C     HAS FAILED.
      ICK=0
C     TEST FOR TYPE OF ALGORITHM.
      IF(DABS(T)-CUT0S)1,1,3
C***********************************************************************
C        FILL F(M) FOR ARGUMENT OF ZERO.
C************************************************************************
    1 DO 2 I=1,M
    2 F(I)=FMZERO(I)
      RETURN
C     TEST FOR EVALUATION OF THE EXP.
    3 TEXP=ZERO
      IF(DABS(T)-CUTML)4,150,150
    4 TEXP=DEXP(-T)
      IF(DABS(T)-CUTSM)10,80,80
C***********************************************************************
C        0 .LT. T .LT. 10
C***********************************************************************
   10 A=DFLOAT(M-1)+HALF
      TERM=ONE/A
      SUM=TERM
      DO 20 IX=2,200
      A=A+ONE
      TERM=TERM*T/A
      SUM=SUM+TERM
      IF(DABS(TERM/SUM)-TOL)30,20,20
   20 CONTINUE
      WRITE(IOUT,2001)IX,T
      STOP
   30 F(M)=HALF*SUM*TEXP
      GO TO 160
C
C***********************************************************************
C        10 .LE. T .LT. 42
C***********************************************************************
   80 A=DFLOAT(M-1)
      B=A+HALF
      A=A-HALF
      TX=ONE/T
      MM1=M-1
      APPROX=RPITWO*DSQRT(TX)*(TX**MM1)
      IF(MM1)90,110,90
   90 DO 100 IX=1,MM1
      B=B-ONE
  100 APPROX=APPROX*B
  110 FIMULT=HALF*TEXP*TX
      SUM=ZERO
      IF(FIMULT)120,140,120
  120 FIPROP=FIMULT/APPROX
      TERM=ONE
      SUM =ONE
      NOTRMS=IDINT(T)+MM1
      DO 130 IX=2,NOTRMS
      TERM=TERM*A*TX
      SUM=SUM+TERM
      IF(DABS(TERM*FIPROP/SUM)-TOL)140,140,130
  130 A=A-ONE
      WRITE(IOUT,2002)T
      ICK=1
      RETURN
  140 F(M)=APPROX-FIMULT*SUM
      GO TO 160
C***********************************************************************
C        T .GE. 42
C***********************************************************************
  150 TX=DFLOAT(M)-HALF
      F(M)=HALF*GA(M)/(T**TX)
C***********************************************************************
C        RECUR DOWNWARDS TO F(1)
C***********************************************************************
  160 TX=T+T
      SUM=DFLOAT(M+M-3)
      MM1=M-1
      IF(MM1)170,190,170
  170 DO 180 IX=1,MM1
      F(M-IX)=(TX*F(M-IX+1)+TEXP)/SUM
  180 SUM=SUM-TWO
  190 RETURN
      END

      IF(MM1)170,190,170
  170 DO 180 IX=1,MM1
      F(M-IX)=(TX*F(M-IX+1)+TEXP)/SUM
  180 SUM=SUM-TWO
  190 RETURN
      END