#include "myf2c.h"


int ribesldll( double x[], double *rx, double y[], double *nb,
                double *alpha, double *ize, double *err );

long ribesl_(double *x, double *alpha, long *nb, long *ize, double *b, long *ncalc);



/*
**  ribesldll
**
**  This is a "wrapper" function that is called by GAUSS' dllcall function
**  and which itself calls the C source code for the Bessel function.
**  All of the input arguments are pointers to doubles.  GAUSS can also
**  pass pointers to strings, but nothing else.  Thus all arguments to
**  wrapper functions to be called using dllcall will have to be either
**  pointers to doubles or pointers to strings.
**
**  Generally, it will be necessary to pass in both column and row
**  information about the input arrays being passed from GAUSS
**  or at least their lengths because the wrapper function won't 
**  otherwise know.  In this case, the type of calculation going on
**  required that x be a vector and y a matrix with the same number
**  of rows as x but with nb columns.  This is an unusual requirement,
**  though, which arises out of the "sequential" nature of the
**  calculation.
**
**  x, y -
**
**    x is a column vector with one column; y is the output vector
**    with the same number of rows as x and nb columns.
**    These dimensions are enforced by the calling routine
**
**    For each row in x, nb modified Bessels are computed and stored in y
**    in the following order:
**
**    x[0] -> y[0], y[1],...,y[nb-1]
**    x[1] -> y[nb], y[nb+1],...,y[2*nb-1]
**    x[2] -> y[2*nb], y[2*nb+1],...,y[3*nb-1]
**     .
**     .
**     .
**
**  nb - scalar integer,   nb > 0
**
**     number of elements in sequence, also columns of y
**
**  alpha - scalar,    0 <= alpha < 1
**
**     A sequence of modified Bessels are computed:
**
**     I_(alpha), I(1+alpha),...,I(nb-1+alpha)
**
**
**  ize - scalar, if 1, unscaled Bessel's, if 2, exponential scaled
**
**  err - GAUSS error value (i.e., what is returned by error(0))
**
*/



int ribesldll( double x[], double *rx, double y[], double *nb,
                double *alpha, double *ize, double *err )
{

   long i,j,i2,rx0,ry0,cy0,nb0,ize0,ncalc;

   rx0 = (long)*rx;   /* rows of input matrix    */
   nb0 = (long)*nb;   /* columns of output matrix */
   ize0 = (long)*ize;

   for (i = 0; i < rx0; i++) {
       i2 = nb0 * i;
       ribesl_(&x[i],alpha,&nb0,&ize0,&y[i2],&ncalc);
       if (ncalc != nb0) {
           for (j = 0; j < nb0; j++) y[i2+j] = *err;
       }
   }
   return 0;
}




/*  -- translated by f2c (version 19960513).
   You must link the resulting object file with the libraries:
        -lf2c -lm   (in that order)
*/

/* Table of constant values */

static doublereal c_b4 = 1.585;

