/* src/slamch.f -- translated by f2c (version 20050501). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib; on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */ #include "sphinxbase/f2c.h" #ifdef _MSC_VER #pragma warning (disable: 4244) #endif /* Table of constant values */ static integer c__1 = 1; static real c_b32 = 0.f; doublereal slamch_(char *cmach, ftnlen cmach_len) { /* Initialized data */ static logical first = TRUE_; /* System generated locals */ integer i__1; real ret_val; /* Builtin functions */ double pow_ri(real *, integer *); /* Local variables */ static real t; static integer it; static real rnd, eps, base; static integer beta; static real emin, prec, emax; static integer imin, imax; static logical lrnd; static real rmin, rmax, rmach; extern logical lsame_(char *, char *, ftnlen, ftnlen); static real small, sfmin; extern /* Subroutine */ int slamc2_(integer *, integer *, logical *, real *, integer *, real *, integer *, real *); /* -- LAPACK auxiliary routine (version 3.0) -- */ /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ /* Courant Institute, Argonne National Lab, and Rice University */ /* October 31, 1992 */ /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* SLAMCH determines single precision machine parameters. */ /* Arguments */ /* ========= */ /* CMACH (input) CHARACTER*1 */ /* Specifies the value to be returned by SLAMCH: */ /* = 'E' or 'e', SLAMCH := eps */ /* = 'S' or 's , SLAMCH := sfmin */ /* = 'B' or 'b', SLAMCH := base */ /* = 'P' or 'p', SLAMCH := eps*base */ /* = 'N' or 'n', SLAMCH := t */ /* = 'R' or 'r', SLAMCH := rnd */ /* = 'M' or 'm', SLAMCH := emin */ /* = 'U' or 'u', SLAMCH := rmin */ /* = 'L' or 'l', SLAMCH := emax */ /* = 'O' or 'o', SLAMCH := rmax */ /* where */ /* eps = relative machine precision */ /* sfmin = safe minimum, such that 1/sfmin does not overflow */ /* base = base of the machine */ /* prec = eps*base */ /* t = number of (base) digits in the mantissa */ /* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */ /* emin = minimum exponent before (gradual) underflow */ /* rmin = underflow threshold - base**(emin-1) */ /* emax = largest exponent before overflow */ /* rmax = overflow threshold - (base**emax)*(1-eps) */ /* ===================================================================== */ /* .. Parameters .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. External Functions .. */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Save statement .. */ /* .. */ /* .. Data statements .. */ /* .. */ /* .. Executable Statements .. */ if (first) { first = FALSE_; slamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax); base = (real) beta; t = (real) it; if (lrnd) { rnd = 1.f; i__1 = 1 - it; eps = pow_ri(&base, &i__1) / 2; } else { rnd = 0.f; i__1 = 1 - it; eps = pow_ri(&base, &i__1); } prec = eps * base; emin = (real) imin; emax = (real) imax; sfmin = rmin; small = 1.f / rmax; if (small >= sfmin) { /* Use SMALL plus a bit, to avoid the possibility of rounding */ /* causing overflow when computing 1/sfmin. */ sfmin = small * (eps + 1.f); } } if (lsame_(cmach, "E", (ftnlen) 1, (ftnlen) 1)) { rmach = eps; } else if (lsame_(cmach, "S", (ftnlen) 1, (ftnlen) 1)) { rmach = sfmin; } else if (lsame_(cmach, "B", (ftnlen) 1, (ftnlen) 1)) { rmach = base; } else if (lsame_(cmach, "P", (ftnlen) 1, (ftnlen) 1)) { rmach = prec; } else if (lsame_(cmach, "N", (ftnlen) 1, (ftnlen) 1)) { rmach = t; } else if (lsame_(cmach, "R", (ftnlen) 1, (ftnlen) 1)) { rmach = rnd; } else if (lsame_(cmach, "M", (ftnlen) 1, (ftnlen) 1)) { rmach = emin; } else if (lsame_(cmach, "U", (ftnlen) 1, (ftnlen) 1)) { rmach = rmin; } else if (lsame_(cmach, "L", (ftnlen) 1, (ftnlen) 1)) { rmach = emax; } else if (lsame_(cmach, "O", (ftnlen) 1, (ftnlen) 1)) { rmach = rmax; } ret_val = rmach; return ret_val; /* End of SLAMCH */ } /* slamch_ */ /* *********************************************************************** */ /* Subroutine */ int slamc1_(integer * beta, integer * t, logical * rnd, logical * ieee1) { /* Initialized data */ static logical first = TRUE_; /* System generated locals */ real r__1, r__2; /* Local variables */ static real a, b, c__, f, t1, t2; static integer lt; static real one, qtr; static logical lrnd; static integer lbeta; static real savec; static logical lieee1; extern doublereal slamc3_(real *, real *); /* -- LAPACK auxiliary routine (version 3.0) -- */ /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ /* Courant Institute, Argonne National Lab, and Rice University */ /* October 31, 1992 */ /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* SLAMC1 determines the machine parameters given by BETA, T, RND, and */ /* IEEE1. */ /* Arguments */ /* ========= */ /* BETA (output) INTEGER */ /* The base of the machine. */ /* T (output) INTEGER */ /* The number of ( BETA ) digits in the mantissa. */ /* RND (output) LOGICAL */ /* Specifies whether proper rounding ( RND = .TRUE. ) or */ /* chopping ( RND = .FALSE. ) occurs in addition. This may not */ /* be a reliable guide to the way in which the machine performs */ /* its arithmetic. */ /* IEEE1 (output) LOGICAL */ /* Specifies whether rounding appears to be done in the IEEE */ /* 'round to nearest' style. */ /* Further Details */ /* =============== */ /* The routine is based on the routine ENVRON by Malcolm and */ /* incorporates suggestions by Gentleman and Marovich. See */ /* Malcolm M. A. (1972) Algorithms to reveal properties of */ /* floating-point arithmetic. Comms. of the ACM, 15, 949-951. */ /* Gentleman W. M. and Marovich S. B. (1974) More on algorithms */ /* that reveal properties of floating point arithmetic units. */ /* Comms. of the ACM, 17, 276-277. */ /* ===================================================================== */ /* .. Local Scalars .. */ /* .. */ /* .. External Functions .. */ /* .. */ /* .. Save statement .. */ /* .. */ /* .. Data statements .. */ /* .. */ /* .. Executable Statements .. */ if (first) { first = FALSE_; one = 1.f; /* LBETA, LIEEE1, LT and LRND are the local values of BETA, */ /* IEEE1, T and RND. */ /* Throughout this routine we use the function SLAMC3 to ensure */ /* that relevant values are stored and not held in registers, or */ /* are not affected by optimizers. */ /* Compute a = 2.0**m with the smallest positive integer m such */ /* that */ /* fl( a + 1.0 ) = a. */ a = 1.f; c__ = 1.f; /* + WHILE( C.EQ.ONE )LOOP */ L10: if (c__ == one) { a *= 2; c__ = slamc3_(&a, &one); r__1 = -a; c__ = slamc3_(&c__, &r__1); goto L10; } /* + END WHILE */ /* Now compute b = 2.0**m with the smallest positive integer m */ /* such that */ /* fl( a + b ) .gt. a. */ b = 1.f; c__ = slamc3_(&a, &b); /* + WHILE( C.EQ.A )LOOP */ L20: if (c__ == a) { b *= 2; c__ = slamc3_(&a, &b); goto L20; } /* + END WHILE */ /* Now compute the base. a and c are neighbouring floating point */ /* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so */ /* their difference is beta. Adding 0.25 to c is to ensure that it */ /* is truncated to beta and not ( beta - 1 ). */ qtr = one / 4; savec = c__; r__1 = -a; c__ = slamc3_(&c__, &r__1); lbeta = c__ + qtr; /* Now determine whether rounding or chopping occurs, by adding a */ /* bit less than beta/2 and a bit more than beta/2 to a. */ b = (real) lbeta; r__1 = b / 2; r__2 = -b / 100; f = slamc3_(&r__1, &r__2); c__ = slamc3_(&f, &a); if (c__ == a) { lrnd = TRUE_; } else { lrnd = FALSE_; } r__1 = b / 2; r__2 = b / 100; f = slamc3_(&r__1, &r__2); c__ = slamc3_(&f, &a); if (lrnd && c__ == a) { lrnd = FALSE_; } /* Try and decide whether rounding is done in the IEEE 'round to */ /* nearest' style. B/2 is half a unit in the last place of the two */ /* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */ /* zero, and SAVEC is odd. Thus adding B/2 to A should not change */ /* A, but adding B/2 to SAVEC should change SAVEC. */ r__1 = b / 2; t1 = slamc3_(&r__1, &a); r__1 = b / 2; t2 = slamc3_(&r__1, &savec); lieee1 = t1 == a && t2 > savec && lrnd; /* Now find the mantissa, t. It should be the integer part of */ /* log to the base beta of a, however it is safer to determine t */ /* by powering. So we find t as the smallest positive integer for */ /* which */ /* fl( beta**t + 1.0 ) = 1.0. */ lt = 0; a = 1.f; c__ = 1.f; /* + WHILE( C.EQ.ONE )LOOP */ L30: if (c__ == one) { ++lt; a *= lbeta; c__ = slamc3_(&a, &one); r__1 = -a; c__ = slamc3_(&c__, &r__1); goto L30; } /* + END WHILE */ } *beta = lbeta; *t = lt; *rnd = lrnd; *ieee1 = lieee1; return 0; /* End of SLAMC1 */ } /* slamc1_ */ /* *********************************************************************** */ /* Subroutine */ int slamc2_(integer * beta, integer * t, logical * rnd, real * eps, integer * emin, real * rmin, integer * emax, real * rmax) { /* Initialized data */ static logical first = TRUE_; static logical iwarn = FALSE_; /* Format strings */ static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre" "ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va" "lue EMIN looks\002,\002 acceptable please comment out \002,/\002" " the IF block as marked within the code of routine\002,\002 SLAM" "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)"; /* System generated locals */ integer i__1; real r__1, r__2, r__3, r__4, r__5; /* Builtin functions */ double pow_ri(real *, integer *); integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static real a, b, c__; static integer i__, lt; static real one, two; static logical ieee; static real half; static logical lrnd; static real leps, zero; static integer lbeta; static real rbase; static integer lemin, lemax, gnmin; static real small; static integer gpmin; static real third, lrmin, lrmax, sixth; static logical lieee1; extern /* Subroutine */ int slamc1_(integer *, integer *, logical *, logical *); extern doublereal slamc3_(real *, real *); extern /* Subroutine */ int slamc4_(integer *, real *, integer *), slamc5_(integer *, integer *, integer *, logical *, integer *, real *); static integer ngnmin, ngpmin; /* Fortran I/O blocks */ static cilist io___58 = { 0, 6, 0, fmt_9999, 0 }; /* -- LAPACK auxiliary routine (version 3.0) -- */ /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ /* Courant Institute, Argonne National Lab, and Rice University */ /* October 31, 1992 */ /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* SLAMC2 determines the machine parameters specified in its argument */ /* list. */ /* Arguments */ /* ========= */ /* BETA (output) INTEGER */ /* The base of the machine. */ /* T (output) INTEGER */ /* The number of ( BETA ) digits in the mantissa. */ /* RND (output) LOGICAL */ /* Specifies whether proper rounding ( RND = .TRUE. ) or */ /* chopping ( RND = .FALSE. ) occurs in addition. This may not */ /* be a reliable guide to the way in which the machine performs */ /* its arithmetic. */ /* EPS (output) REAL */ /* The smallest positive number such that */ /* fl( 1.0 - EPS ) .LT. 1.0, */ /* where fl denotes the computed value. */ /* EMIN (output) INTEGER */ /* The minimum exponent before (gradual) underflow occurs. */ /* RMIN (output) REAL */ /* The smallest normalized number for the machine, given by */ /* BASE**( EMIN - 1 ), where BASE is the floating point value */ /* of BETA. */ /* EMAX (output) INTEGER */ /* The maximum exponent before overflow occurs. */ /* RMAX (output) REAL */ /* The largest positive number for the machine, given by */ /* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */ /* value of BETA. */ /* Further Details */ /* =============== */ /* The computation of EPS is based on a routine PARANOIA by */ /* W. Kahan of the University of California at Berkeley. */ /* ===================================================================== */ /* .. Local Scalars .. */ /* .. */ /* .. External Functions .. */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Save statement .. */ /* .. */ /* .. Data statements .. */ /* .. */ /* .. Executable Statements .. */ if (first) { first = FALSE_; zero = 0.f; one = 1.f; two = 2.f; /* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */ /* BETA, T, RND, EPS, EMIN and RMIN. */ /* Throughout this routine we use the function SLAMC3 to ensure */ /* that relevant values are stored and not held in registers, or */ /* are not affected by optimizers. */ /* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */ slamc1_(&lbeta, <, &lrnd, &lieee1); /* Start to find EPS. */ b = (real) lbeta; i__1 = -lt; a = pow_ri(&b, &i__1); leps = a; /* Try some tricks to see whether or not this is the correct EPS. */ b = two / 3; half = one / 2; r__1 = -half; sixth = slamc3_(&b, &r__1); third = slamc3_(&sixth, &sixth); r__1 = -half; b = slamc3_(&third, &r__1); b = slamc3_(&b, &sixth); b = dabs(b); if (b < leps) { b = leps; } leps = 1.f; /* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */ L10: if (leps > b && b > zero) { leps = b; r__1 = half * leps; /* Computing 5th power */ r__3 = two, r__4 = r__3, r__3 *= r__3; /* Computing 2nd power */ r__5 = leps; r__2 = r__4 * (r__3 * r__3) * (r__5 * r__5); c__ = slamc3_(&r__1, &r__2); r__1 = -c__; c__ = slamc3_(&half, &r__1); b = slamc3_(&half, &c__); r__1 = -b; c__ = slamc3_(&half, &r__1); b = slamc3_(&half, &c__); goto L10; } /* + END WHILE */ if (a < leps) { leps = a; } /* Computation of EPS complete. */ /* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */ /* Keep dividing A by BETA until (gradual) underflow occurs. This */ /* is detected when we cannot recover the previous A. */ rbase = one / lbeta; small = one; for (i__ = 1; i__ <= 3; ++i__) { r__1 = small * rbase; small = slamc3_(&r__1, &zero); /* L20: */ } a = slamc3_(&one, &small); slamc4_(&ngpmin, &one, &lbeta); r__1 = -one; slamc4_(&ngnmin, &r__1, &lbeta); slamc4_(&gpmin, &a, &lbeta); r__1 = -a; slamc4_(&gnmin, &r__1, &lbeta); ieee = FALSE_; if (ngpmin == ngnmin && gpmin == gnmin) { if (ngpmin == gpmin) { lemin = ngpmin; /* ( Non twos-complement machines, no gradual underflow; */ /* e.g., VAX ) */ } else if (gpmin - ngpmin == 3) { lemin = ngpmin - 1 + lt; ieee = TRUE_; /* ( Non twos-complement machines, with gradual underflow; */ /* e.g., IEEE standard followers ) */ } else { lemin = min(ngpmin, gpmin); /* ( A guess; no known machine ) */ iwarn = TRUE_; } } else if (ngpmin == gpmin && ngnmin == gnmin) { if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) { lemin = max(ngpmin, ngnmin); /* ( Twos-complement machines, no gradual underflow; */ /* e.g., CYBER 205 ) */ } else { lemin = min(ngpmin, ngnmin); /* ( A guess; no known machine ) */ iwarn = TRUE_; } } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin) { if (gpmin - min(ngpmin, ngnmin) == 3) { lemin = max(ngpmin, ngnmin) - 1 + lt; /* ( Twos-complement machines with gradual underflow; */ /* no known machine ) */ } else { lemin = min(ngpmin, ngnmin); /* ( A guess; no known machine ) */ iwarn = TRUE_; } } else { /* Computing MIN */ i__1 = min(ngpmin, ngnmin), i__1 = min(i__1, gpmin); lemin = min(i__1, gnmin); /* ( A guess; no known machine ) */ iwarn = TRUE_; } /* ** */ /* Comment out this if block if EMIN is ok */ if (iwarn) { first = TRUE_; s_wsfe(&io___58); do_fio(&c__1, (char *) &lemin, (ftnlen) sizeof(integer)); e_wsfe(); } /* ** */ /* Assume IEEE arithmetic if we found denormalised numbers above, */ /* or if arithmetic seems to round in the IEEE style, determined */ /* in routine SLAMC1. A true IEEE machine should have both things */ /* true; however, faulty machines may have one or the other. */ ieee = ieee || lieee1; /* Compute RMIN by successive division by BETA. We could compute */ /* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */ /* this computation. */ lrmin = 1.f; i__1 = 1 - lemin; for (i__ = 1; i__ <= i__1; ++i__) { r__1 = lrmin * rbase; lrmin = slamc3_(&r__1, &zero); /* L30: */ } /* Finally, call SLAMC5 to compute EMAX and RMAX. */ slamc5_(&lbeta, <, &lemin, &ieee, &lemax, &lrmax); } *beta = lbeta; *t = lt; *rnd = lrnd; *eps = leps; *emin = lemin; *rmin = lrmin; *emax = lemax; *rmax = lrmax; return 0; /* End of SLAMC2 */ } /* slamc2_ */ /* *********************************************************************** */ doublereal slamc3_(real * a, real * b) { /* System generated locals */ real ret_val; /* -- LAPACK auxiliary routine (version 3.0) -- */ /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ /* Courant Institute, Argonne National Lab, and Rice University */ /* October 31, 1992 */ /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* SLAMC3 is intended to force A and B to be stored prior to doing */ /* the addition of A and B , for use in situations where optimizers */ /* might hold one of these in a register. */ /* Arguments */ /* ========= */ /* A, B (input) REAL */ /* The values A and B. */ /* ===================================================================== */ /* .. Executable Statements .. */ ret_val = *a + *b; return ret_val; /* End of SLAMC3 */ } /* slamc3_ */ /* *********************************************************************** */ /* Subroutine */ int slamc4_(integer * emin, real * start, integer * base) { /* System generated locals */ integer i__1; real r__1; /* Local variables */ static real a; static integer i__; static real b1, b2, c1, c2, d1, d2, one, zero, rbase; extern doublereal slamc3_(real *, real *); /* -- LAPACK auxiliary routine (version 3.0) -- */ /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ /* Courant Institute, Argonne National Lab, and Rice University */ /* October 31, 1992 */ /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* SLAMC4 is a service routine for SLAMC2. */ /* Arguments */ /* ========= */ /* EMIN (output) EMIN */ /* The minimum exponent before (gradual) underflow, computed by */ /* setting A = START and dividing by BASE until the previous A */ /* can not be recovered. */ /* START (input) REAL */ /* The starting point for determining EMIN. */ /* BASE (input) INTEGER */ /* The base of the machine. */ /* ===================================================================== */ /* .. Local Scalars .. */ /* .. */ /* .. External Functions .. */ /* .. */ /* .. Executable Statements .. */ a = *start; one = 1.f; rbase = one / *base; zero = 0.f; *emin = 1; r__1 = a * rbase; b1 = slamc3_(&r__1, &zero); c1 = a; c2 = a; d1 = a; d2 = a; /* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */ /* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */ L10: if (c1 == a && c2 == a && d1 == a && d2 == a) { --(*emin); a = b1; r__1 = a / *base; b1 = slamc3_(&r__1, &zero); r__1 = b1 * *base; c1 = slamc3_(&r__1, &zero); d1 = zero; i__1 = *base; for (i__ = 1; i__ <= i__1; ++i__) { d1 += b1; /* L20: */ } r__1 = a * rbase; b2 = slamc3_(&r__1, &zero); r__1 = b2 / rbase; c2 = slamc3_(&r__1, &zero); d2 = zero; i__1 = *base; for (i__ = 1; i__ <= i__1; ++i__) { d2 += b2; /* L30: */ } goto L10; } /* + END WHILE */ return 0; /* End of SLAMC4 */ } /* slamc4_ */ /* *********************************************************************** */ /* Subroutine */ int slamc5_(integer * beta, integer * p, integer * emin, logical * ieee, integer * emax, real * rmax) { /* System generated locals */ integer i__1; real r__1; /* Local variables */ static integer i__; static real y, z__; static integer try__, lexp; static real oldy; static integer uexp, nbits; extern doublereal slamc3_(real *, real *); static real recbas; static integer exbits, expsum; /* -- LAPACK auxiliary routine (version 3.0) -- */ /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ /* Courant Institute, Argonne National Lab, and Rice University */ /* October 31, 1992 */ /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* SLAMC5 attempts to compute RMAX, the largest machine floating-point */ /* number, without overflow. It assumes that EMAX + abs(EMIN) sum */ /* approximately to a power of 2. It will fail on machines where this */ /* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */ /* EMAX = 28718). It will also fail if the value supplied for EMIN is */ /* too large (i.e. too close to zero), probably with overflow. */ /* Arguments */ /* ========= */ /* BETA (input) INTEGER */ /* The base of floating-point arithmetic. */ /* P (input) INTEGER */ /* The number of base BETA digits in the mantissa of a */ /* floating-point value. */ /* EMIN (input) INTEGER */ /* The minimum exponent before (gradual) underflow. */ /* IEEE (input) LOGICAL */ /* A logical flag specifying whether or not the arithmetic */ /* system is thought to comply with the IEEE standard. */ /* EMAX (output) INTEGER */ /* The largest exponent before overflow */ /* RMAX (output) REAL */ /* The largest machine floating-point number. */ /* ===================================================================== */ /* .. Parameters .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. External Functions .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* First compute LEXP and UEXP, two powers of 2 that bound */ /* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */ /* approximately to the bound that is closest to abs(EMIN). */ /* (EMAX is the exponent of the required number RMAX). */ lexp = 1; exbits = 1; L10: try__ = lexp << 1; if (try__ <= -(*emin)) { lexp = try__; ++exbits; goto L10; } if (lexp == -(*emin)) { uexp = lexp; } else { uexp = try__; ++exbits; } /* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */ /* than or equal to EMIN. EXBITS is the number of bits needed to */ /* store the exponent. */ if (uexp + *emin > -lexp - *emin) { expsum = lexp << 1; } else { expsum = uexp << 1; } /* EXPSUM is the exponent range, approximately equal to */ /* EMAX - EMIN + 1 . */ *emax = expsum + *emin - 1; nbits = exbits + 1 + *p; /* NBITS is the total number of bits needed to store a */ /* floating-point number. */ if (nbits % 2 == 1 && *beta == 2) { /* Either there are an odd number of bits used to store a */ /* floating-point number, which is unlikely, or some bits are */ /* not used in the representation of numbers, which is possible, */ /* (e.g. Cray machines) or the mantissa has an implicit bit, */ /* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */ /* most likely. We have to assume the last alternative. */ /* If this is true, then we need to reduce EMAX by one because */ /* there must be some way of representing zero in an implicit-bit */ /* system. On machines like Cray, we are reducing EMAX by one */ /* unnecessarily. */ --(*emax); } if (*ieee) { /* Assume we are on an IEEE machine which reserves one exponent */ /* for infinity and NaN. */ --(*emax); } /* Now create RMAX, the largest machine number, which should */ /* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */ /* First compute 1.0 - BETA**(-P), being careful that the */ /* result is less than 1.0 . */ recbas = 1.f / *beta; z__ = *beta - 1.f; y = 0.f; i__1 = *p; for (i__ = 1; i__ <= i__1; ++i__) { z__ *= recbas; if (y < 1.f) { oldy = y; } y = slamc3_(&y, &z__); /* L20: */ } if (y >= 1.f) { y = oldy; } /* Now multiply by BETA**EMAX to get RMAX. */ i__1 = *emax; for (i__ = 1; i__ <= i__1; ++i__) { r__1 = y * *beta; y = slamc3_(&r__1, &c_b32); /* L30: */ } *rmax = y; return 0; /* End of SLAMC5 */ } /* slamc5_ */