diff options
Diffstat (limited to 'media/sphinxbase/src/libsphinxbase/util/slamch.c')
| -rw-r--r-- | media/sphinxbase/src/libsphinxbase/util/slamch.c | 1029 |
1 files changed, 0 insertions, 1029 deletions
diff --git a/media/sphinxbase/src/libsphinxbase/util/slamch.c b/media/sphinxbase/src/libsphinxbase/util/slamch.c deleted file mode 100644 index 229458470..000000000 --- a/media/sphinxbase/src/libsphinxbase/util/slamch.c +++ /dev/null @@ -1,1029 +0,0 @@ -/* 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_ */ |