/*     ALGORITHM 597, COLLECTED ALGORITHMS FROM ACM. */
/*     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2, */
/*     JUN., 1983, P. 242-245. */
/* Subroutine */ integer ribesl_(x, alpha, nb, ize, b, ncalc)
doublereal *x, *alpha;
integer *nb, *ize;
doublereal *b;
integer *ncalc;
{
    /* Initialized data */

    static doublereal one = 1.;
    static doublereal rtnsig = 1e-4;
    static doublereal enmten = 1.2e-37;
    static doublereal two = 2.;
    static doublereal zero = 0.;
    static doublereal half = .5;
    static integer nsig = 17;
    static doublereal xlarge = 1e4;
    static doublereal exparg = 88.;
    static doublereal enten = 1e38;
    static doublereal ensig = 1e17;

    /* System generated locals */
    integer i__1;
    doublereal d__1, d__2, d__3;

     /* Local variables */
    static integer nend, magx;
    static doublereal pold;
    static integer nbmx;
    static doublereal test;
    static integer l, n;
    static doublereal p, empal, halfx, tempa, tempb, tempc, psave, plast,
            tover, emp2al, em, en;
    extern doublereal dgamma_();
    static doublereal psavel;
    static integer nstart;
    static doublereal sum;

/* ------------------------------------------------------------------- */

/*  THIS ROUTINE CALCULATES BESSEL FUNCTIONS I SUB(N+ALPHA) (X) */
/*  FOR NON-NEGATIVE ARGUMENT X, AND NON-NEGATIVE ORDER N+ALPHA, */
/*  WITH OR WITHOUT EXPONENTIAL SCALING. */


/* EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE */

/* X     - WORKING PRECISION NON-NEGATIVE REAL ARGUMENT FOR WHICH */
/*         I'S OR EXPONENTIALLY SCALED I'S (I*EXP(-X)) */
/*         ARE TO BE CALCULATED.  IF I'S ARE TO BE CALCULATED, */
/*         X MUST BE LESS THAN EXPARG (SEE BELOW). */
/* ALPHA - WORKING PRECISION FRACTIONAL PART OF ORDER FOR WHICH */
/*         I'S OR EXPONENTIALLY SCALED I'S (I*EXP(-X)) ARE */
/*         TO BE CALCULATED.  0 .LE. ALPHA .LT. 1.0. */
/* NB    - INTEGER NUMBER OF FUNCTIONS TO BE CALCULATED, NB .GT. 0. */
/*         THE FIRST FUNCTION CALCULATED IS OF ORDER ALPHA, AND THE */
/*         LAST IS OF ORDER (NB - 1 + ALPHA). */
/* IZE   - INTEGER TYPE.  IZE = 1 IF UNSCALED I'S ARE TO CALCULATED, */
/*         AND 2 IF EXPONENTIALLY SCALED I'S ARE TO BE CALCULATED. */
/* B     - WORKING PRECISION OUTPUT VECTOR OF LENGTH NB.  IF THE ROUTINE
*/
/*         TERMINATES NORMALLY (NCALC=NB), THE VECTOR B CONTAINS THE */
/*         FUNCTIONS I/ALPHA/(X) THROUGH I/NB-1+ALPHA/(X), OR THE */
/*         CORRESPONDING EXPONENTIALLY SCALED FUNCTIONS. */
/* NCALC - INTEGER OUTPUT VARIABLE INDICATING POSSIBLE ERRORS. */
/*         BEFORE USING THE VECTOR B, THE USER SHOULD CHECK THAT */
/*         NCALC=NB, I.E., ALL ORDERS HAVE BEEN CALCULATED TO */
/*         THE DESIRED ACCURACY.  SEE ERROR RETURNS BELOW. */


/* ******************************************************************* */
/* ******************************************************************* */

/* EXPLANATION OF MACHINE-DEPENDENT CONSTANTS */

/* NSIG   - DECIMAL SIGNIFICANCE DESIRED.  SHOULD BE SET TO */
/*          IFIX(ALOG10(2)*NBIT+1), WHERE NBIT IS THE NUMBER OF */
/*          BITS IN THE MANTISSA OF A WORKING PRECISION VARIABLE. */
/*          SETTING NSIG LOWER WILL RESULT IN DECREASED ACCURACY */
/*          WHILE SETTING NSIG HIGHER WILL INCREASE CPU TIME */
/*          WITHOUT INCREASING ACCURACY.  THE TRUNCATION ERROR */
/*          IS LIMITED TO A RELATIVE ERROR OF T=.5*10**(-NSIG). */
/* ENTEN  - 10.0 ** K, WHERE K IS THE LARGEST INTEGER SUCH THAT */
/*          ENTEN IS MACHINE-REPRESENTABLE IN WORKING PRECISION. */
/* ENSIG  - 10.0 ** NSIG. */
/* RTNSIG - 10.0 ** (-K) FOR THE SMALLEST INTEGER K SUCH THAT */
/*          K .GE. NSIG/4. */
/* ENMTEN - THE SMALLEST ABS(X) SUCH THAT X/4 DOES NOT UNDERFLOW. */
/* XLARGE - UPPER LIMIT ON THE MAGNITUDE OF X WHEN IZE=2.  BEAR */
/*          IN MIND THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS */
/*          OF THE BACKWARD RECURSION WILL BE EXECUTED. */
/* EXPARG - LARGEST WORKING PRECISION ARGUMENT THAT THE LIBRARY */
/*          EXP ROUTINE CAN HANDLE AND UPPER LIMIT ON THE */
/*          MAGNITUDE OF X WHEN IZE=1. */

/*     APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: */

/*         IBM/195    CDC/7600  UNIVAC/1108   VAX 11/780 (UNIX) */
/*          (D.P.)  (S.P.,RNDG)    (D.P.)     (S.P.)     (D.P.) */

/* NSIG      16         14          18           8        17 */
/* ENTEN   1.0D75     1.0E322     1.0D307     1.0E38    1.0D38 */
/* ENSIG   1.0D16     1.0E14      1.0D18      1.0E8     1.0D17 */
/* RTNSIG  1.0D-4     1.0E-4      1.0D-5      1.0E-2    1.0D-4 */
/* ENMTEN  2.2D-78    1.0E-290    1.2D-308    1.2E-37   1.2D-37 */
/* XLARGE  1.0D4      1.0E4       1.0D4       1.0E4     1.0D4 */
/* EXPARG  174.0D0    740.0E0     709.0D0     88.0E0    88.0D0 */

/* ******************************************************************* */
/* ******************************************************************* */


/* ERROR RETURNS */

/*  IN CASE OF AN ERROR,  NCALC .NE. NB,  AND NOT ALL I'S ARE */
/*  CALCULATED TO THE DESIRED ACCURACY. */

/*  NCALC .LT. 0:  AN ARGUMENT IS OUT OF RANGE. FOR EXAMPLE, */
/*     NB .LE. 0, IZE IS NOT 1 OR 2, OR IZE=1 AND ABS(X) .GE. EXPARG. */
/*     IN THIS CASE, THE B-VECTOR IS NOT CALCULATED, AND NCALC IS */
/*     SET TO  MIN0(NB,0)-1  SO THAT NCALC .NE. NB. */

/*  NB .GT. NCALC .GT. 0: NOT ALL REQUESTED FUNCTION VALUES COULD */
/*     BE CALCULATED ACCURATELY.  THIS USUALLY OCCURS BECAUSE NB IS */
/*     MUCH LARGER THAN ABS(X).  IN THIS CASE, B(N) IS CALCULATED */
/*     TO THE DESIRED ACCURACY FOR  N .LE. NCALC,  BUT PRECISION */
/*     IS LOST FOR NCALC .LT. N .LE. NB.  IF B(N) DOES NOT VANISH */
/*     FOR  N .GT. NCALC  (BECAUSE IT IS TOO SMALL TO BE REPRESENTED), */
/*     AND  B(N)/B(NCALC) = 10**(-K), THEN ONLY THE FIRST NSIG-K */
/*     SIGNIFICANT FIGURES OF B(N) CAN BE TRUSTED. */


/* OTHER SUBPROGRAMS REQUIRED (SINGLE PRECISION VERSION) */

/*     EXP,GAMMA,AMAX1,SQRT,FLOAT,IFIX,MIN0 */

/* OTHER SUBPROGRAMS REQUIRED (DOUBLE PRECISION VERSION) */

/*     DBLE,DEXP,DGAMMA,DMAX1,DSQRT,FLOAT,IFIX,MIN0,SNGL */


/* ACKNOWLEDGEMENT */

/*  THIS PROGRAM IS BASED ON A PROGRAM WRITTEN BY DAVID J. */
/*  SOOKNE (2) THAT COMPUTES VALUES OF THE BESSEL FUNCTIONS J OR */
/*  I OF REAL ARGUMENT AND INTEGER ORDER.  MODIFICATIONS INCLUDE */
/*  THE RESTRICTION OF THE COMPUTATION TO THE I BESSEL FUNCTION */
/*  OF NON-NEGATIVE REAL ARGUMENT, THE EXTENSION OF THE COMPUTATION */
/*  TO ARBITRARY POSITIVE ORDER, THE INCLUSION OF OPTIONAL */
/*  EXPONENTIAL SCALING, AND THE ELIMINATION OF MOST UNDERFLOW. */

/* REFERENCES */

/*  (1) OLVER, F. W. J., AND SOOKNE, D. J., "A NOTE ON BACKWARD */
/*        RECURRENCE ALGORITHMS," MATH. COMP. 26, 1972, PP 941-947. */

/*  (2) SOOKNE, D. J., "BESSEL FUNCTIONS OF REAL ARGUMENT AND */
/*        INTEGER ORDER," NBS JOUR. OF RES. B. 77B, 1973, PP 125-132. */


/*      MODIFIED BY: W. J. CODY */
/*                   APPLIED MATHEMATICS DIVISION */
/*                   ARGONNE NATIONAL LABORATORY */
/*                   ARGONNE, IL  60439 */

/*      LATEST MODIFICATION: MAY 18, 1982 */

/* ------------------------------------------------------------------- */
/* S    REAL              ALPHA,B,EM,EMPAL,EMP2AL,EN,ENMTEN,ENSIG, */
/* S   2 ENTEN,EXPARG,GAMMA,HALF,HALFX,ONE,P,PLAST,POLD,PSAVE,PSAVEL, */
/* S   3 RTNSIG,SUM,TEMPA,TEMPB,TEMPC,TEST,TOVER,TWO,X,XLARGE,ZERO */
/* ------------------------------------------------------------------- */
/*  MATHEMATICAL CONSTANTS */
/* ------------------------------------------------------------------- */
/* S    DATA ONE,TWO,ZERO,HALF/1.0E0,2.0E0,0.0E0,0.5E0/ */
    /* Parameter adjustments */
    --b;

    /* Function Body */
/* ------------------------------------------------------------------- */
/*  MACHINE DEPENDENT PARAMETERS */
/* ------------------------------------------------------------------- */
/* S    DATA NSIG,XLARGE,EXPARG / 7,1.0E4,88.0E0/ */
/* S    DATA ENTEN,ENSIG,RTNSIG/1.0E38,1.0E7,1.0E-2/ */
/* S    DATA ENMTEN/1.2E-37/ */
/* ------------------------------------------------------------------- */
/* S    MAGX = IFIX(X) */
    magx = (integer) ((real) (*x));
    if (*nb > 0 && *x >= zero && *alpha >= zero && *alpha < one && (*ize == 1
            && *x <= exparg || *ize == 2 && *x <= xlarge)) {
        goto L10;
    }
/* ------------------------------------------------------------------- */
/* ERROR RETURN -- X,NB,OR IZE IS OUT OF RANGE */
/* ------------------------------------------------------------------- */
    *ncalc = min(*nb,0) - 1;
    return 0;
/* ------------------------------------------------------------------- */
/* USE 2-TERM ASCENDING SERIES FOR SMALL X */
/* ------------------------------------------------------------------- */
L10:
    *ncalc = *nb;
    if (*x < rtnsig) {
        goto L210;
    }
/* ------------------------------------------------------------------- */
/* INITIALIZE THE FORWARD SWEEP, THE P-SEQUENCE OF OLVER */
/* ------------------------------------------------------------------- */
    nbmx = *nb - magx;
    n = magx + 1;
/* S    EN = FLOAT(N+N) + (ALPHA+ALPHA) */
    en = (doublereal) ((real) (n + n)) + (*alpha + *alpha);
    plast = one;
    p = en / *x;
/* ------------------------------------------------------------------- */
/* CALCULATE GENERAL SIGNIFICANCE TEST */
/* ------------------------------------------------------------------- */
    test = ensig + ensig;
/* S    IF (2*MAGX .GT. 5*NSIG) TEST = SQRT(TEST*P) */
    if (magx << 1 > nsig * 5) {
        test = sqrt(test * p);
    }
/* S    IF (2*MAGX .LE. 5*NSIG) TEST = TEST / 1.585E0**MAGX */
    if (magx << 1 <= nsig * 5) {
        test /= pow_di(&c_b4, &magx);
    }
    if (nbmx < 3) {
        goto L30;
    }
/* ------------------------------------------------------------------- */
/* CALCULATE P-SEQUENCE UNTIL N = NB-1.  CHECK FOR POSSIBLE OVERFLOW. */
/* ------------------------------------------------------------------- */
    tover = enten / ensig;
    nstart = magx + 2;
    nend = *nb - 1;
    i__1 = nend;
    for (n = nstart; n <= i__1; ++n) {
        en += two;
        pold = plast;
        plast = p;
        p = en * plast / *x + pold;
        if (p > tover) {
            goto L40;
        }
/* L20: */
    }
    n = nend;
/* S    EN = FLOAT(N+N) + (ALPHA+ALPHA) */
    en = (doublereal) ((real) (n + n)) + (*alpha + *alpha);
/* ------------------------------------------------------------------- */
/* CALCULATE SPECIAL SIGNIFICANCE TEST FOR NBMX.GT.2. */
/* ------------------------------------------------------------------- */
/* S    TEST = AMAX1(TEST,SQRT(PLAST*ENSIG)*SQRT(P+P)) */
/* Computing MAX */
    d__1 = test, d__2 = sqrt(plast * ensig) * sqrt(p + p);
    test = max(d__1,d__2);
/* ------------------------------------------------------------------- */
/* CALCULATE P-SEQUENCE UNTIL SIGNIFICANCE TEST PASSES */
/* ------------------------------------------------------------------- */
L30:
    ++n;
    en += two;
    pold = plast;
    plast = p;
    p = en * plast / *x + pold;
    if (p < test) {
        goto L30;
    }
    goto L80;
/* ------------------------------------------------------------------- */
/* TO AVOID OVERFLOW, DIVIDE P-SEQUENCE BY TOVER.  CALCULATE */
/* P-SEQUENCE UNTIL ABS(P).GT.1. */
/* ------------------------------------------------------------------- */
L40:
    tover = enten;
    p /= tover;
    plast /= tover;
    psave = p;
    psavel = plast;
    nstart = n + 1;
L50:
    ++n;
    en += two;
    pold = plast;
    plast = p;
    p = en * plast / *x + pold;
    if (p <= one) {
        goto L50;
    }
    tempb = en / *x;
/* ------------------------------------------------------------------- */
/* CALCULATE BACKWARD TEST, AND FIND NCALC, THE HIGHEST N */
/* SUCH THAT THE TEST IS PASSED. */
/* ------------------------------------------------------------------- */
    test = pold * plast * (half - half / (tempb * tempb)) / ensig;
    p = plast * tover;
    --n;
    en -= two;
    nend = min(*nb,n);
    i__1 = nend;
    for (l = nstart; l <= i__1; ++l) {
        *ncalc = l;
        pold = psavel;
        psavel = psave;
        psave = en * psavel / *x + pold;
        if (psave * psavel > test) {
            goto L70;
        }
/* L60: */
    }
    *ncalc = nend + 1;
L70:
    --(*ncalc);
/* ------------------------------------------------------------------- */
/* INITIALIZE THE BACKWARD RECURSION AND THE NORMALIZATION SUM */
/* ------------------------------------------------------------------- */
L80:
    ++n;
    en += two;
    tempb = zero;
    tempa = one / p;
/* S    EM = FLOAT(N) - ONE */
    em = (doublereal) ((real) n) - one;
    empal = em + *alpha;
    emp2al = em - one + (*alpha + *alpha);
    sum = tempa * empal * emp2al / em;
    nend = n - *nb;
    if (nend < 0) {
        goto L130;
    } else if (nend == 0) {
        goto L110;
    } else {
        goto L90;
    }
/* ------------------------------------------------------------------- */
/* RECUR BACKWARD VIA DIFFERENCE EQUATION, CALCULATING (BUT */
/* NOT STORING) B(N), UNTIL N = NB. */
/* ------------------------------------------------------------------- */
L90:
    i__1 = nend;
    for (l = 1; l <= i__1; ++l) {
        --n;
        en -= two;
        tempc = tempb;
        tempb = tempa;
        tempa = en * tempb / *x + tempc;
        em -= one;
        emp2al -= one;
        if (n == 1) {
            goto L110;
        }
        if (n == 2) {
            emp2al = one;
        }
        empal -= one;
        sum = (sum + tempa * empal) * emp2al / em;
/* L100: */
    }
/* ------------------------------------------------------------------- */
/* STORE B(NB) */
/* ------------------------------------------------------------------- */
L110:
    b[n] = tempa;
    if (*nb > 1) {
        goto L120;
    }
    sum = sum + sum + tempa;
    goto L190;
/* ------------------------------------------------------------------- */
/* CALCULATE AND STORE B(NB-1) */
/* ------------------------------------------------------------------- */
L120:
    --n;
    en -= two;
    b[n] = en * tempa / *x + tempb;
    if (n == 1) {
        goto L180;
    }
    em -= one;
    emp2al -= one;
    if (n == 2) {
        emp2al = one;
    }
    empal -= one;
    sum = (sum + b[n] * empal) * emp2al / em;
    goto L150;
/* ------------------------------------------------------------------- */
/* N.LT.NB, SO STORE B(N) AND SET HIGHER ORDERS TO ZERO */
/* ------------------------------------------------------------------- */
L130:
    b[n] = tempa;
    nend = -nend;
    i__1 = nend;
    for (l = 1; l <= i__1; ++l) {
        b[n + l] = zero;
/* L140: */
    }
L150:
    nend = n - 2;
    if (nend == 0) {
        goto L170;
    }
/* ------------------------------------------------------------------- */
/* CALCULATE VIA DIFFERENCE EQUATION AND STORE B(N), UNTIL N = 2 */
/* ------------------------------------------------------------------- */
    i__1 = nend;
    for (l = 1; l <= i__1; ++l) {
        --n;
        en -= two;
        b[n] = en * b[n + 1] / *x + b[n + 2];
        em -= one;
        emp2al -= one;
        if (n == 2) {
            emp2al = one;
        }
        empal -= one;
        sum = (sum + b[n] * empal) * emp2al / em;
/* L160: */
    }
/* ------------------------------------------------------------------- */
/* CALCULATE B(1) */
/* ------------------------------------------------------------------- */
L170:
    b[1] = two * empal * b[2] / *x + b[3];
L180:
    sum = sum + sum + b[1];
/* ------------------------------------------------------------------- */
/* NORMALIZE.  DIVIDE ALL B(N) BY SUM. */
/* ------------------------------------------------------------------- */
L190:
/* S    IF (ALPHA.NE.ZERO)SUM=SUM*GAMMA(ONE+ALPHA)*(X*HALF)**(-ALPHA) */
    if (*alpha != zero) {
        d__1 = one + *alpha;
        d__2 = *x * half;
        d__3 = -(*alpha);
        sum = sum * dgamma_(&d__1) * pow_dd(&d__2, &d__3);
    }
/* S    IF (IZE .EQ. 1) SUM = SUM * EXP(-X) */
    if (*ize == 1) {
        sum *= exp(-(*x));
    }
    tempa = enmten;
    if (sum > one) {
        tempa *= sum;
    }
    i__1 = *nb;
    for (n = 1; n <= i__1; ++n) {
        if (b[n] < tempa) {
            b[n] = zero;
        }
        b[n] /= sum;
/* L200: */
    }
    return 0;
/* ------------------------------------------------------------------- */
/* TWO-TERM ASCENDING SERIES FOR SMALL X */
/* ------------------------------------------------------------------- */
L210:
    tempa = one;
    empal = one + *alpha;
    halfx = zero;
    if (*x > enmten) {
        halfx = half * *x;
    }
/* S    IF (ALPHA .NE. ZERO) TEMPA = HALFX ** ALPHA / GAMMA(EMPAL) */
    if (*alpha != zero) {
        tempa = pow_dd(&halfx, alpha) / dgamma_(&empal);
    }
/* S    IF (IZE .EQ. 2) TEMPA = TEMPA * EXP(-X) */
    if (*ize == 2) {
        tempa *= exp(-(*x));
    }
    tempb = zero;
    if (*x + one > one) {
        tempb = halfx * halfx;
    }
    b[1] = tempa + tempa * tempb / empal;
    if (*x != zero && b[1] == zero) {
        *ncalc = 0;
    }
    if (*nb == 1) {
        goto L250;
    }
    if (*x > zero) {
        goto L230;
    }
    i__1 = *nb;
    for (n = 2; n <= i__1; ++n) {
        b[n] = zero;
/* L220: */
    }
    goto L250;
/* ------------------------------------------------------------------- */
/* CALCULATE HIGHER ORDER FUNCTIONS */
/* ------------------------------------------------------------------- */
L230:
    tempc = halfx;
    tover = (enmten + enmten) / *x;
    if (tempb != zero) {
        tover = enmten / tempb;
    }
    i__1 = *nb;
    for (n = 2; n <= i__1; ++n) {
        tempa /= empal;
        empal += one;
        tempa *= tempc;
        if (tempa <= tover * empal) {
            tempa = zero;
        }
        b[n] = tempa + tempa * tempb / empal;
        if (b[n] == zero && *ncalc > n) {
            *ncalc = n - 1;
        }
/* L240: */
    }
L250:
    return 0;
/* ---------- LAST CARD OF RIBESL ---------- */
} /* ribesl_ */

