/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:55 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dgam1.h" #include #include double /*FUNCTION*/ dgam1( double x) { long int _l0; double d, dgam1_v, t, w, z; static double round; static double xmax = -1.0e0; static double a0 = .611609510448141581788e-08; static double a1 = .624730830116465516210e-08; static double b1 = .203610414066806987300e00; static double b2 = .266205348428949217746e-01; static double b3 = .493944979382446875238e-03; static double b4 = -.851419432440314906588e-05; static double b5 = -.643045481779353022248e-05; static double b6 = .992641840672773722196e-06; static double b7 = -.607761895722825260739e-07; static double b8 = .195755836614639731882e-09; static double f0 = .6116095104481415817861e-08; static double f1 = .6871674113067198736152e-08; static double f2 = .6820161668496170657918e-09; static double f3 = .4686843322948848031080e-10; static double f4 = .1572833027710446286995e-11; static double f5 = -.1249441572276366213222e-12; static double f6 = .4343529937408594255178e-14; static double g1 = .3056961078365221025009e00; static double g2 = .5464213086042296536016e-01; static double g3 = .4956830093825887312020e-02; static double g4 = .2692369466186361192876e-03; static double c = -.422784335098467139393487909917598e00; static double c0 = .577215664901532860606512090082402e00; static double c1 = -.655878071520253881077019515145390e00; static double c2 = -.420026350340952355290039348754298e-01; static double c3 = .166538611382291489501700795102105e00; static double c4 = -.421977345555443367482083012891874e-01; static double c5 = -.962197152787697356211492167234820e-02; static double c6 = .721894324666309954239501034044657e-02; static double c7 = -.116516759185906511211397108401839e-02; static double c8 = -.215241674114950972815729963053648e-03; static double c9 = .128050282388116186153198626328164e-03; static double c10 = -.201348547807882386556893914210218e-04; static double c11 = -.125049348214267065734535947383309e-05; static double c12 = .113302723198169588237412962033074e-05; static double c13 = -.205633841697760710345015413002057e-06; static double p1 = .577215664901533e00; static double p2 = -.409078193005776e00; static double p3 = -.230975380857675e00; static double p4 = .597275330452234e-01; static double p5 = .766968181649490e-02; static double p6 = -.514889771323592e-02; static double p7 = .589597428611429e-03; static double q2 = .427569613095214e00; static double q3 = .158451672430138e00; static double q4 = .261132021441447e-01; static double q5 = .423244297896961e-02; static double r1 = -.422784335098468e00; static double r2 = -.771330383816272e00; static double r3 = -.244757765222226e00; static double r4 = .118378989872749e00; static double r5 = .930357293360349e-03; static double r6 = -.118290993445146e-01; static double r7 = .223047661158249e-02; static double r8 = .266505979058923e-03; static double r9 = -.132674909766242e-03; static double s1 = .273076135303957e00; static double s2 = .559398236957378e-01; /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1998-10-29 DGAM1 Krogh Moved external statement up for mangle. *>> 1994-10-20 DGAM1 Krogh Changes to use M77CON *>> 1994-05-18 DGAM1 WVS make DP and SP versions alike using CHGTYP *>> 1993-08-04 DGAM1 CLL Add type stmt for DGAMMA. *>> 1993-05-06 DGAM1 WVS JPL Convert from NSWC to Math 77 *--D replaces "?": ?GAM1, ?GAMMA * ---------------------------------------------------------------------- * EVALUATION OF 1/GAMMA(1 + X) - 1 FOR -0.5 .LE. X .LE. 1.5 * ---------------------------------------------------------------------- * * THE FOLLOWING ARE THE FIRST 49 COEFFICIENTS OF THE MACLAURIN * EXPANSION FOR 1/GAMMA(1 + X) - 1. THE COEFFICIENTS ARE * CORRECT TO 40 DIGITS. THE COEFFICIENTS WERE OBTAINED BY * ALFRED H. MORRIS JR. (NAVAL SURFACE WARFARE CENTER) AND ARE * GIVEN HERE FOR REFERENCE. ONLY THE FIRST 14 COEFFICIENTS ARE * USED IN THIS CODE. * * ----------- * * DATA A(1) / .5772156649015328606065120900824024310422D+00/, * * A(2) /-.6558780715202538810770195151453904812798D+00/, * * A(3) /-.4200263503409523552900393487542981871139D-01/, * * A(4) / .1665386113822914895017007951021052357178D+00/, * * A(5) /-.4219773455554433674820830128918739130165D-01/, * * A(6) /-.9621971527876973562114921672348198975363D-02/, * * A(7) / .7218943246663099542395010340446572709905D-02/, * * A(8) /-.1165167591859065112113971084018388666809D-02/, * * A(9) /-.2152416741149509728157299630536478064782D-03/, * * A(10) / .1280502823881161861531986263281643233949D-03/ * DATA A(11) /-.2013485478078823865568939142102181838229D-04/, * * A(12) /-.1250493482142670657345359473833092242323D-05/, * * A(13) / .1133027231981695882374129620330744943324D-05/, * * A(14) /-.2056338416977607103450154130020572836513D-06/, * * A(15) / .6116095104481415817862498682855342867276D-08/, * * A(16) / .5002007644469222930055665048059991303045D-08/, * * A(17) /-.1181274570487020144588126565436505577739D-08/, * * A(18) / .1043426711691100510491540332312250191401D-09/, * * A(19) / .7782263439905071254049937311360777226068D-11/, * * A(20) /-.3696805618642205708187815878085766236571D-11/ * DATA A(21) / .5100370287454475979015481322863231802727D-12/, * * A(22) /-.2058326053566506783222429544855237419746D-13/, * * A(23) /-.5348122539423017982370017318727939948990D-14/, * * A(24) / .1226778628238260790158893846622422428165D-14/, * * A(25) /-.1181259301697458769513764586842297831212D-15/, * * A(26) / .1186692254751600332579777242928674071088D-17/, * * A(27) / .1412380655318031781555803947566709037086D-17/, * * A(28) /-.2298745684435370206592478580633699260285D-18/, * * A(29) / .1714406321927337433383963370267257066813D-19/, * * A(30) / .1337351730493693114864781395122268022875D-21/ * DATA A(31) /-.2054233551766672789325025351355733796682D-21/, * * A(32) / .2736030048607999844831509904330982014865D-22/, * * A(33) /-.1732356445910516639057428451564779799070D-23/, * * A(34) /-.2360619024499287287343450735427531007926D-25/, * * A(35) / .1864982941717294430718413161878666898946D-25/, * * A(36) /-.2218095624207197204399716913626860379732D-26/, * * A(37) / .1297781974947993668824414486330594165619D-27/, * * A(38) / .1180697474966528406222745415509971518560D-29/, * * A(39) /-.1124584349277088090293654674261439512119D-29/, * * A(40) / .1277085175140866203990206677751124647749D-30/ * DATA A(41) /-.7391451169615140823461289330108552823711D-32/, * * A(42) / .1134750257554215760954165259469306393009D-34/, * * A(43) / .4639134641058722029944804907952228463058D-34/, * * A(44) /-.5347336818439198875077418196709893320905D-35/, * * A(45) / .3207995923613352622861237279082794391090D-36/, * * A(46) /-.4445829736550756882101590352124643637401D-38/, * * A(47) /-.1311174518881988712901058494389922190237D-38/, * * A(48) / .1647033352543813886818259327906394145400D-39/, * * A(49) /-.1056233178503581218600561071538285049997D-40/ * * ----------- * * C = A(1) - 1 IS ALSO FREQUENTLY NEEDED. C HAS THE VALUE ... * * DATA C /-.4227843350984671393934879099175975689578D+00/ * * ---------------------------------------------------------------------- */ /* --------------------------- */ /* --------------------------- */ /* --------------------------- */ /* --------------------------- * C = C0 - 1 * --------------------------- */ /* --------------------------- */ /* -------------------------- */ /* -------------------------- */ /* -------------------------- */ /* -------------------------- */ /* --------------------------- */ if (xmax < 0.0e0) { round = DBL_EPSILON; xmax = 1.0e0/round; } t = x; d = x - 0.5e0; if (d > 0.e0) t = d - 0.5e0; if (round < 1.0e-14) { /* CASE WHEN PRECISION IS GREATER THAN 14 DIGITS: * */ switch (ARITHIF(t)) { case -1: goto L_30; case 0: goto L_10; case 1: goto L_20; } L_10: dgam1_v = 0.e0; return( dgam1_v ); /* ----------- * * CASE WHEN 0 .LT. T .LE. 0.5 * * W IS A MINIMAX APPROXIMATION FOR * THE SERIES A(15) + A(16)*T + ... * * ----------- */ L_20: if (t > 0.5) { if (x > xmax) { dgam1_v = -1.0e0; return( dgam1_v ); } dgam1_v = 1.0e0/dgamma( 1.0e0 + x ) - 1.0e0; return( dgam1_v ); } 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; if (d <= 0.e0) { dgam1_v = x*z; return( dgam1_v ); } dgam1_v = (t/x)*((z - 0.5e0) - 0.5e0); return( dgam1_v ); /* ----------- * * CASE WHEN -0.5 .LE. T .LT. 0 * * W IS A MINIMAX APPROXIMATION FOR * THE SERIES A(15) + A(16)*T + ... * * ----------- */ L_30: if (t < -0.5) { if (x > xmax) { dgam1_v = -1.0e0; return( dgam1_v ); } dgam1_v = 1.0e0/dgamma( 1.0e0 + x ) - 1.0e0; return( dgam1_v ); } 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; if (d <= 0.e0) { dgam1_v = x*((z + 0.5e0) + 0.5e0); return( dgam1_v ); } dgam1_v = t*z/x; return( dgam1_v ); } else { /* CASE WHEN PRECISION IS LESS THAN 14 DIGITS * */ switch (ARITHIF(t)) { case -1: goto L_300; case 0: goto L_100; case 1: goto L_200; } L_100: dgam1_v = 0.0; return( dgam1_v ); /* ----------- * * CASE WHEN 0 .LT. T .LE. 0.5 * * W IS A MINIMAX APPROXIMATION FOR * THE SERIES A(15) + A(16)*T + ... * * ----------- */ L_200: if (t > 0.5e0) { if (x > xmax) { dgam1_v = -1.0e0; return( dgam1_v ); } dgam1_v = 1.0e0/dgamma( 1.0e0 + x ) - 1.0e0; return( dgam1_v ); } 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 <= 0.0) { dgam1_v = x*z; return( dgam1_v ); } dgam1_v = (t/x)*((z - 0.5) - 0.5); return( dgam1_v ); /* ----------- * * CASE WHEN -0.5 .LE. T .LT. 0 * * W IS A MINIMAX APPROXIMATION FOR * THE SERIES A(15) + A(16)*T + ... * * ----------- */ L_300: if (t < -0.5e0) { if (x > xmax) { dgam1_v = -1.0e0; return( dgam1_v ); } dgam1_v = 1.0e0/dgamma( 1.0e0 + x ) - 1.0e0; return( dgam1_v ); } 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 <= 0.0) { dgam1_v = x*((z + 0.5) + 0.5); return( dgam1_v ); } dgam1_v = t*z/x; return( dgam1_v ); } return( dgam1_v ); } /* end of function */