      DOUBLE PRECISION FUNCTION DGAM1 (X)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 1998-10-29 DGAM1 Krogh  Moved external statement up for mangle.
c>> 1994-10-20 DGAM1 Krogh  Changes to use M77CON
c>> 1994-05-18 DGAM1 WVS make DP and SP versions alike using CHGTYP
C>> 1993-08-04 DGAM1 CLL Add type stmt for DGAMMA.
C>> 1993-05-06 DGAM1 WVS JPL Convert from NSWC to Math 77
c--D replaces "?": ?GAM1, ?GAMMA
C ----------------------------------------------------------------------
C     EVALUATION OF 1/GAMMA(1 + X) - 1  FOR -0.5 .LE. X .LE. 1.5
C ----------------------------------------------------------------------
C
C     THE FOLLOWING ARE THE FIRST 49 COEFFICIENTS OF THE MACLAURIN
C     EXPANSION FOR 1/GAMMA(1 + X) - 1. THE COEFFICIENTS ARE
C     CORRECT TO 40 DIGITS. THE COEFFICIENTS WERE OBTAINED BY
C     ALFRED H. MORRIS JR. (NAVAL SURFACE WARFARE CENTER) AND ARE
C     GIVEN HERE FOR REFERENCE. ONLY THE FIRST 14 COEFFICIENTS ARE
C     USED IN THIS CODE.
C
C                           -----------
C
C     DATA A(1)  / .5772156649015328606065120900824024310422D+00/,
C    *     A(2)  /-.6558780715202538810770195151453904812798D+00/,
C    *     A(3)  /-.4200263503409523552900393487542981871139D-01/,
C    *     A(4)  / .1665386113822914895017007951021052357178D+00/,
C    *     A(5)  /-.4219773455554433674820830128918739130165D-01/,
C    *     A(6)  /-.9621971527876973562114921672348198975363D-02/,
C    *     A(7)  / .7218943246663099542395010340446572709905D-02/,
C    *     A(8)  /-.1165167591859065112113971084018388666809D-02/,
C    *     A(9)  /-.2152416741149509728157299630536478064782D-03/,
C    *     A(10) / .1280502823881161861531986263281643233949D-03/
C     DATA A(11) /-.2013485478078823865568939142102181838229D-04/,
C    *     A(12) /-.1250493482142670657345359473833092242323D-05/,
C    *     A(13) / .1133027231981695882374129620330744943324D-05/,
C    *     A(14) /-.2056338416977607103450154130020572836513D-06/,
C    *     A(15) / .6116095104481415817862498682855342867276D-08/,
C    *     A(16) / .5002007644469222930055665048059991303045D-08/,
C    *     A(17) /-.1181274570487020144588126565436505577739D-08/,
C    *     A(18) / .1043426711691100510491540332312250191401D-09/,
C    *     A(19) / .7782263439905071254049937311360777226068D-11/,
C    *     A(20) /-.3696805618642205708187815878085766236571D-11/
C     DATA A(21) / .5100370287454475979015481322863231802727D-12/,
C    *     A(22) /-.2058326053566506783222429544855237419746D-13/,
C    *     A(23) /-.5348122539423017982370017318727939948990D-14/,
C    *     A(24) / .1226778628238260790158893846622422428165D-14/,
C    *     A(25) /-.1181259301697458769513764586842297831212D-15/,
C    *     A(26) / .1186692254751600332579777242928674071088D-17/,
C    *     A(27) / .1412380655318031781555803947566709037086D-17/,
C    *     A(28) /-.2298745684435370206592478580633699260285D-18/,
C    *     A(29) / .1714406321927337433383963370267257066813D-19/,
C    *     A(30) / .1337351730493693114864781395122268022875D-21/
C     DATA A(31) /-.2054233551766672789325025351355733796682D-21/,
C    *     A(32) / .2736030048607999844831509904330982014865D-22/,
C    *     A(33) /-.1732356445910516639057428451564779799070D-23/,
C    *     A(34) /-.2360619024499287287343450735427531007926D-25/,
C    *     A(35) / .1864982941717294430718413161878666898946D-25/,
C    *     A(36) /-.2218095624207197204399716913626860379732D-26/,
C    *     A(37) / .1297781974947993668824414486330594165619D-27/,
C    *     A(38) / .1180697474966528406222745415509971518560D-29/,
C    *     A(39) /-.1124584349277088090293654674261439512119D-29/,
C    *     A(40) / .1277085175140866203990206677751124647749D-30/
C     DATA A(41) /-.7391451169615140823461289330108552823711D-32/,
C    *     A(42) / .1134750257554215760954165259469306393009D-34/,
C    *     A(43) / .4639134641058722029944804907952228463058D-34/,
C    *     A(44) /-.5347336818439198875077418196709893320905D-35/,
C    *     A(45) / .3207995923613352622861237279082794391090D-36/,
C    *     A(46) /-.4445829736550756882101590352124643637401D-38/,
C    *     A(47) /-.1311174518881988712901058494389922190237D-38/,
C    *     A(48) / .1647033352543813886818259327906394145400D-39/,
C    *     A(49) /-.1056233178503581218600561071538285049997D-40/
C
C                           -----------
C
C     C = A(1) - 1 IS ALSO FREQUENTLY NEEDED. C HAS THE VALUE ...
C
C     DATA C /-.4227843350984671393934879099175975689578D+00/
C
C ----------------------------------------------------------------------
      EXTERNAL D1MACH, DGAMMA
      DOUBLE PRECISION D1MACH, DGAMMA
      DOUBLE PRECISION D, ROUND, T, W, X, XMAX, Z
      SAVE ROUND, XMAX
      DOUBLE PRECISION A0, A1, B1, B2, B3, B4, B5, B6, B7, B8
      DOUBLE PRECISION C, C0, C1, C2, C3, C4, C5, C6, C7, C8, C9,
     *                 C10, C11, C12, C13
      DOUBLE PRECISION F0, F1, F2, F3, F4, F5, F6, G1, G2, G3, G4
      DOUBLE PRECISION P1, P2, P3, P4, P5, P6, P7
      DOUBLE PRECISION Q2, Q3, Q4, Q5
      DOUBLE PRECISION R1, R2, R3, R4, R5, R6, R7, R8, R9
      DOUBLE PRECISION S1, S2