doublereal dgamma_(x)
doublereal *x;
{
    /* Initialized data */

    static doublereal one = 1.;
    static doublereal half = .5;
    static doublereal twelve = 12.;
    static doublereal zero = 0.;
    static doublereal sqrtpi = .9189385332046727417803297;
    static doublereal pi = 3.1415926535897932384626434;
    static doublereal xbig = 34.844;
    static doublereal xminin = 5.883e-39;
    static doublereal eps = 1.388e-16;
    static doublereal xinf = 1.7014e38;
    static doublereal p[8] = { -1.71618513886549492533811,
            24.7656508055759199108314,-379.804256470945635097577,
            629.331155312818442661052,866.966202790413211295064,
            -31451.2729688483675254357,-36144.4134186911729807069,
            66456.1438202405440627855 };
    static doublereal q[8] = { -30.8402300119738975254353,
            315.350626979604161529144,-1015.15636749021914166146,
            -3107.77167157231109440444,22538.1184209801510330112,
            4755.84627752788110767815,-134659.959864969306392456,
            -115132.259675553483497211 };
    static doublereal c__[7] = { -.001910444077728,8.4171387781295e-4,
            -5.952379913043012e-4,7.93650793500350248e-4,
            -.002777777777777681622553,.08333333333333333331554247,
            .0057083835261 };

    /* System generated locals */
    integer i__1;
    doublereal ret_val;

    /* Builtin functions */
    doublereal sin(), log(), exp();

    /* Local variables */
    static doublereal fact, xden, xnum;
    static integer i__, j, n;
    static doublereal y, z__, y1;
    static logical parity;
    static doublereal res, sum, ysq;

/* S    REAL FUNCTION GAMMA(X) */
/* ----------------------------------------------------------------------
*/

/* THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X. */
/*     COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN W. J. CODY, */
/*     'AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL FUNCTIONS', */
/*     LECTURE NOTES IN MATHEMATICS, 506, NUMERICAL ANALYSIS DUNDEE, */
/*     1975, G. A. WATSON (ED.), SPRINGER VERLAG, BERLIN, 1976.  THE */
/*     PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA */
/*     FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS.  COEFFICIENTS
*/
/*     FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. */
/*     THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM HART, ET. AL.,
*/
/*     COMPUTER APPROXIMATIONS, WILEY AND SONS, NEW YORK, 1968.  LOWER */
/*     ORDER APPROXIMATIONS CAN BE SUBSTITUTED FOR THESE ON MACHINES */
/*     WITH LESS PRECISE ARITHMETIC. */


/* ******************************************************************* */
/* ******************************************************************* */

/* EXPLANATION OF MACHINE-DEPENDENT CONSTANTS */

/* EPS    - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT */
/*          1.0 + EPS .GT. 1.0 */
/* XBIG   - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE */
/*          IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION */
/*                  GAMMA(XBIG) = XINF. */
/* XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT */
/*          1/XMININ IS MACHINE REPRESENTABLE. */
/* XINF   - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER. */

/*     APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: */

/*         IBM/195    CDC/7600  UNIVAC/1108   VAX 11/780 (UNIX) */
/*          (D.P.)  (S.P.,RNDG)    (D.P.)     (S.P.)     (D.P.) */

/* EPS     2.221D-16  3.553E-15   1.735D-18   5.961E-08  1.388D-16 */
/* XBIG    57.574     177.802     171.489     34.844     34.844 */
/* XMININ  1.382D-76  3.132E-294  1.113D-308  5.883E-39  5.883D-39 */
/* XINF    7.23D+75   1.26E+322   8.98D+307   1.70E+38   1.70D+38 */

/* ******************************************************************* */
/* ******************************************************************* */


/* ERROR RETURNS */

/*  THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR */
/*     WHEN OVERFLOW WOULD OCCUR.  THE COMPUTATION IS BELIEVED */
/*     TO BE FREE OF UNDERFLOW AND OVERFLOW. */



/* OTHER SUBPROGRAMS REQUIRED (SINGLE PRECISION VERSION) */

/*     ALOG,EXP,FLOAT,IFIX,SIN */

/* OTHER SUBPROGRAMS REQUIRED (DOUBLE PRECISION VERSION) */

/*     DBLE,DEXP,DLOG,DSIN,FLOAT,IFIX,SNGL */



/*  AUTHOR: W. J. CODY */
/*          APPLIED MATHEMATICS DIVISION */
/*          ARGONNE NATIONAL LABORATORY */
/*          ARGONNE, IL 60439 */

/*  LATEST MODIFICATION: MAY 18, 1982 */

/* ----------------------------------------------------------------------
*/
/* S    REAL C,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI, */
/* S   1     SUM,TWELVE,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO */
/* ----------------------------------------------------------------------
*/
/*  MATHEMATICAL CONSTANTS */
/* ----------------------------------------------------------------------
*/
/* S    DATA ONE,HALF,TWELVE,ZERO/1.0E0,0.5E0,12.0E0,0.0E0/ */
/* S    DATA SQRTPI/0.9189385332046727417803297E0/ */
/* S    DATA PI/3.1415926535897932384626434E0/ */
/* ----------------------------------------------------------------------
*/
/*  MACHINE DEPENDENT PARAMETERS */
/* ----------------------------------------------------------------------
*/
/* S    DATA XBIG,XMININ,EPS/34.844E0,5.883E-39,5.961E-08/ */
/* S    DATA XINF/1.7014E38/ */
/* ----------------------------------------------------------------------
*/
/*  NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX */
/*     APPROXIMATION OVER (1,2). */
/* ----------------------------------------------------------------------
*/
/* S    DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1,
 */
/* S   1       -3.79804256470945635097577E+2,6.29331155312818442661052E+2,
 */
/* S   2       8.66966202790413211295064E+2,-3.14512729688483675254357E+4,
 */
/* S   3       -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
 */
/* S    DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2,
 */
/* S   1      -1.01515636749021914166146E+3,-3.10777167157231109440444E+3,
 */
/* S   2        2.25381184209801510330112E+4,4.75584627752788110767815E+3,
 */
/* S   3      -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
 */
/* ----------------------------------------------------------------------
*/
/*  COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). */
/* ----------------------------------------------------------------------
*/
/* S    DATA C/-1.910444077728E-03,8.4171387781295E-04, */
/* S   1     -5.952379913043012E-04,7.93650793500350248E-04, */
/* S   2     -2.777777777777681622553E-03,8.333333333333333331554247E-02,
*/
/* S   3      5.7083835261E-03/ */
/* ----------------------------------------------------------------------
*/
    parity = FALSE_;
    fact = one;
    n = 0;
    y = *x;
    if (y > zero) {
        goto L10;
    }
/* ----------------------------------------------------------------------
*/
/*  ARGUMENT IS NEGATIVE */
/* ----------------------------------------------------------------------
*/
    y = -(*x);
/* S    J = IFIX(Y) */
    j = (integer) ((real) y);
/* S    RES = Y - FLOAT(J) */
    res = y - (doublereal) ((real) j);
    if (res == zero) {
        goto L100;
    }
    if (j != j / 2 << 1) {
        parity = TRUE_;
    }
/* S    FACT = -PI / SIN(PI*RES) */
    fact = -pi / sin(pi * res);
    y += one;
/* ----------------------------------------------------------------------
*/
/*  ARGUMENT IS POSITIVE */
/* ----------------------------------------------------------------------
*/
L10:
    if (y < eps) {
        goto L90;
    }
    if (y >= twelve) {
        goto L70;
    }
    y1 = y;
    if (y >= one) {
        goto L20;
    }
/* ----------------------------------------------------------------------
*/
/*  0.0 .LT. ARGUMENT .LT. 1.0 */
/* ----------------------------------------------------------------------
*/
    z__ = y;
    y += one;
    goto L30;
/* ----------------------------------------------------------------------
*/
/*  1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY */
/* ----------------------------------------------------------------------
*/
/* S210 N = IFIX(Y) - 1 */
L20:
    n = (integer) ((real) y) - 1;
/* S    Y = Y - FLOAT(N) */
    y -= (doublereal) ((real) n);
    z__ = y - one;
/* ----------------------------------------------------------------------
*/
/*  EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 */
/* ----------------------------------------------------------------------
*/
L30:
    xnum = zero;
    xden = one;
    for (i__ = 1; i__ <= 8; ++i__) {
        xnum = (xnum + p[i__ - 1]) * z__;
        xden = xden * z__ + q[i__ - 1];
/* L40: */
    }
    res = xnum / xden + one;
    if (y == y1) {
        goto L110;
    }
    if (y1 > y) {
        goto L50;
    }
/* ----------------------------------------------------------------------
*/
/*  ADJUST RESULT FOR CASE  0.0 .LT. ARGUMENT .LT. 1.0 */
/* ----------------------------------------------------------------------
*/
    res /= y1;
    goto L110;
/* ----------------------------------------------------------------------
*/
/*  ADJUST RESULT FOR CASE  2.0 .LT. ARGUMENT .LT. 12.0 */
/* ----------------------------------------------------------------------
*/
L50:
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        res *= y;
        y += one;
/* L60: */
    }
    goto L110;
