SphinxBase
0.6
|
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 "sphinxbase/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, <, &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, <, &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_ */