C ---------------------------
      DATA XMAX /-1.0D0/
C ---------------------------
      DATA A0 / .611609510448141581788D-08/,
     *     A1 / .624730830116465516210D-08/
      DATA B1 / .203610414066806987300D+00/,
     *     B2 / .266205348428949217746D-01/,
     *     B3 / .493944979382446875238D-03/,
     *     B4 /-.851419432440314906588D-05/,
     *     B5 /-.643045481779353022248D-05/,
     *     B6 / .992641840672773722196D-06/,
     *     B7 /-.607761895722825260739D-07/,
     *     B8 / .195755836614639731882D-09/
C ---------------------------
      DATA F0 /.6116095104481415817861D-08/,
     *     F1 /.6871674113067198736152D-08/,
     *     F2 /.6820161668496170657918D-09/,
     *     F3 /.4686843322948848031080D-10/,
     *     F4 /.1572833027710446286995D-11/,
     *     F5/-.1249441572276366213222D-12/,
     *     F6 /.4343529937408594255178D-14/
      DATA G1 /.3056961078365221025009D+00/,
     *     G2 /.5464213086042296536016D-01/,
     *     G3 /.4956830093825887312020D-02/,
     *     G4 /.2692369466186361192876D-03/
C ---------------------------
C     C = C0 - 1
C ---------------------------
      DATA C /-.422784335098467139393487909917598D+00/
C ---------------------------
      DATA C0  / .577215664901532860606512090082402D+00/,
     *     C1  /-.655878071520253881077019515145390D+00/,
     *     C2  /-.420026350340952355290039348754298D-01/,
     *     C3  / .166538611382291489501700795102105D+00/,
     *     C4  /-.421977345555443367482083012891874D-01/,
     *     C5  /-.962197152787697356211492167234820D-02/,
     *     C6  / .721894324666309954239501034044657D-02/,
     *     C7  /-.116516759185906511211397108401839D-02/,
     *     C8  /-.215241674114950972815729963053648D-03/,
     *     C9  / .128050282388116186153198626328164D-03/
      DATA C10 /-.201348547807882386556893914210218D-04/,
     *     C11 /-.125049348214267065734535947383309D-05/,
     *     C12 / .113302723198169588237412962033074D-05/,
     *     C13 /-.205633841697760710345015413002057D-06/
C --------------------------
      DATA P1/ .577215664901533D+00/, P2/-.409078193005776D+00/,
     *     P3/-.230975380857675D+00/, P4/ .597275330452234D-01/,
     *     P5/ .766968181649490D-02/, P6/-.514889771323592D-02/,
     *     P7/ .589597428611429D-03/
C --------------------------
      DATA Q2/ .427569613095214D+00/,
     *     Q3/ .158451672430138D+00/, Q4/ .261132021441447D-01/,
     *     Q5/ .423244297896961D-02/
C --------------------------
      DATA R1/-.422784335098468D+00/, R2/-.771330383816272D+00/,
     *     R3/-.244757765222226D+00/, R4/ .118378989872749D+00/,
     *     R5/ .930357293360349D-03/, R6/-.118290993445146D-01/,
     *     R7/ .223047661158249D-02/, R8/ .266505979058923D-03/,
     *     R9/-.132674909766242D-03/
C --------------------------
      DATA S1/ .273076135303957D+00/, S2/ .559398236957378D-01/
