?? slamc1.c
字號:
/** ======================================================================* NIST Guide to Available Math Software.* Fullsource for module SSYEVX.C from package CLAPACK.* Retrieved from NETLIB on Fri Mar 10 14:23:44 2000.* ======================================================================*/#include <f2c.h>/* Subroutine */ int slamc1_(integer *beta, integer *t, logical *rnd, logical *ieee1){/* -- LAPACK auxiliary routine (version 2.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 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. ===================================================================== */ /* Initialized data */ static logical first = TRUE_; /* System generated locals */ real r__1, r__2; /* Local variables */ static logical lrnd; static real a, b, c, f; static integer lbeta; static real savec; static logical lieee1; static real t1, t2; extern doublereal slamc3_(real *, real *); static integer lt; static real one, qtr; 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 = -(doublereal)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 = -(doublereal)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 = -(doublereal)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 = -(doublereal)a; c = slamc3_(&c, &r__1); goto L30; }/* + END WHILE */ } *beta = lbeta; *t = lt; *rnd = lrnd; *ieee1 = lieee1; return 0;/* End of SLAMC1 */} /* slamc1_ */
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -