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 .