/* ----------------------------------------------------------------------
*/
/*  EVALUATE FOR ARGUMENT .GE. 12.0, */
/* ----------------------------------------------------------------------
*/
L70:
    if (y > xbig) {
        goto L100;
    }
    ysq = y * y;
    sum = c__[6];
    for (i__ = 1; i__ <= 6; ++i__) {
        sum = sum / ysq + c__[i__ - 1];
/* L80: */
    }
    sum = sum / y - y + sqrtpi;
/* S    SUM = SUM + (Y-HALF)*ALOG(Y) */
    sum += (y - half) * log(y);
/* S    RES = EXP(SUM) */
    res = exp(sum);
    goto L110;
/* ----------------------------------------------------------------------
*/
/*  ARGUMENT .LT. EPS */
/* ----------------------------------------------------------------------
*/
L90:
    if (y < xminin) {
        goto L100;
    }
    res = one / y;
    goto L110;
/* ----------------------------------------------------------------------
*/
/*  RETURN FOR SINGULARITIES, EXTREME ARGUMENTS, ETC. */
/* ----------------------------------------------------------------------
*/
/* S700 GAMMA = XINF */
L100:
    ret_val = xinf;
    goto L120;
/* ----------------------------------------------------------------------
*/
/*  FINAL ADJUSTMENTS AND RETURN */
/* ----------------------------------------------------------------------
*/
L110:
    if (parity) {
        res = -res;
    }
    if (fact != one) {
        res = fact / res;
    }
/* S    GAMMA = RES */
    ret_val = res;
L120:
    return ret_val;
/* ---------- LAST CARD OF GAMMA ---------- */
} /* dgamma_ */

/* Subroutine */ integer machar_(ibeta, it, irnd, ngrd, machep, negep, iexp,
        minexp, maxexp, eps, epsneg, xmin, xmax)
integer *ibeta, *it, *irnd, *ngrd, *machep, *negep, *iexp, *minexp, *maxexp;
doublereal *eps, *epsneg, *xmin, *xmax;
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    static doublereal beta, zero, a, b;
    static integer i__, j, k;
    static doublereal y, z__, betam1, betain;
    static integer iz, mx;
    static doublereal one;


/* S    REAL A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,XMIN,Y,Z,ZERO */

/*     THIS SUBROUTINE IS INTENDED TO DETERMINE THE CHARACTERISTICS */
/*     OF THE FLOATING-POINT ARITHMETIC SYSTEM THAT ARE SPECIFIED */
/*     BELOW.  THE FIRST THREE ARE DETERMINED ACCORDING TO AN */
/*     ALGORITHM DUE TO M. MALCOLM, CACM 15 (1972), PP. 949-951, */
/*     INCORPORATING SOME, BUT NOT ALL, OF THE IMPROVEMENTS */
/*     SUGGESTED BY M. GENTLEMAN AND S. MAROVICH, CACM 17 (1974), */
/*     PP. 276-277.  THE VERSION GIVEN HERE IS FOR SINGLE PRECISION. */
/*     CARDS CONTAINING  CD  IN COLUMNS 1 AND 2 CAN BE USED TO */
/*     CONVERT THE SUBROUTINE TO DOUBLE PRECISION BY REPLACING */
/*     EXISTING CARDS IN THE OBVIOUS MANNER. */


