/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:07 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dgami.h" #include #include /* PARAMETER translations */ #define CLOG10 2.30258509299404568401799145468e0 #define D00 (-1.0e0/3.0e0) #define D01 (1.0e0/12.0e0) #define D02 (-2.0e0/135.0e0) #define D03 (1.0e0/864.0e0) #define D04 (1.0e0/2835.0e0) #define D10 (-1.0e0/540.0e0) #define D11 (-1.0e0/288.0e0) #define D12 (1.0e0/378.0e0) #define RT2PIN .398942280401432677939946059935e0 #define RTPI 1.77245385090551602729816748334e0 #define RTPIN (1.0e0/RTPI) /* end of PARAMETER translations */ /* COMMON translations */ struct t_dgamic { double ptol, qtol, xerr, pqerr, round; long int msgoff; } dgamic; /* end of COMMON translations */ void /*FUNCTION*/ dgami( double a, double x, double *pans, double *qans, long *ierr) { long int i, iop, m, n, _i, _r; double a2n, a2nm1, acc, amn, apn, b2n, b2nm1, c, c0, c1, c2, c3, c4, c5, c6, c7, c8, g, h, j, l, r, rta, rtx, s, sum, t, tol, twoa, u, w, wk[20], y, z; static double a0[4], a1[4], a2[2], a3[2], a4[2], a5[2], a6[2], a7[2], a8[2], b0[6], b1[4], b2[5], b3[5], b4[4], b5[3], b6[2], b7[2], d0[6], d1[4], d2[2], d3[2], d4[1], d5[1], d6[1], e8, r30; static LOGICAL32 first = TRUE; static double acc0[5]={5.0e-30,5.0e-17,5.0e-15,5.0e-7,5.0e-4}; static double big[5]={50.0e0,30.0e0,25.0e0,14.0e0,10.0e0}; static double e0[5]={0.0e0,0.0e0,0.25e-3,0.25e-1,0.14e0}; static double x0[5]={68.0e0,45.0e0,31.0e0,17.0e0,9.7e0}; static double x1[5]={2.0e0,2.0e0,1.1e0,1.1e0,1.1e0}; static double d20 = .413359788359788e-02; static double d30 = .649434156378601e-03; static double d40 = -.861888290916712e-03; static double d50 = -.336798553366358e-03; static double d60 = .531307936463992e-03; static double d70 = .344367606892378e-03; static double d80 = -.652623918595309e-03; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const A0 = &a0[0] - 1; double *const A1 = &a1[0] - 1; double *const A2 = &a2[0] - 1; double *const A3 = &a3[0] - 1; double *const A4 = &a4[0] - 1; double *const A5 = &a5[0] - 1; double *const A6 = &a6[0] - 1; double *const A7 = &a7[0] - 1; double *const A8 = &a8[0] - 1; double *const Acc0 = &acc0[0] - 1; double *const B0 = &b0[0] - 1; double *const B1 = &b1[0] - 1; double *const B2 = &b2[0] - 1; double *const B3 = &b3[0] - 1; double *const B4 = &b4[0] - 1; double *const B5 = &b5[0] - 1; double *const B6 = &b6[0] - 1; double *const B7 = &b7[0] - 1; double *const Big = &big[0] - 1; double *const D0 = &d0[0] - 1; double *const D1 = &d1[0] - 1; double *const D2 = &d2[0] - 1; double *const D3 = &d3[0] - 1; double *const D4 = &d4[0] - 1; double *const D5 = &d5[0] - 1; double *const D6 = &d6[0] - 1; double *const E0 = &e0[0] - 1; double *const Wk = &wk[0] - 1; double *const X0 = &x0[0] - 1; double *const X1 = &x1[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ A0[1] = -.231272501940775e-02; A0[2] = -.335378520024220e-01; A0[3] = -.159840143443990e00; A0[4] = -.333333333333333e00; B0[1] = .633763414209504e-06; B0[2] = -.939001940478355e-05; B0[3] = .239521354917408e-02; B0[4] = .376245718289389e-01; B0[5] = .238549219145773e00; B0[6] = .729520430331981e00; A1[1] = -.398783924370770e-05; A1[2] = -.587926036018402e-03; A1[3] = -.491687131726920e-02; A1[4] = -.185185185184291e-02; B1[1] = .386325038602125e-02; B1[2] = .506042559238939e-01; B1[3] = .283344278023803e00; B1[4] = .780110511677243e00; A2[1] = .669564126155663e-03; A2[2] = .413359788442192e-02; B2[1] = -.421924263980656e-03; B2[2] = .650837693041777e-02; B2[3] = .682034997401259e-01; B2[4] = .339173452092224e00; B2[5] = .810647620703045e00; A3[1] = .810586158563431e-03; A3[2] = .649434157619770e-03; B3[1] = -.632276587352120e-03; B3[2] = .905375887385478e-02; B3[3] = .906610359762969e-01; B3[4] = .406288930253881e00; B3[5] = .894800593794972e00; A4[1] = -.105014537920131e-03; A4[2] = -.861888301199388e-03; B4[1] = .322609381345173e-01; B4[2] = .178295773562970e00; B4[3] = .591353097931237e00; B4[4] = .103151890792185e01; A5[1] = -.435211415445014e-03; A5[2] = -.336806989710598e-03; B5[1] = .178716720452422e00; B5[2] = .600380376956324e00; B5[3] = .108515217314415e01; A6[1] = -.182503596367782e-03; A6[2] = .531279816209452e-03; B6[1] = .345608222411837e00; B6[2] = .770341682526774e00; A7[1] = .443219646726422e-03; A7[2] = .344430064306926e-03; B7[1] = .821824741357866e00; B7[2] = .115029088777769e01; A8[1] = .878371203603888e-03; A8[2] = -.686013280418038e-03; D0[1] = D01; D0[2] = D02; D0[3] = D03; D0[4] = D04; D0[5] = -.178755144032922e-03; D0[6] = .391926317852244e-04; D1[1] = D11; D1[2] = D12; D1[3] = -.990226337448560e-03; D1[4] = .205761316872428e-03; D2[1] = -.268132716049383e-02; D2[2] = .771604938271605e-03; D3[1] = .229472093621399e-03; D3[2] = -.469189494395256e-03; D4[1] = .784039221720067e-03; D5[1] = -.697281375836586e-04; D6[1] = -.592166437353694e-03; _aini = 0; } /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. * *>> 1998-10-29 DGAMI Krogh Moved external statement up for mangle. *>> 1996-03-30 DGAMI Krogh ALOG10 => CLOG10 *>> 1994-11-23 DGAMI CLL Added ?GR17 & ?GR29 to list of "c--" names. *>> 1994-11-18 DGAMI WV Snyder del block data. Change error reporting. *>> 1994-11-03 DGAMI Krogh Changed F90 do's to F77 do's. *>> 1994-10-20 DGAMI Krogh Changes to use M77CON *>> 1994-05-09 DGAMI WVS JPL Make SP and DP from same source *>> 1994-03-14 DGAMI WVS JPL Repair incomplete error message *>> 1993-07-21 DGAMI WVS JPL Conversion from NSWC to Math 77 *--D replaces "?": ?GAMI, ?GAMIC, ?GAMIE, ?GAMIK, ?ERF, ?ERFC, ?ERFCE, *--& ?GAM1, ?RCOMP, ?RLOG, ?RLOG1, ?REXP, ?GAMIB, *--& ?GR17, ?GR29, ?ERM1, ?ERV1 * ---------------------------------------------------------------------- * Evaluation of the incomplete gamma function ratios * P(A,X) and Q(A,X), where P(A,X) = gamma(A,X)/GAMMA(A), * and P(A,X) + Q(A,X) = 1. * * ---------- * * IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X * ARE NOT BOTH 0. * * PANS AND QANS ARE VARIABLES. DGAMI ASSIGNS PANS THE VALUE * P(A,X) AND QANS THE VALUE Q(A,X). * * NORMAL RETURN ... * * UPON NORMAL RETURN, BOTH OF PANS AND QANS WILL BE IN THE RANGE * [0.0 .. 1.0], PANS + QANS WILL BE WITHIN ROUND-OFF OF 1.0, AND * IERR = 0. * * ERROR RETURN ... * * PANS IS ASSIGNED THE VALUE 3 WHEN A AND X ARE BOTH ZERO AND THE * VALUE 4 WHEN ONE OR BOTH OF A OR X IS NEGATIVE. IERR = INT(PANS) * WHEN PANS > 1. IERR is assigned the value 2 if the desired error * was not achieved. * * To set tolerances, CALL DGAMIK, q.v. To discover the error * committed on the last call to DGAMI, CALL DGAMIE, q.v. * * ---------------------------------------------------------------------- * WRITTEN BY ALFRED H. MORRIS, JR. * NAVAL SURFACE WARFARE CENTER, DAHLGREN, VIRGINIA * REVISED ... DEC 1991 * See "Computation of the Incomplete Gamma Function", ACM TOMS 12, 4 * (Dec 1986) 377-393. * Adapted to Math 77 by W. Van Snyder, JPL, 1993 May 5. * ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* CLOG10 = LN(10) * RT2PIN = 1/SQRT(2*PI) * RTPI = SQRT(PI) * RTPIN = 1/SQRT(PI) */ /* ------------------------ * SOME OF THE TEMME COEFFICIENTS: */ /* --- COMMON /DGAMIC/ ---- (see DGAMIK for discussion) * PQERR is approximately X * ABS(derivative of P(a,x) with respect to * X). PQERR is used to estimate the error in P and Q due to * uncertainty in X: Error in P or Q = PQERR * (relative error in * X). When A and X are equal, PQERR = sqrt(A / (2*Pi) ). Since P * and Q are bounded by 1, there can be no precision in the result if * A or X has relative error eps, and A > 2*Pi/eps**2. If A and X * are not nearly equal, PQERR is small. */ /* ------------------------ */ /* ------------------------ * * COEFFICIENTS FOR MINIMAX APPROXIMATIONS * FOR C0,...,C8 * * ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ * * COEFFICIENTS FOR THE TEMME EXPANSION * * ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ */ /* ------------------------ * */ if (first) { first = FALSE; dgamib(); dgamik( 0.0e0, 0.0e0, -1.0e0, 0 ); e8 = 8.0e0*dgamic.round; r30 = fmax( dgamic.round, 1.0e-30 ); } if (a == 0.5e0) { if (x < 0.0e0) goto L_510; /* Equation 8 in ref. */ rtx = sqrt( x ); dgamic.pqerr = rtx*RTPIN*exp( -x ); if (x < 0.25e0) { *pans = derf( rtx ); goto L_395; } *qans = derfc( rtx ); goto L_390; } acc = fmax( r30, fmin( fabs( dgamic.ptol ), fabs( dgamic.qtol ) ) ); iop = 1; for (n = 2; n <= 5; n++) { if (acc >= Acc0[n]) iop = n; } acc = fmin( Acc0[5], acc ); /* SELECT THE APPROPRIATE ALGORITHM * */ if (a >= 1.0e0) { if (x < 0.0e0) goto L_510; if (x == 0.0e0) goto L_300; if (a < Big[iop]) { if (a <= x && x < X0[iop]) { twoa = a + a; m = twoa; if (twoa == m) { /* FINITE SUMS FOR Q WHEN A .GE. 1 AND 2*A IS AN INTEGER * Equation 14 in ref. * */ i = m/2; t = exp( -x ); c = a - i; if (c == 0.0e0) { sum = t; } else { rtx = sqrt( x ); t = 2.0e0*t*rtx/RTPI; sum = derfc( rtx ) + t; } for (n = 2; n <= i; n++) { c += 1.0e0; t = (x*t)/c; sum += t; } dgamic.pqerr = x*t; *qans = sum; goto L_390; } } } else { l = x/a; if (l == 0.0e0) goto L_300; if (l <= 0.5e0) { s = (l - 0.5e0) - 0.5e0; z = drlog( l ); } else { s = x - a; if (s + a != x) { /* Exponents of S and A are different */ if (x > a) { s = (0.5e0*x - a) + 0.5e0*x; } else { s = (x - 0.5e0*a) - 0.5e0*a; } } s /= a; z = drlog1( s ); } rta = sqrt( a ); if (z >= 700.0e0/a) goto L_330; y = a*z; if (fabs( s ) <= 0.4e0) { /* MINIMAX APPROXIMATIONS * */ u = 1.0e0/a; z = sqrt( z + z ); if (l < 1.0e0) z = -z; c = exp( -y ); w = 0.5e0*derfce( sqrt( y ) ); if (rta*fabs( s ) > E0[iop]) { /* Equations 17 in ref. */ switch (iop) { case 1: goto L_210; case 2: goto L_220; case 3: goto L_230; case 4: goto L_240; case 5: goto L_250; } L_210: ; dgr29( z, u, &t ); goto L_290; L_220: ; dgr17( z, u, &t ); goto L_290; L_230: if (fabs( s ) <= 1.0e-3) goto L_260; /* USING THE MINIMAX APPROXIMATIONS * */ c0 = (((A0[1]*z + A0[2])*z + A0[3])*z + A0[4])/ ((((((B0[1]*z + B0[2])*z + B0[3])*z + B0[4])* z + B0[5])*z + B0[6])*z + 1.0); c1 = (((A1[1]*z + A1[2])*z + A1[3])*z + A1[4])/ ((((B1[1]*z + B1[2])*z + B1[3])*z + B1[4])*z + 1.0); c2 = (A2[1]*z + A2[2])/(((((B2[1]*z + B2[2])*z + B2[3])*z + B2[4])*z + B2[5])*z + 1.0); c3 = (A3[1]*z + A3[2])/(((((B3[1]*z + B3[2])*z + B3[3])*z + B3[4])*z + B3[5])*z + 1.0); c4 = (A4[1]*z + A4[2])/((((B4[1]*z + B4[2])*z + B4[3])*z + B4[4])*z + 1.0); c5 = (A5[1]*z + A5[2])/(((B5[1]*z + B5[2])*z + B5[3])*z + 1.0); c6 = (A6[1]*z + A6[2])/((B6[1]*z + B6[2])*z + 1.0); c7 = (A7[1]*z + A7[2])/((B7[1]*z + B7[2])*z + 1.0); c8 = A8[1]*z + A8[2]; t = (((((((c8*u + c7)*u + c6)*u + c5)*u + c4)* u + c3)*u + c2)*u + c1)*u + c0; goto L_290; /* TEMME EXPANSION * */ L_240: ; c0 = (((((D0[6]*z + D0[5])*z + D0[4])*z + D0[3])* z + D0[2])*z + D0[1])*z + D00; c1 = (((D1[4]*z + D1[3])*z + D1[2])*z + D1[1])* z + D10; c2 = D2[1]*z + d20; t = (c2*u + c1)*u + c0; goto L_290; L_250: ; t = ((D0[3]*z + D0[2])*z + D0[1])*z + D00; goto L_290; } /* TEMME EXPANSION FOR L = 1 (Equations 17-18 in ref). * (Can't get here if IOP = 1 or 2 because E0(1..2) = 0.0). * */ switch (IARITHIF(iop - 4)) { case -1: goto L_260; case 0: goto L_270; case 1: goto L_280; } L_260: c0 = ((D0[3]*z + D0[2])*z + D0[1])*z + D00; c1 = ((D1[3]*z + D1[2])*z + D1[1])*z + D10; c2 = (D2[2]*z + D2[1])*z + d20; c3 = (D3[2]*z + D3[1])*z + d30; c4 = D4[1]*z + d40; c5 = D5[1]*z + d50; c6 = D6[1]*z + d60; t = (((((((d80*u + d70)*u + c6)*u + c5)*u + c4)*u + c3)*u + c2)*u + c1)*u + c0; goto L_290; L_270: c0 = (D0[2]*z + D0[1])*z + D00; c1 = D1[1]*z + D10; t = (d20*u + c1)*u + c0; goto L_290; L_280: t = D0[1]*z + D00; L_290: dgamic.pqerr = c*rta*RT2PIN; if (l >= 1.0e0) { *qans = c*(w + RT2PIN*t/rta); goto L_390; } *pans = c*(w - RT2PIN*t/rta); goto L_395; } } r = drcomp( a, x ); if (r == 0.0e0) goto L_331; dgamic.pqerr = r*sqrt( a )*RT2PIN; if (x <= fmax( a, CLOG10 )) { /* TAYLOR SERIES FOR P/R (Equation 15 in ref.) * */ apn = a + 1.0e0; t = x/apn; for (n = 1; n <= 20; n++) { Wk[n] = t; apn += 1.0e0; t *= x/apn; if (t <= 1.0e-3) goto L_60; } n = 20; L_60: sum = t; tol = 0.5e0*acc; L_61: ; apn += 1.0e0; t *= x/apn; sum += t; if (t > tol) goto L_61; for (m = n; m >= 1; m--) { sum += Wk[m]; } *pans = (r/a)*(1.0e0 + sum); goto L_395; } if (x < X0[iop]) goto L_170; /* ASYMPTOTIC EXPANSION (Equation 15 in ref). * */ amn = a - 1.0e0; t = amn/x; for (n = 1; n <= 20; n++) { Wk[n] = t; amn -= 1.0e0; t *= amn/x; if (fabs( t ) <= 1.0e-3) goto L_90; } n = 20; L_90: sum = t; L_91: if (fabs( t ) >= acc) { amn -= 1.0e0; t *= amn/x; sum += t; goto L_91; } for (m = n; m >= 1; m--) { sum += Wk[m]; } *qans = (r/x)*(1.0e0 + sum); goto L_390; } /* A .LT. 1.0 HERE * */ if (x < X1[iop]) { if (a < 0.0e0 || x < 0.0e0) goto L_510; if (a*x == 0.0e0) goto L_320; /* TAYLOR SERIES FOR P(A,X)/X**A (Equations 9-10 in ref.) * */ l = 3.0e0; c = x; sum = x/(a + 3.0e0); tol = 3.0e0*acc/(a + 1.0e0); L_120: ; l += 1.0e0; c *= -(x/l); t = c/(a + l); sum += t; if (fabs( t ) > tol) goto L_120; j = a*x*((sum/6.0e0 - 0.5e0/(a + 2.0e0))*x + 1.0e0/(a + 1.0e0)); z = a*log( x ); h = dgam1( a ); g = 1.0e0 + h; dgamic.pqerr = a*sqrt( a )*exp( z - x )*RT2PIN; u = exp( z ); /* Equation 9 in ref. */ *pans = u*g*(0.5e0 + (0.5e0 - j)); if (*pans <= 0.9e0) goto L_395; l = drexp( z ); /* Equation 10 in ref. */ *qans = fmax( (u*j - l)*g - h, 0.0e0 ); goto L_390; } /* A .LT. 1.0 AND X .GE. X1(IOP) here (BUT A .NE. 0.5) * */ if (a < 0.0e0) goto L_510; if (a == 0.0e0) goto L_310; r = drcomp( a, x ); dgamic.pqerr = r*sqrt( a )*RT2PIN; if (r == 0.0e0) goto L_310; /* CONTINUED FRACTION EXPANSION (Equation 11 in ref.) * a .le. 1 and x .ge. 1.1, or 1 .le. a .le. big and * a .le. x .lt. x0, or a .gt. big and 7*a/5 .le. x .lt. x0. * */ L_170: tol = fmax( e8, 4.0e0*acc ); b2n = x + (1.0e0 - a); b2nm1 = x/b2n; a2nm1 = 1.0e0/b2n; a2n = a2nm1; c = 1.0e0; L_180: ; a2nm1 = x*a2n + c*a2nm1; b2nm1 = x + c*b2nm1; c += 1.0e0; t = c - a; a2n = a2nm1 + t*a2n; b2n = b2nm1 + t; a2nm1 /= b2n; b2nm1 /= b2n; a2n /= b2n; if (fabs( a2n - a2nm1/b2nm1 ) >= tol*a2n) goto L_180; *qans = r*a2n; goto L_390; /* SPECIAL CASES * */ L_330: if (a == x) { *qans = 0.5e0; dgamic.pqerr = rta*RT2PIN; goto L_390; } L_331: if (x > a) goto L_310; /* The error magnification factor is identically zero along both axes. */ L_300: *pans = 0.0e0; *qans = 1.0e0; dgamic.pqerr = 0.0; goto L_400; L_320: if (a != 0.0e0) goto L_300; if (x == 0.0e0) goto L_500; L_310: *qans = 0.0e0; dgamic.pqerr = 0.0; /* Check the error committed * */ L_390: *pans = 0.5e0 + (0.5e0 - *qans); goto L_400; L_395: *qans = 0.5e0 + (0.5e0 - *pans); L_400: dgamic.pqerr *= dgamic.xerr; if (dgamic.pqerr > 4.0e0*dgamic.round) { if (dgamic.ptol > 0.0) { if (dgamic.pqerr > dgamic.ptol**pans) goto L_520; } else { if (dgamic.pqerr > -dgamic.ptol) goto L_520; } if (dgamic.qtol > 0.0) { if (dgamic.pqerr > dgamic.qtol**qans) goto L_520; } else { if (dgamic.pqerr > -dgamic.qtol) goto L_520; } } *ierr = 0; return; /* ERROR RETURNS * */ L_500: ermsg( "DGAMI", 3, 2 + dgamic.msgoff, "A and X are both zero", '.' ); *pans = 3.0; *ierr = 3; return; L_510: derm1( "DGAMI", 4, 2 + dgamic.msgoff, "One of A or X is negative" , "A", a, ',' ); derv1( "X", x, '.' ); *pans = 4.0; *ierr = 4; return; L_520: derm1( "DGAMI", 2, 2 + dgamic.msgoff, "Error tolerances not achieved" , "Absolute error", dgamic.pqerr, ',' ); derv1( "PANS", *pans, ',' ); derv1( "PTOL", dgamic.ptol, ',' ); derv1( "QANS", *qans, ',' ); derv1( "QTOL", dgamic.qtol, ',' ); derv1( "A", a, ',' ); derv1( "X", x, '.' ); *ierr = 2; return; } /* end of function */ void /*FUNCTION*/ dgamik( double ptola, double qtola, double xerra, long msga) { long int _l0; static LOGICAL32 first = TRUE; /* Control DGAMI tolerances and error action. * * PTOLA and QTOLA are tolerances for PANS and QANS. When negative they * are absolute tolerances; when positive they are relative * tolerances; when zero they mean use the default tolerance. The * default is 4 times the round-off level. * XERRA is the relative error in X. Negative means use the default, * which is one round-off. * MSGA is the message action offset, added onto the message level for * each message produced by DGAMI. * */ /* ROUND is the round-off level. * */ /* /DGAMIC/ is used to remember and communicate the tolerance, error and * message action settings. * PTOL, QTOL are tolerances for PANS and QANS, set from PTOLA and QTOLA. * XERR is the relative error in X, set from XERRA. * PQERR is the absolute error in the last computation of P or Q, if any, * else it is -1.0. * ROUND is the round-off level. * MSGOFF is the message action offset. * FIRST is .TRUE. iff DGAMIK has never been called. * */ if (first) { dgamib(); dgamic.round = DBL_EPSILON; first = FALSE; } if (ptola == 0.0e0) { dgamic.ptol = 4.0e0*dgamic.round; } else { dgamic.ptol = ptola; } if (qtola == 0.0e0) { dgamic.qtol = 4.0e0*dgamic.round; } else { dgamic.qtol = qtola; } if (xerra < 0.0e0) { dgamic.xerr = dgamic.round; } else { dgamic.xerr = xerra; } dgamic.msgoff = msga; return; } /* end of function */ void /*FUNCTION*/ dgamie( double *pqerrv) { static LOGICAL32 first = TRUE; /* Return the error committed on the last call to DGAMI, if any, else * return -1.0 * */ if (first) { dgamib(); first = FALSE; } *pqerrv = dgamic.pqerr; return; } /* end of function */ void /*FUNCTION*/ dgamib() { static LOGICAL32 first = TRUE; if (first) { dgamic.pqerr = -1.0e0; first = FALSE; } return; } /* end of function */