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