/*       IBETA   - THE RADIX OF THE FLOATING-POINT REPRESENTATION */
/*       IT      - THE NUMBER OF BASE IBETA DIGITS IN THE FLOATING-POINT
*/
/*                 SIGNIFICAND */
/*       IRND    - 0 IF FLOATING-POINT ADDITION CHOPS, */
/*                 1 IF FLOATING-POINT ADDITION ROUNDS */
/*       NGRD    - THE NUMBER OF GUARD DIGITS FOR MULTIPLICATION.  IT IS
*/
/*                 0 IF  IRND=1, OR IF  IRND=0  AND ONLY  IT  BASE  IBETA
*/
/*                   DIGITS PARTICIPATE IN THE POST NORMALIZATION SHIFT */
/*                   OF THE FLOATING-POINT SIGNIFICAND IN MULTIPLICATION
*/
/*                 1 IF  IRND=0  AND MORE THAN  IT  BASE  IBETA  DIGITS */
/*                   PARTICIPATE IN THE POST NORMALIZATION SHIFT OF THE */
/*                   FLOATING-POINT SIGNIFICAND IN MULTIPLICATION */
/*       MACHEP  - THE LARGEST NEGATIVE INTEGER SUCH THAT */
/*                 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, EXCEPT THAT */
/*                 MACHEP IS BOUNDED BELOW BY  -(IT+3) */
/*       NEGEPS  - THE LARGEST NEGATIVE INTEGER SUCH THAT */
/*                 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, EXCEPT THAT */
/*                 NEGEPS IS BOUNDED BELOW BY  -(IT+3) */
/*       IEXP    - THE NUMBER OF BITS (DECIMAL PLACES IF IBETA = 10) */
/*                 RESERVED FOR THE REPRESENTATION OF THE EXPONENT */
/*                 (INCLUDING THE BIAS OR SIGN) OF A FLOATING-POINT */
/*                 NUMBER */
/*       MINEXP  - THE LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH THAT */
/*                 FLOAT(IBETA)**MINEXP IS A POSITIVE FLOATING-POINT */
/*                 NUMBER */
/*       MAXEXP  - THE LARGEST POSITIVE INTEGER EXPONENT FOR A FINITE */
/*                 FLOATING-POINT NUMBER */
/*       EPS     - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH */
/*                 THAT  1.0+EPS .NE. 1.0. IN PARTICULAR, IF EITHER */
/*                 IBETA = 2  OR  IRND = 0, EPS = FLOAT(IBETA)**MACHEP. */
/*                 OTHERWISE,  EPS = (FLOAT(IBETA)**MACHEP)/2 */
/*       EPSNEG  - A SMALL POSITIVE FLOATING-POINT NUMBER SUCH THAT */
/*                 1.0-EPSNEG .NE. 1.0. IN PARTICULAR, IF IBETA = 2 */
/*                 OR  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS. */
/*                 OTHERWISE,  EPSNEG = (IBETA**NEGEPS)/2.  BECAUSE */
/*                 NEGEPS IS BOUNDED BELOW BY -(IT+3), EPSNEG MAY NOT */
/*                 BE THE SMALLEST NUMBER WHICH CAN ALTER 1.0 BY */
/*                 SUBTRACTION. */
/*       XMIN    - THE SMALLEST NON-VANISHING FLOATING-POINT POWER OF THE
*/
/*                 RADIX.  IN PARTICULAR,  XMIN = FLOAT(IBETA)**MINEXP */
/*       XMAX    - THE LARGEST FINITE FLOATING-POINT NUMBER.  IN */
/*                 PARTICULAR   XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
*/
/*                 NOTE - ON SOME MACHINES  XMAX  WILL BE ONLY THE */
/*                 SECOND, OR PERHAPS THIRD, LARGEST NUMBER, BEING */
/*                 TOO SMALL BY 1 OR 2 UNITS IN THE LAST DIGIT OF */
/*                 THE SIGNIFICAND. */

