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