• Main Page
  • Related Pages
  • Data Structures
  • Files
  • File List
  • Globals

src/libsphinxbase/util/slamch.c

00001 /* src/slamch.f -- translated by f2c (version 20050501).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 
00015 #ifdef _MSC_VER
00016 #pragma warning (disable: 4244)
00017 #endif
00018 
00019 /* Table of constant values */
00020 
00021 static integer c__1 = 1;
00022 static real c_b32 = 0.f;
00023 
00024 doublereal
00025 slamch_(char *cmach, ftnlen cmach_len)
00026 {
00027     /* Initialized data */
00028 
00029     static logical first = TRUE_;
00030 
00031     /* System generated locals */
00032     integer i__1;
00033     real ret_val;
00034 
00035     /* Builtin functions */
00036     double pow_ri(real *, integer *);
00037 
00038     /* Local variables */
00039     static real t;
00040     static integer it;
00041     static real rnd, eps, base;
00042     static integer beta;
00043     static real emin, prec, emax;
00044     static integer imin, imax;
00045     static logical lrnd;
00046     static real rmin, rmax, rmach;
00047     extern logical lsame_(char *, char *, ftnlen, ftnlen);
00048     static real small, sfmin;
00049     extern /* Subroutine */ int slamc2_(integer *, integer *, logical *, real
00050                                         *, integer *, real *, integer *,
00051                                         real *);
00052 
00053 
00054 /*  -- LAPACK auxiliary routine (version 3.0) -- */
00055 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
00056 /*     Courant Institute, Argonne National Lab, and Rice University */
00057 /*     October 31, 1992 */
00058 
00059 /*     .. Scalar Arguments .. */
00060 /*     .. */
00061 
00062 /*  Purpose */
00063 /*  ======= */
00064 
00065 /*  SLAMCH determines single precision machine parameters. */
00066 
00067 /*  Arguments */
00068 /*  ========= */
00069 
00070 /*  CMACH   (input) CHARACTER*1 */
00071 /*          Specifies the value to be returned by SLAMCH: */
00072 /*          = 'E' or 'e',   SLAMCH := eps */
00073 /*          = 'S' or 's ,   SLAMCH := sfmin */
00074 /*          = 'B' or 'b',   SLAMCH := base */
00075 /*          = 'P' or 'p',   SLAMCH := eps*base */
00076 /*          = 'N' or 'n',   SLAMCH := t */
00077 /*          = 'R' or 'r',   SLAMCH := rnd */
00078 /*          = 'M' or 'm',   SLAMCH := emin */
00079 /*          = 'U' or 'u',   SLAMCH := rmin */
00080 /*          = 'L' or 'l',   SLAMCH := emax */
00081 /*          = 'O' or 'o',   SLAMCH := rmax */
00082 
00083 /*          where */
00084 
00085 /*          eps   = relative machine precision */
00086 /*          sfmin = safe minimum, such that 1/sfmin does not overflow */
00087 /*          base  = base of the machine */
00088 /*          prec  = eps*base */
00089 /*          t     = number of (base) digits in the mantissa */
00090 /*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise */
00091 /*          emin  = minimum exponent before (gradual) underflow */
00092 /*          rmin  = underflow threshold - base**(emin-1) */
00093 /*          emax  = largest exponent before overflow */
00094 /*          rmax  = overflow threshold  - (base**emax)*(1-eps) */
00095 
00096 /* ===================================================================== */
00097 
00098 /*     .. Parameters .. */
00099 /*     .. */
00100 /*     .. Local Scalars .. */
00101 /*     .. */
00102 /*     .. External Functions .. */
00103 /*     .. */
00104 /*     .. External Subroutines .. */
00105 /*     .. */
00106 /*     .. Save statement .. */
00107 /*     .. */
00108 /*     .. Data statements .. */
00109 /*     .. */
00110 /*     .. Executable Statements .. */
00111 
00112     if (first) {
00113         first = FALSE_;
00114         slamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
00115         base = (real) beta;
00116         t = (real) it;
00117         if (lrnd) {
00118             rnd = 1.f;
00119             i__1 = 1 - it;
00120             eps = pow_ri(&base, &i__1) / 2;
00121         }
00122         else {
00123             rnd = 0.f;
00124             i__1 = 1 - it;
00125             eps = pow_ri(&base, &i__1);
00126         }
00127         prec = eps * base;
00128         emin = (real) imin;
00129         emax = (real) imax;
00130         sfmin = rmin;
00131         small = 1.f / rmax;
00132         if (small >= sfmin) {
00133 
00134 /*           Use SMALL plus a bit, to avoid the possibility of rounding */
00135 /*           causing overflow when computing  1/sfmin. */
00136 
00137             sfmin = small * (eps + 1.f);
00138         }
00139     }
00140 
00141     if (lsame_(cmach, "E", (ftnlen) 1, (ftnlen) 1)) {
00142         rmach = eps;
00143     }
00144     else if (lsame_(cmach, "S", (ftnlen) 1, (ftnlen) 1)) {
00145         rmach = sfmin;
00146     }
00147     else if (lsame_(cmach, "B", (ftnlen) 1, (ftnlen) 1)) {
00148         rmach = base;
00149     }
00150     else if (lsame_(cmach, "P", (ftnlen) 1, (ftnlen) 1)) {
00151         rmach = prec;
00152     }
00153     else if (lsame_(cmach, "N", (ftnlen) 1, (ftnlen) 1)) {
00154         rmach = t;
00155     }
00156     else if (lsame_(cmach, "R", (ftnlen) 1, (ftnlen) 1)) {
00157         rmach = rnd;
00158     }
00159     else if (lsame_(cmach, "M", (ftnlen) 1, (ftnlen) 1)) {
00160         rmach = emin;
00161     }
00162     else if (lsame_(cmach, "U", (ftnlen) 1, (ftnlen) 1)) {
00163         rmach = rmin;
00164     }
00165     else if (lsame_(cmach, "L", (ftnlen) 1, (ftnlen) 1)) {
00166         rmach = emax;
00167     }
00168     else if (lsame_(cmach, "O", (ftnlen) 1, (ftnlen) 1)) {
00169         rmach = rmax;
00170     }
00171 
00172     ret_val = rmach;
00173     return ret_val;
00174 
00175 /*     End of SLAMCH */
00176 
00177 }                               /* slamch_ */
00178 
00179 
00180 /* *********************************************************************** */
00181 
00182 /* Subroutine */ int
00183 slamc1_(integer * beta, integer * t, logical * rnd, logical * ieee1)
00184 {
00185     /* Initialized data */
00186 
00187     static logical first = TRUE_;
00188 
00189     /* System generated locals */
00190     real r__1, r__2;
00191 
00192     /* Local variables */
00193     static real a, b, c__, f, t1, t2;
00194     static integer lt;
00195     static real one, qtr;
00196     static logical lrnd;
00197     static integer lbeta;
00198     static real savec;
00199     static logical lieee1;
00200     extern doublereal slamc3_(real *, real *);
00201 
00202 
00203 /*  -- LAPACK auxiliary routine (version 3.0) -- */
00204 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
00205 /*     Courant Institute, Argonne National Lab, and Rice University */
00206 /*     October 31, 1992 */
00207 
00208 /*     .. Scalar Arguments .. */
00209 /*     .. */
00210 
00211 /*  Purpose */
00212 /*  ======= */
00213 
00214 /*  SLAMC1 determines the machine parameters given by BETA, T, RND, and */
00215 /*  IEEE1. */
00216 
00217 /*  Arguments */
00218 /*  ========= */
00219 
00220 /*  BETA    (output) INTEGER */
00221 /*          The base of the machine. */
00222 
00223 /*  T       (output) INTEGER */
00224 /*          The number of ( BETA ) digits in the mantissa. */
00225 
00226 /*  RND     (output) LOGICAL */
00227 /*          Specifies whether proper rounding  ( RND = .TRUE. )  or */
00228 /*          chopping  ( RND = .FALSE. )  occurs in addition. This may not */
00229 /*          be a reliable guide to the way in which the machine performs */
00230 /*          its arithmetic. */
00231 
00232 /*  IEEE1   (output) LOGICAL */
00233 /*          Specifies whether rounding appears to be done in the IEEE */
00234 /*          'round to nearest' style. */
00235 
00236 /*  Further Details */
00237 /*  =============== */
00238 
00239 /*  The routine is based on the routine  ENVRON  by Malcolm and */
00240 /*  incorporates suggestions by Gentleman and Marovich. See */
00241 
00242 /*     Malcolm M. A. (1972) Algorithms to reveal properties of */
00243 /*        floating-point arithmetic. Comms. of the ACM, 15, 949-951. */
00244 
00245 /*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms */
00246 /*        that reveal properties of floating point arithmetic units. */
00247 /*        Comms. of the ACM, 17, 276-277. */
00248 
00249 /* ===================================================================== */
00250 
00251 /*     .. Local Scalars .. */
00252 /*     .. */
00253 /*     .. External Functions .. */
00254 /*     .. */
00255 /*     .. Save statement .. */
00256 /*     .. */
00257 /*     .. Data statements .. */
00258 /*     .. */
00259 /*     .. Executable Statements .. */
00260 
00261     if (first) {
00262         first = FALSE_;
00263         one = 1.f;
00264 
00265 /*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA, */
00266 /*        IEEE1, T and RND. */
00267 
00268 /*        Throughout this routine  we use the function  SLAMC3  to ensure */
00269 /*        that relevant values are  stored and not held in registers,  or */
00270 /*        are not affected by optimizers. */
00271 
00272 /*        Compute  a = 2.0**m  with the  smallest positive integer m such */
00273 /*        that */
00274 
00275 /*           fl( a + 1.0 ) = a. */
00276 
00277         a = 1.f;
00278         c__ = 1.f;
00279 
00280 /* +       WHILE( C.EQ.ONE )LOOP */
00281       L10:
00282         if (c__ == one) {
00283             a *= 2;
00284             c__ = slamc3_(&a, &one);
00285             r__1 = -a;
00286             c__ = slamc3_(&c__, &r__1);
00287             goto L10;
00288         }
00289 /* +       END WHILE */
00290 
00291 /*        Now compute  b = 2.0**m  with the smallest positive integer m */
00292 /*        such that */
00293 
00294 /*           fl( a + b ) .gt. a. */
00295 
00296         b = 1.f;
00297         c__ = slamc3_(&a, &b);
00298 
00299 /* +       WHILE( C.EQ.A )LOOP */
00300       L20:
00301         if (c__ == a) {
00302             b *= 2;
00303             c__ = slamc3_(&a, &b);
00304             goto L20;
00305         }
00306 /* +       END WHILE */
00307 
00308 /*        Now compute the base.  a and c  are neighbouring floating point */
00309 /*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so */
00310 /*        their difference is beta. Adding 0.25 to c is to ensure that it */
00311 /*        is truncated to beta and not ( beta - 1 ). */
00312 
00313         qtr = one / 4;
00314         savec = c__;
00315         r__1 = -a;
00316         c__ = slamc3_(&c__, &r__1);
00317         lbeta = c__ + qtr;
00318 
00319 /*        Now determine whether rounding or chopping occurs,  by adding a */
00320 /*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a. */
00321 
00322         b = (real) lbeta;
00323         r__1 = b / 2;
00324         r__2 = -b / 100;
00325         f = slamc3_(&r__1, &r__2);
00326         c__ = slamc3_(&f, &a);
00327         if (c__ == a) {
00328             lrnd = TRUE_;
00329         }
00330         else {
00331             lrnd = FALSE_;
00332         }
00333         r__1 = b / 2;
00334         r__2 = b / 100;
00335         f = slamc3_(&r__1, &r__2);
00336         c__ = slamc3_(&f, &a);
00337         if (lrnd && c__ == a) {
00338             lrnd = FALSE_;
00339         }
00340 
00341 /*        Try and decide whether rounding is done in the  IEEE  'round to */
00342 /*        nearest' style. B/2 is half a unit in the last place of the two */
00343 /*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit */
00344 /*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change */
00345 /*        A, but adding B/2 to SAVEC should change SAVEC. */
00346 
00347         r__1 = b / 2;
00348         t1 = slamc3_(&r__1, &a);
00349         r__1 = b / 2;
00350         t2 = slamc3_(&r__1, &savec);
00351         lieee1 = t1 == a && t2 > savec && lrnd;
00352 
00353 /*        Now find  the  mantissa, t.  It should  be the  integer part of */
00354 /*        log to the base beta of a,  however it is safer to determine  t */
00355 /*        by powering.  So we find t as the smallest positive integer for */
00356 /*        which */
00357 
00358 /*           fl( beta**t + 1.0 ) = 1.0. */
00359 
00360         lt = 0;
00361         a = 1.f;
00362         c__ = 1.f;
00363 
00364 /* +       WHILE( C.EQ.ONE )LOOP */
00365       L30:
00366         if (c__ == one) {
00367             ++lt;
00368             a *= lbeta;
00369             c__ = slamc3_(&a, &one);
00370             r__1 = -a;
00371             c__ = slamc3_(&c__, &r__1);
00372             goto L30;
00373         }
00374 /* +       END WHILE */
00375 
00376     }
00377 
00378     *beta = lbeta;
00379     *t = lt;
00380     *rnd = lrnd;
00381     *ieee1 = lieee1;
00382     return 0;
00383 
00384 /*     End of SLAMC1 */
00385 
00386 }                               /* slamc1_ */
00387 
00388 
00389 /* *********************************************************************** */
00390 
00391 /* Subroutine */ int
00392 slamc2_(integer * beta, integer * t, logical * rnd, real *
00393         eps, integer * emin, real * rmin, integer * emax, real * rmax)
00394 {
00395     /* Initialized data */
00396 
00397     static logical first = TRUE_;
00398     static logical iwarn = FALSE_;
00399 
00400     /* Format strings */
00401     static char fmt_9999[] =
00402         "(//\002 WARNING. The value EMIN may be incorre"
00403         "ct:-\002,\002  EMIN = \002,i8,/\002 If, after inspection, the va"
00404         "lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
00405         " the IF block as marked within the code of routine\002,\002 SLAM"
00406         "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
00407 
00408     /* System generated locals */
00409     integer i__1;
00410     real r__1, r__2, r__3, r__4, r__5;
00411 
00412     /* Builtin functions */
00413     double pow_ri(real *, integer *);
00414     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen),
00415         e_wsfe(void);
00416 
00417     /* Local variables */
00418     static real a, b, c__;
00419     static integer i__, lt;
00420     static real one, two;
00421     static logical ieee;
00422     static real half;
00423     static logical lrnd;
00424     static real leps, zero;
00425     static integer lbeta;
00426     static real rbase;
00427     static integer lemin, lemax, gnmin;
00428     static real small;
00429     static integer gpmin;
00430     static real third, lrmin, lrmax, sixth;
00431     static logical lieee1;
00432     extern /* Subroutine */ int slamc1_(integer *, integer *, logical *,
00433                                         logical *);
00434     extern doublereal slamc3_(real *, real *);
00435     extern /* Subroutine */ int slamc4_(integer *, real *, integer *),
00436         slamc5_(integer *, integer *, integer *, logical *, integer *,
00437                 real *);
00438     static integer ngnmin, ngpmin;
00439 
00440     /* Fortran I/O blocks */
00441     static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
00442 
00443 
00444 
00445 /*  -- LAPACK auxiliary routine (version 3.0) -- */
00446 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
00447 /*     Courant Institute, Argonne National Lab, and Rice University */
00448 /*     October 31, 1992 */
00449 
00450 /*     .. Scalar Arguments .. */
00451 /*     .. */
00452 
00453 /*  Purpose */
00454 /*  ======= */
00455 
00456 /*  SLAMC2 determines the machine parameters specified in its argument */
00457 /*  list. */
00458 
00459 /*  Arguments */
00460 /*  ========= */
00461 
00462 /*  BETA    (output) INTEGER */
00463 /*          The base of the machine. */
00464 
00465 /*  T       (output) INTEGER */
00466 /*          The number of ( BETA ) digits in the mantissa. */
00467 
00468 /*  RND     (output) LOGICAL */
00469 /*          Specifies whether proper rounding  ( RND = .TRUE. )  or */
00470 /*          chopping  ( RND = .FALSE. )  occurs in addition. This may not */
00471 /*          be a reliable guide to the way in which the machine performs */
00472 /*          its arithmetic. */
00473 
00474 /*  EPS     (output) REAL */
00475 /*          The smallest positive number such that */
00476 
00477 /*             fl( 1.0 - EPS ) .LT. 1.0, */
00478 
00479 /*          where fl denotes the computed value. */
00480 
00481 /*  EMIN    (output) INTEGER */
00482 /*          The minimum exponent before (gradual) underflow occurs. */
00483 
00484 /*  RMIN    (output) REAL */
00485 /*          The smallest normalized number for the machine, given by */
00486 /*          BASE**( EMIN - 1 ), where  BASE  is the floating point value */
00487 /*          of BETA. */
00488 
00489 /*  EMAX    (output) INTEGER */
00490 /*          The maximum exponent before overflow occurs. */
00491 
00492 /*  RMAX    (output) REAL */
00493 /*          The largest positive number for the machine, given by */
00494 /*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point */
00495 /*          value of BETA. */
00496 
00497 /*  Further Details */
00498 /*  =============== */
00499 
00500 /*  The computation of  EPS  is based on a routine PARANOIA by */
00501 /*  W. Kahan of the University of California at Berkeley. */
00502 
00503 /* ===================================================================== */
00504 
00505 /*     .. Local Scalars .. */
00506 /*     .. */
00507 /*     .. External Functions .. */
00508 /*     .. */
00509 /*     .. External Subroutines .. */
00510 /*     .. */
00511 /*     .. Intrinsic Functions .. */
00512 /*     .. */
00513 /*     .. Save statement .. */
00514 /*     .. */
00515 /*     .. Data statements .. */
00516 /*     .. */
00517 /*     .. Executable Statements .. */
00518 
00519     if (first) {
00520         first = FALSE_;
00521         zero = 0.f;
00522         one = 1.f;
00523         two = 2.f;
00524 
00525 /*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of */
00526 /*        BETA, T, RND, EPS, EMIN and RMIN. */
00527 
00528 /*        Throughout this routine  we use the function  SLAMC3  to ensure */
00529 /*        that relevant values are stored  and not held in registers,  or */
00530 /*        are not affected by optimizers. */
00531 
00532 /*        SLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1. */
00533 
00534         slamc1_(&lbeta, &lt, &lrnd, &lieee1);
00535 
00536 /*        Start to find EPS. */
00537 
00538         b = (real) lbeta;
00539         i__1 = -lt;
00540         a = pow_ri(&b, &i__1);
00541         leps = a;
00542 
00543 /*        Try some tricks to see whether or not this is the correct  EPS. */
00544 
00545         b = two / 3;
00546         half = one / 2;
00547         r__1 = -half;
00548         sixth = slamc3_(&b, &r__1);
00549         third = slamc3_(&sixth, &sixth);
00550         r__1 = -half;
00551         b = slamc3_(&third, &r__1);
00552         b = slamc3_(&b, &sixth);
00553         b = dabs(b);
00554         if (b < leps) {
00555             b = leps;
00556         }
00557 
00558         leps = 1.f;
00559 
00560 /* +       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
00561       L10:
00562         if (leps > b && b > zero) {
00563             leps = b;
00564             r__1 = half * leps;
00565 /* Computing 5th power */
00566             r__3 = two, r__4 = r__3, r__3 *= r__3;
00567 /* Computing 2nd power */
00568             r__5 = leps;
00569             r__2 = r__4 * (r__3 * r__3) * (r__5 * r__5);
00570             c__ = slamc3_(&r__1, &r__2);
00571             r__1 = -c__;
00572             c__ = slamc3_(&half, &r__1);
00573             b = slamc3_(&half, &c__);
00574             r__1 = -b;
00575             c__ = slamc3_(&half, &r__1);
00576             b = slamc3_(&half, &c__);
00577             goto L10;
00578         }
00579 /* +       END WHILE */
00580 
00581         if (a < leps) {
00582             leps = a;
00583         }
00584 
00585 /*        Computation of EPS complete. */
00586 
00587 /*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)). */
00588 /*        Keep dividing  A by BETA until (gradual) underflow occurs. This */
00589 /*        is detected when we cannot recover the previous A. */
00590 
00591         rbase = one / lbeta;
00592         small = one;
00593         for (i__ = 1; i__ <= 3; ++i__) {
00594             r__1 = small * rbase;
00595             small = slamc3_(&r__1, &zero);
00596 /* L20: */
00597         }
00598         a = slamc3_(&one, &small);
00599         slamc4_(&ngpmin, &one, &lbeta);
00600         r__1 = -one;
00601         slamc4_(&ngnmin, &r__1, &lbeta);
00602         slamc4_(&gpmin, &a, &lbeta);
00603         r__1 = -a;
00604         slamc4_(&gnmin, &r__1, &lbeta);
00605         ieee = FALSE_;
00606 
00607         if (ngpmin == ngnmin && gpmin == gnmin) {
00608             if (ngpmin == gpmin) {
00609                 lemin = ngpmin;
00610 /*            ( Non twos-complement machines, no gradual underflow; */
00611 /*              e.g.,  VAX ) */
00612             }
00613             else if (gpmin - ngpmin == 3) {
00614                 lemin = ngpmin - 1 + lt;
00615                 ieee = TRUE_;
00616 /*            ( Non twos-complement machines, with gradual underflow; */
00617 /*              e.g., IEEE standard followers ) */
00618             }
00619             else {
00620                 lemin = min(ngpmin, gpmin);
00621 /*            ( A guess; no known machine ) */
00622                 iwarn = TRUE_;
00623             }
00624 
00625         }
00626         else if (ngpmin == gpmin && ngnmin == gnmin) {
00627             if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
00628                 lemin = max(ngpmin, ngnmin);
00629 /*            ( Twos-complement machines, no gradual underflow; */
00630 /*              e.g., CYBER 205 ) */
00631             }
00632             else {
00633                 lemin = min(ngpmin, ngnmin);
00634 /*            ( A guess; no known machine ) */
00635                 iwarn = TRUE_;
00636             }
00637 
00638         }
00639         else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1
00640                  && gpmin == gnmin) {
00641             if (gpmin - min(ngpmin, ngnmin) == 3) {
00642                 lemin = max(ngpmin, ngnmin) - 1 + lt;
00643 /*            ( Twos-complement machines with gradual underflow; */
00644 /*              no known machine ) */
00645             }
00646             else {
00647                 lemin = min(ngpmin, ngnmin);
00648 /*            ( A guess; no known machine ) */
00649                 iwarn = TRUE_;
00650             }
00651 
00652         }
00653         else {
00654 /* Computing MIN */
00655             i__1 = min(ngpmin, ngnmin), i__1 = min(i__1, gpmin);
00656             lemin = min(i__1, gnmin);
00657 /*         ( A guess; no known machine ) */
00658             iwarn = TRUE_;
00659         }
00660 /* ** */
00661 /* Comment out this if block if EMIN is ok */
00662         if (iwarn) {
00663             first = TRUE_;
00664             s_wsfe(&io___58);
00665             do_fio(&c__1, (char *) &lemin, (ftnlen) sizeof(integer));
00666             e_wsfe();
00667         }
00668 /* ** */
00669 
00670 /*        Assume IEEE arithmetic if we found denormalised  numbers above, */
00671 /*        or if arithmetic seems to round in the  IEEE style,  determined */
00672 /*        in routine SLAMC1. A true IEEE machine should have both  things */
00673 /*        true; however, faulty machines may have one or the other. */
00674 
00675         ieee = ieee || lieee1;
00676 
00677 /*        Compute  RMIN by successive division by  BETA. We could compute */
00678 /*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during */
00679 /*        this computation. */
00680 
00681         lrmin = 1.f;
00682         i__1 = 1 - lemin;
00683         for (i__ = 1; i__ <= i__1; ++i__) {
00684             r__1 = lrmin * rbase;
00685             lrmin = slamc3_(&r__1, &zero);
00686 /* L30: */
00687         }
00688 
00689 /*        Finally, call SLAMC5 to compute EMAX and RMAX. */
00690 
00691         slamc5_(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);
00692     }
00693 
00694     *beta = lbeta;
00695     *t = lt;
00696     *rnd = lrnd;
00697     *eps = leps;
00698     *emin = lemin;
00699     *rmin = lrmin;
00700     *emax = lemax;
00701     *rmax = lrmax;
00702 
00703     return 0;
00704 
00705 
00706 /*     End of SLAMC2 */
00707 
00708 }                               /* slamc2_ */
00709 
00710 
00711 /* *********************************************************************** */
00712 
00713 doublereal
00714 slamc3_(real * a, real * b)
00715 {
00716     /* System generated locals */
00717     real ret_val;
00718 
00719 
00720 /*  -- LAPACK auxiliary routine (version 3.0) -- */
00721 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
00722 /*     Courant Institute, Argonne National Lab, and Rice University */
00723 /*     October 31, 1992 */
00724 
00725 /*     .. Scalar Arguments .. */
00726 /*     .. */
00727 
00728 /*  Purpose */
00729 /*  ======= */
00730 
00731 /*  SLAMC3  is intended to force  A  and  B  to be stored prior to doing */
00732 /*  the addition of  A  and  B ,  for use in situations where optimizers */
00733 /*  might hold one of these in a register. */
00734 
00735 /*  Arguments */
00736 /*  ========= */
00737 
00738 /*  A, B    (input) REAL */
00739 /*          The values A and B. */
00740 
00741 /* ===================================================================== */
00742 
00743 /*     .. Executable Statements .. */
00744 
00745     ret_val = *a + *b;
00746 
00747     return ret_val;
00748 
00749 /*     End of SLAMC3 */
00750 
00751 }                               /* slamc3_ */
00752 
00753 
00754 /* *********************************************************************** */
00755 
00756 /* Subroutine */ int
00757 slamc4_(integer * emin, real * start, integer * base)
00758 {
00759     /* System generated locals */
00760     integer i__1;
00761     real r__1;
00762 
00763     /* Local variables */
00764     static real a;
00765     static integer i__;
00766     static real b1, b2, c1, c2, d1, d2, one, zero, rbase;
00767     extern doublereal slamc3_(real *, real *);
00768 
00769 
00770 /*  -- LAPACK auxiliary routine (version 3.0) -- */
00771 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
00772 /*     Courant Institute, Argonne National Lab, and Rice University */
00773 /*     October 31, 1992 */
00774 
00775 /*     .. Scalar Arguments .. */
00776 /*     .. */
00777 
00778 /*  Purpose */
00779 /*  ======= */
00780 
00781 /*  SLAMC4 is a service routine for SLAMC2. */
00782 
00783 /*  Arguments */
00784 /*  ========= */
00785 
00786 /*  EMIN    (output) EMIN */
00787 /*          The minimum exponent before (gradual) underflow, computed by */
00788 /*          setting A = START and dividing by BASE until the previous A */
00789 /*          can not be recovered. */
00790 
00791 /*  START   (input) REAL */
00792 /*          The starting point for determining EMIN. */
00793 
00794 /*  BASE    (input) INTEGER */
00795 /*          The base of the machine. */
00796 
00797 /* ===================================================================== */
00798 
00799 /*     .. Local Scalars .. */
00800 /*     .. */
00801 /*     .. External Functions .. */
00802 /*     .. */
00803 /*     .. Executable Statements .. */
00804 
00805     a = *start;
00806     one = 1.f;
00807     rbase = one / *base;
00808     zero = 0.f;
00809     *emin = 1;
00810     r__1 = a * rbase;
00811     b1 = slamc3_(&r__1, &zero);
00812     c1 = a;
00813     c2 = a;
00814     d1 = a;
00815     d2 = a;
00816 /* +    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
00817 /*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP */
00818   L10:
00819     if (c1 == a && c2 == a && d1 == a && d2 == a) {
00820         --(*emin);
00821         a = b1;
00822         r__1 = a / *base;
00823         b1 = slamc3_(&r__1, &zero);
00824         r__1 = b1 * *base;
00825         c1 = slamc3_(&r__1, &zero);
00826         d1 = zero;
00827         i__1 = *base;
00828         for (i__ = 1; i__ <= i__1; ++i__) {
00829             d1 += b1;
00830 /* L20: */
00831         }
00832         r__1 = a * rbase;
00833         b2 = slamc3_(&r__1, &zero);
00834         r__1 = b2 / rbase;
00835         c2 = slamc3_(&r__1, &zero);
00836         d2 = zero;
00837         i__1 = *base;
00838         for (i__ = 1; i__ <= i__1; ++i__) {
00839             d2 += b2;
00840 /* L30: */
00841         }
00842         goto L10;
00843     }
00844 /* +    END WHILE */
00845 
00846     return 0;
00847 
00848 /*     End of SLAMC4 */
00849 
00850 }                               /* slamc4_ */
00851 
00852 
00853 /* *********************************************************************** */
00854 
00855 /* Subroutine */ int
00856 slamc5_(integer * beta, integer * p, integer * emin,
00857         logical * ieee, integer * emax, real * rmax)
00858 {
00859     /* System generated locals */
00860     integer i__1;
00861     real r__1;
00862 
00863     /* Local variables */
00864     static integer i__;
00865     static real y, z__;
00866     static integer try__, lexp;
00867     static real oldy;
00868     static integer uexp, nbits;
00869     extern doublereal slamc3_(real *, real *);
00870     static real recbas;
00871     static integer exbits, expsum;
00872 
00873 
00874 /*  -- LAPACK auxiliary routine (version 3.0) -- */
00875 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
00876 /*     Courant Institute, Argonne National Lab, and Rice University */
00877 /*     October 31, 1992 */
00878 
00879 /*     .. Scalar Arguments .. */
00880 /*     .. */
00881 
00882 /*  Purpose */
00883 /*  ======= */
00884 
00885 /*  SLAMC5 attempts to compute RMAX, the largest machine floating-point */
00886 /*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum */
00887 /*  approximately to a power of 2.  It will fail on machines where this */
00888 /*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
00889 /*  EMAX = 28718).  It will also fail if the value supplied for EMIN is */
00890 /*  too large (i.e. too close to zero), probably with overflow. */
00891 
00892 /*  Arguments */
00893 /*  ========= */
00894 
00895 /*  BETA    (input) INTEGER */
00896 /*          The base of floating-point arithmetic. */
00897 
00898 /*  P       (input) INTEGER */
00899 /*          The number of base BETA digits in the mantissa of a */
00900 /*          floating-point value. */
00901 
00902 /*  EMIN    (input) INTEGER */
00903 /*          The minimum exponent before (gradual) underflow. */
00904 
00905 /*  IEEE    (input) LOGICAL */
00906 /*          A logical flag specifying whether or not the arithmetic */
00907 /*          system is thought to comply with the IEEE standard. */
00908 
00909 /*  EMAX    (output) INTEGER */
00910 /*          The largest exponent before overflow */
00911 
00912 /*  RMAX    (output) REAL */
00913 /*          The largest machine floating-point number. */
00914 
00915 /* ===================================================================== */
00916 
00917 /*     .. Parameters .. */
00918 /*     .. */
00919 /*     .. Local Scalars .. */
00920 /*     .. */
00921 /*     .. External Functions .. */
00922 /*     .. */
00923 /*     .. Intrinsic Functions .. */
00924 /*     .. */
00925 /*     .. Executable Statements .. */
00926 
00927 /*     First compute LEXP and UEXP, two powers of 2 that bound */
00928 /*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
00929 /*     approximately to the bound that is closest to abs(EMIN). */
00930 /*     (EMAX is the exponent of the required number RMAX). */
00931 
00932     lexp = 1;
00933     exbits = 1;
00934   L10:
00935     try__ = lexp << 1;
00936     if (try__ <= -(*emin)) {
00937         lexp = try__;
00938         ++exbits;
00939         goto L10;
00940     }
00941     if (lexp == -(*emin)) {
00942         uexp = lexp;
00943     }
00944     else {
00945         uexp = try__;
00946         ++exbits;
00947     }
00948 
00949 /*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
00950 /*     than or equal to EMIN. EXBITS is the number of bits needed to */
00951 /*     store the exponent. */
00952 
00953     if (uexp + *emin > -lexp - *emin) {
00954         expsum = lexp << 1;
00955     }
00956     else {
00957         expsum = uexp << 1;
00958     }
00959 
00960 /*     EXPSUM is the exponent range, approximately equal to */
00961 /*     EMAX - EMIN + 1 . */
00962 
00963     *emax = expsum + *emin - 1;
00964     nbits = exbits + 1 + *p;
00965 
00966 /*     NBITS is the total number of bits needed to store a */
00967 /*     floating-point number. */
00968 
00969     if (nbits % 2 == 1 && *beta == 2) {
00970 
00971 /*        Either there are an odd number of bits used to store a */
00972 /*        floating-point number, which is unlikely, or some bits are */
00973 /*        not used in the representation of numbers, which is possible, */
00974 /*        (e.g. Cray machines) or the mantissa has an implicit bit, */
00975 /*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
00976 /*        most likely. We have to assume the last alternative. */
00977 /*        If this is true, then we need to reduce EMAX by one because */
00978 /*        there must be some way of representing zero in an implicit-bit */
00979 /*        system. On machines like Cray, we are reducing EMAX by one */
00980 /*        unnecessarily. */
00981 
00982         --(*emax);
00983     }
00984 
00985     if (*ieee) {
00986 
00987 /*        Assume we are on an IEEE machine which reserves one exponent */
00988 /*        for infinity and NaN. */
00989 
00990         --(*emax);
00991     }
00992 
00993 /*     Now create RMAX, the largest machine number, which should */
00994 /*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
00995 
00996 /*     First compute 1.0 - BETA**(-P), being careful that the */
00997 /*     result is less than 1.0 . */
00998 
00999     recbas = 1.f / *beta;
01000     z__ = *beta - 1.f;
01001     y = 0.f;
01002     i__1 = *p;
01003     for (i__ = 1; i__ <= i__1; ++i__) {
01004         z__ *= recbas;
01005         if (y < 1.f) {
01006             oldy = y;
01007         }
01008         y = slamc3_(&y, &z__);
01009 /* L20: */
01010     }
01011     if (y >= 1.f) {
01012         y = oldy;
01013     }
01014 
01015 /*     Now multiply by BETA**EMAX to get RMAX. */
01016 
01017     i__1 = *emax;
01018     for (i__ = 1; i__ <= i__1; ++i__) {
01019         r__1 = y * *beta;
01020         y = slamc3_(&r__1, &c_b32);
01021 /* L30: */
01022     }
01023 
01024     *rmax = y;
01025     return 0;
01026 
01027 /*     End of SLAMC5 */
01028 
01029 }                               /* slamc5_ */

Generated on Fri Jan 14 2011 for SphinxBase by  doxygen 1.7.1