/*     LATEST REVISION - OCTOBER 22, 1979 */

/*     AUTHOR - W. J. CODY */
/*              ARGONNE NATIONAL LABORATORY */

/* ----------------------------------------------------------------- */
/* S    ONE = FLOAT(1) */
    one = 1.;
/* S    ZERO = 0.0E0 */
    zero = 0.;
/* ----------------------------------------------------------------- */
/*     DETERMINE IBETA,BETA ALA MALCOLM */
/* ----------------------------------------------------------------- */
    a = one;
L10:
    a += a;
    if (a + one - a - one == zero) {
        goto L10;
    }
    b = one;
L20:
    b += b;
    if (a + b - a == zero) {
        goto L20;
    }
/* S    IBETA = INT((A+B)-A) */
    *ibeta = (integer) ((real) (a + b - a));
/* S    BETA = FLOAT(IBETA) */
    beta = (doublereal) ((real) (*ibeta));
/* ----------------------------------------------------------------- */
/*     DETERMINE IT, IRND */
/* ----------------------------------------------------------------- */
    *it = 0;
    b = one;
L30:
    ++(*it);
    b *= beta;
    if (b + one - b - one == zero) {
        goto L30;
    }
    *irnd = 0;
    betam1 = beta - one;
    if (a + betam1 - a != zero) {
        *irnd = 1;
    }