C ---------------------------
      IF (XMAX .LT. 0.0D0) THEN
         ROUND = D1MACH(4)
         XMAX = 1.0D0 / ROUND
      END IF
      T = X
      D = X - 0.5D0
      IF (D .GT. 0.D0) T = D - 0.5D0
      IF (ROUND .LT. 1.0D-14) THEN
C
C        CASE WHEN PRECISION IS GREATER THAN 14 DIGITS:
C
         IF (T) 30,10,20
C
   10    DGAM1 = 0.D0
         RETURN
C -----------
C
C                CASE WHEN 0 .LT. T .LE. 0.5
C
C              W IS A MINIMAX APPROXIMATION FOR
C              THE SERIES A(15) + A(16)*T + ...
C
C -----------
   20    IF (T .GT. 0.5) THEN
            IF (X .GT. XMAX) THEN
               DGAM1 = -1.0D0
               RETURN
            END IF
            DGAM1 = 1.0D0 / DGAMMA(1.0D0 + X) - 1.0D0
            RETURN
         END IF
         W = ((((((F6*T + F5)*T + F4)*T + F3)*T + F2)*T + F1)*T + F0)/
     *         ((((G4*T + G3)*T + G2)*T + G1)*T + 1.D0)
         Z = (((((((((((((W*T + C13)*T + C12)*T + C11)*T + C10)*T +
     *            C9)*T + C8)*T + C7)*T + C6)*T + C5)*T + C4)*T +
     *            C3)*T + C2)*T + C1)*T + C0
C
         IF (D .LE. 0.D0) THEN
            DGAM1 = X*Z
            RETURN
         END IF
         DGAM1 = (T/X)*((Z - 0.5D0) - 0.5D0)
         RETURN
C -----------
C
C                CASE WHEN -0.5 .LE. T .LT. 0
C
C              W IS A MINIMAX APPROXIMATION FOR
C              THE SERIES A(15) + A(16)*T + ...
C
C -----------
   30    IF (T .LT. -0.5) THEN
            IF (X .GT. XMAX) THEN
               DGAM1 = -1.0D0
               RETURN
            END IF
            DGAM1 = 1.0D0 / DGAMMA(1.0D0 + X) - 1.0D0
            RETURN
         END IF
         W = (A1*T + A0)/((((((((B8*T + B7)*T + B6)*T + B5)*T +
     *          B4)*T + B3)*T + B2)*T + B1)*T + 1.D0)
         Z = (((((((((((((W*T + C13)*T + C12)*T + C11)*T + C10)*T +
     *            C9)*T + C8)*T + C7)*T + C6)*T + C5)*T + C4)*T +
     *            C3)*T + C2)*T + C1)*T + C
C
         IF (D .LE. 0.D0) THEN
            DGAM1 = X*((Z + 0.5D0) + 0.5D0)
            RETURN
         END IF
         DGAM1 = T*Z/X
         RETURN
      ELSE
C
C        CASE WHEN PRECISION IS LESS THAN 14 DIGITS
C
         IF (T) 300,100,200
C
  100    DGAM1 = 0.0
         RETURN
C -----------
C
C                CASE WHEN 0 .LT. T .LE. 0.5
C
C              W IS A MINIMAX APPROXIMATION FOR
C              THE SERIES A(15) + A(16)*T + ...
C
C -----------
  200    IF (T .GT. 0.5D0) THEN
            IF (X .GT. XMAX) THEN
               DGAM1 = -1.0D0
               RETURN
            END IF
            DGAM1 = 1.0D0 / DGAMMA(1.0D0 + X) - 1.0D0
            RETURN
         END IF
         Z = ((((((P7*T + P6)*T + P5)*T + P4)*T + P3)*T + P2)*T + P1) /
     *       ((((Q5*T + Q4)*T + Q3)*T + Q2)*T + 1.0)
         IF (D .LE. 0.0) THEN
            DGAM1 = X*Z
            RETURN
         END IF
         DGAM1 = (T/X)*((Z - 0.5) - 0.5)
         RETURN
C
C -----------
C
C                CASE WHEN -0.5 .LE. T .LT. 0
C
C              W IS A MINIMAX APPROXIMATION FOR
C              THE SERIES A(15) + A(16)*T + ...
C
C -----------
  300    IF (T .LT. -0.5D0) THEN
            IF (X .GT. XMAX) THEN
               DGAM1 = -1.0D0
               RETURN
            END IF
            DGAM1 = 1.0D0 / DGAMMA(1.0D0 + X) - 1.0D0
            RETURN
         END IF
         Z = ((((((((R9*T + R8)*T + R7)*T + R6)*T + R5)*T
     *                       + R4)*T + R3)*T + R2)*T + R1) /
     *       ((S2*T + S1)*T + 1.0)
         IF (D .LE. 0.0) THEN
            DGAM1 = X*((Z + 0.5) + 0.5)
            RETURN
         END IF
         DGAM1 = T*Z/X
         RETURN
      END IF
      END