/* ----------------------------------------------------------------- */
/*     DETERMINE NEGEP, EPSNEG */
/* ----------------------------------------------------------------- */
    *negep = *it + 3;
    betain = one / beta;
    a = one;

    i__1 = *negep;
    for (i__ = 1; i__ <= i__1; ++i__) {
        a *= betain;
/* L40: */
    }

    b = a;
L50:
    if (one - a - one != zero) {
        goto L60;
    }
    a *= beta;
    --(*negep);
    goto L50;
L60:
    *negep = -(*negep);
    *epsneg = a;
    if (*ibeta == 2 || *irnd == 0) {
        goto L70;
    }
    a = a * (one + a) / (one + one);
    if (one - a - one != zero) {
        *epsneg = a;
    }
/* ----------------------------------------------------------------- */
/*     DETERMINE MACHEP, EPS */
/* ----------------------------------------------------------------- */
L70:
    *machep = -(*it) - 3;
    a = b;
L80:
    if (one + a - one != zero) {
        goto L90;
    }
    a *= beta;
    ++(*machep);
    goto L80;
L90:
    *eps = a;
    if (*ibeta == 2 || *irnd == 0) {
        goto L100;
    }
    a = a * (one + a) / (one + one);
    if (one + a - one != zero) {
        *eps = a;
    }
/* ----------------------------------------------------------------- */
/*     DETERMINE NGRD */
/* ----------------------------------------------------------------- */
L100:
    *ngrd = 0;
    if (*irnd == 0 && (one + *eps) * one - one != zero) {
        *ngrd = 1;
    }
/* ----------------------------------------------------------------- */
/*     DETERMINE IEXP, MINEXP, XMIN */

/*     LOOP TO DETERMINE LARGEST I AND K = 2**I SUCH THAT */
/*         (1/BETA) ** (2**(I)) */
/*     DOES NOT UNDERFLOW */
/*     EXIT FROM LOOP IS SIGNALED BY AN UNDERFLOW. */
/* ----------------------------------------------------------------- */
    i__ = 0;
    k = 1;
    z__ = betain;
L110:
    y = z__;
    z__ = y * y;
/* ----------------------------------------------------------------- */
/*        CHECK FOR UNDERFLOW HERE */
/* ----------------------------------------------------------------- */
    a = z__ * one;
/* S       IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410 */
    if (a + a == zero || abs(z__) >= y) {
        goto L120;
    }
    ++i__;
    k += k;
    goto L110;
L120:
    if (*ibeta == 10) {
        goto L130;
    }
    *iexp = i__ + 1;
    mx = k + k;
    goto L160;
/* ----------------------------------------------------------------- */
/*     FOR DECIMAL MACHINES ONLY */
/* ----------------------------------------------------------------- */
L130:
    *iexp = 2;
    iz = *ibeta;
L140:
    if (k < iz) {
        goto L150;
    }
    iz *= *ibeta;
    ++(*iexp);
    goto L140;
L150:
    mx = iz + iz - 1;
/* ----------------------------------------------------------------- */
/*     LOOP TO DETERMINE MINEXP, XMIN */
/*     EXIT FROM LOOP IS SIGNALED BY AN UNDERFLOW. */
/* ----------------------------------------------------------------- */
L160:
    *xmin = y;
    y *= betain;
/* ----------------------------------------------------------------- */
/*        CHECK FOR UNDERFLOW HERE */
/* ----------------------------------------------------------------- */
    a = y * one;
/* S       IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460 */
    if (a + a == zero || abs(y) >= *xmin) {
        goto L170;
    }
    ++k;
    goto L160;
L170:
    *minexp = -k;
/* ----------------------------------------------------------------- */
/*     DETERMINE MAXEXP, XMAX */
/* ----------------------------------------------------------------- */
    if (mx > k + k - 3 || *ibeta == 10) {
        goto L180;
    }
    mx += mx;
    ++(*iexp);
L180:
    *maxexp = mx + *minexp;
/* ----------------------------------------------------------------- */
/*     ADJUST FOR MACHINES WITH IMPLICIT LEADING */
/*     BIT IN BINARY SIGNIFICAND AND MACHINES WITH */
/*     RADIX POINT AT EXTREME RIGHT OF SIGNIFICAND */
/* ----------------------------------------------------------------- */
    i__ = *maxexp + *minexp;
    if (*ibeta == 2 && i__ == 0) {
        --(*maxexp);
    }
    if (i__ > 20) {
        --(*maxexp);
    }
    if (a != y) {
        *maxexp += -2;
    }
    *xmax = one - *epsneg;
    if (*xmax * one != *xmax) {
        *xmax = one - beta * *epsneg;
    }
    *xmax /= beta * beta * beta * *xmin;
    i__ = *maxexp + *minexp + 3;
    if (i__ <= 0) {
        goto L200;
    }

    i__1 = i__;
    for (j = 1; j <= i__1; ++j) {
        if (*ibeta == 2) {
            *xmax += *xmax;
        }
        if (*ibeta != 2) {
            *xmax *= beta;
        }
/* L190: */
    }

L200:
    return 0;
/*     ---------- LAST CARD OF MACHAR ---------- */
} /* machar_ */

doublereal ren_(k)
integer *k;
{
    /* Initialized data */

    static integer iy = 100001;

    /* System generated locals */
    doublereal ret_val;

    /* Local variables */
    static integer j;


/*     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND */
/*      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM, */
/*      VOL. 8, NO. 10, OCTOBER 1965. */

/*     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH */
/*      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS */
/*      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST */
/*      29 BITS. */


    j = *k;
    iy *= 125;
    iy -= iy / 2796203 * 2796203;
    ret_val = (doublereal) ((real) iy) / 2796203. * 1.000001000001;
    return ret_val;
/*     ---------- LAST CARD OF REN ---------- */
} /* ren_ */

#ifdef uNdEfInEd
comments from the converter:  (stderr from f2c)
   ribesl:
   dgamma:
   machar:
   ren:
#endif