|
kusano |
7d535a |
#include <stdio.h></stdio.h>
|
|
kusano |
7d535a |
#define TRUE_ (1)
|
|
kusano |
7d535a |
#define FALSE_ (0)
|
|
kusano |
7d535a |
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
|
kusano |
7d535a |
#define min(a,b) ((a) <= (b) ? (a) : (b))
|
|
kusano |
7d535a |
#define max(a,b) ((a) >= (b) ? (a) : (b))
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
double dlamch_(char *cmach)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
/* -- LAPACK auxiliary routine (version 2.0) --
|
|
kusano |
7d535a |
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
kusano |
7d535a |
Courant Institute, Argonne National Lab, and Rice University
|
|
kusano |
7d535a |
October 31, 1992
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Purpose
|
|
kusano |
7d535a |
=======
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DLAMCH determines double precision machine parameters.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Arguments
|
|
kusano |
7d535a |
=========
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
CMACH (input) CHARACTER*1
|
|
kusano |
7d535a |
Specifies the value to be returned by DLAMCH:
|
|
kusano |
7d535a |
= 'E' or 'e', DLAMCH := eps
|
|
kusano |
7d535a |
= 'S' or 's , DLAMCH := sfmin
|
|
kusano |
7d535a |
= 'B' or 'b', DLAMCH := base
|
|
kusano |
7d535a |
= 'P' or 'p', DLAMCH := eps*base
|
|
kusano |
7d535a |
= 'N' or 'n', DLAMCH := t
|
|
kusano |
7d535a |
= 'R' or 'r', DLAMCH := rnd
|
|
kusano |
7d535a |
= 'M' or 'm', DLAMCH := emin
|
|
kusano |
7d535a |
= 'U' or 'u', DLAMCH := rmin
|
|
kusano |
7d535a |
= 'L' or 'l', DLAMCH := emax
|
|
kusano |
7d535a |
= 'O' or 'o', DLAMCH := rmax
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
where
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
eps = relative machine precision
|
|
kusano |
7d535a |
sfmin = safe minimum, such that 1/sfmin does not overflow
|
|
kusano |
7d535a |
base = base of the machine
|
|
kusano |
7d535a |
prec = eps*base
|
|
kusano |
7d535a |
t = number of (base) digits in the mantissa
|
|
kusano |
7d535a |
rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
|
|
kusano |
7d535a |
emin = minimum exponent before (gradual) underflow
|
|
kusano |
7d535a |
rmin = underflow threshold - base**(emin-1)
|
|
kusano |
7d535a |
emax = largest exponent before overflow
|
|
kusano |
7d535a |
rmax = overflow threshold - (base**emax)*(1-eps)
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
=====================================================================
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
static int first = TRUE_;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* System generated locals */
|
|
kusano |
7d535a |
int i__1;
|
|
kusano |
7d535a |
double ret_val;
|
|
kusano |
7d535a |
/* Builtin functions */
|
|
kusano |
7d535a |
double pow_di(double *, int *);
|
|
kusano |
7d535a |
/* Local variables */
|
|
kusano |
7d535a |
static double base;
|
|
kusano |
7d535a |
static int beta;
|
|
kusano |
7d535a |
static double emin, prec, emax;
|
|
kusano |
7d535a |
static int imin, imax;
|
|
kusano |
7d535a |
static int lrnd;
|
|
kusano |
7d535a |
static double rmin, rmax, t, rmach;
|
|
kusano |
7d535a |
/* extern int lsame_(char *, char *);*/
|
|
kusano |
7d535a |
static double small, sfmin;
|
|
kusano |
7d535a |
extern /* Subroutine */ int dlamc2_(int *, int *, int *,
|
|
kusano |
7d535a |
double *, int *, double *, int *, double *);
|
|
kusano |
7d535a |
static int it;
|
|
kusano |
7d535a |
static double rnd, eps;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (first) {
|
|
kusano |
7d535a |
first = FALSE_;
|
|
kusano |
7d535a |
dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
|
|
kusano |
7d535a |
base = (double) beta;
|
|
kusano |
7d535a |
t = (double) it;
|
|
kusano |
7d535a |
if (lrnd) {
|
|
kusano |
7d535a |
rnd = 1.;
|
|
kusano |
7d535a |
i__1 = 1 - it;
|
|
kusano |
7d535a |
eps = pow_di(&base, &i__1) / 2;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
rnd = 0.;
|
|
kusano |
7d535a |
i__1 = 1 - it;
|
|
kusano |
7d535a |
eps = pow_di(&base, &i__1);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
prec = eps * base;
|
|
kusano |
7d535a |
emin = (double) imin;
|
|
kusano |
7d535a |
emax = (double) imax;
|
|
kusano |
7d535a |
sfmin = rmin;
|
|
kusano |
7d535a |
small = 1. / rmax;
|
|
kusano |
7d535a |
if (small >= sfmin) {
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Use SMALL plus a bit, to avoid the possibility of rounding
|
|
kusano |
7d535a |
causing overflow when computing 1/sfmin. */
|
|
kusano |
7d535a |
sfmin = small * (eps + 1.);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (lsame_(cmach, "E")) {
|
|
kusano |
7d535a |
rmach = eps;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "S")) {
|
|
kusano |
7d535a |
rmach = sfmin;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "B")) {
|
|
kusano |
7d535a |
rmach = base;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "P")) {
|
|
kusano |
7d535a |
rmach = prec;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "N")) {
|
|
kusano |
7d535a |
rmach = t;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "R")) {
|
|
kusano |
7d535a |
rmach = rnd;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "M")) {
|
|
kusano |
7d535a |
rmach = emin;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "U")) {
|
|
kusano |
7d535a |
rmach = rmin;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "L")) {
|
|
kusano |
7d535a |
rmach = emax;
|
|
kusano |
7d535a |
} else if (lsame_(cmach, "O")) {
|
|
kusano |
7d535a |
rmach = rmax;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
ret_val = rmach;
|
|
kusano |
7d535a |
return ret_val;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* End of DLAMCH */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} /* dlamch_ */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Subroutine */ int dlamc1_(int *beta, int *t, int *rnd, int
|
|
kusano |
7d535a |
*ieee1)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
/* -- LAPACK auxiliary routine (version 2.0) --
|
|
kusano |
7d535a |
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
kusano |
7d535a |
Courant Institute, Argonne National Lab, and Rice University
|
|
kusano |
7d535a |
October 31, 1992
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Purpose
|
|
kusano |
7d535a |
=======
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DLAMC1 determines the machine parameters given by BETA, T, RND, and
|
|
kusano |
7d535a |
IEEE1.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Arguments
|
|
kusano |
7d535a |
=========
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
BETA (output) INT
|
|
kusano |
7d535a |
The base of the machine.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
T (output) INT
|
|
kusano |
7d535a |
The number of ( BETA ) digits in the mantissa.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
RND (output) INT
|
|
kusano |
7d535a |
Specifies whether proper rounding ( RND = .TRUE. ) or
|
|
kusano |
7d535a |
chopping ( RND = .FALSE. ) occurs in addition. This may not
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
be a reliable guide to the way in which the machine performs
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
its arithmetic.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
IEEE1 (output) INT
|
|
kusano |
7d535a |
Specifies whether rounding appears to be done in the IEEE
|
|
kusano |
7d535a |
'round to nearest' style.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Further Details
|
|
kusano |
7d535a |
===============
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
The routine is based on the routine ENVRON by Malcolm and
|
|
kusano |
7d535a |
incorporates suggestions by Gentleman and Marovich. See
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Malcolm M. A. (1972) Algorithms to reveal properties of
|
|
kusano |
7d535a |
floating-point arithmetic. Comms. of the ACM, 15, 949-951.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Gentleman W. M. and Marovich S. B. (1974) More on algorithms
|
|
kusano |
7d535a |
that reveal properties of floating point arithmetic units.
|
|
kusano |
7d535a |
Comms. of the ACM, 17, 276-277.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
=====================================================================
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
/* Initialized data */
|
|
kusano |
7d535a |
static int first = TRUE_;
|
|
kusano |
7d535a |
/* System generated locals */
|
|
kusano |
7d535a |
double d__1, d__2;
|
|
kusano |
7d535a |
/* Local variables */
|
|
kusano |
7d535a |
static int lrnd;
|
|
kusano |
7d535a |
static double a, b, c, f;
|
|
kusano |
7d535a |
static int lbeta;
|
|
kusano |
7d535a |
static double savec;
|
|
kusano |
7d535a |
extern double dlamc3_(double *, double *);
|
|
kusano |
7d535a |
static int lieee1;
|
|
kusano |
7d535a |
static double t1, t2;
|
|
kusano |
7d535a |
static int lt;
|
|
kusano |
7d535a |
static double one, qtr;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (first) {
|
|
kusano |
7d535a |
first = FALSE_;
|
|
kusano |
7d535a |
one = 1.;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* LBETA, LIEEE1, LT and LRND are the local values of BE
|
|
kusano |
7d535a |
TA,
|
|
kusano |
7d535a |
IEEE1, T and RND.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Throughout this routine we use the function DLAMC3 to ens
|
|
kusano |
7d535a |
ure
|
|
kusano |
7d535a |
that relevant values are stored and not held in registers,
|
|
kusano |
7d535a |
or
|
|
kusano |
7d535a |
are not affected by optimizers.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Compute a = 2.0**m with the smallest positive integer m s
|
|
kusano |
7d535a |
uch
|
|
kusano |
7d535a |
that
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
fl( a + 1.0 ) = a. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
a = 1.;
|
|
kusano |
7d535a |
c = 1.;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* + WHILE( C.EQ.ONE )LOOP */
|
|
kusano |
7d535a |
L10:
|
|
kusano |
7d535a |
if (c == one) {
|
|
kusano |
7d535a |
a *= 2;
|
|
kusano |
7d535a |
c = dlamc3_(&a, &one);
|
|
kusano |
7d535a |
d__1 = -a;
|
|
kusano |
7d535a |
c = dlamc3_(&c, &d__1);
|
|
kusano |
7d535a |
goto L10;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* + END WHILE
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Now compute b = 2.0**m with the smallest positive integer
|
|
kusano |
7d535a |
m
|
|
kusano |
7d535a |
such that
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
fl( a + b ) .gt. a. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
b = 1.;
|
|
kusano |
7d535a |
c = dlamc3_(&a, &b);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* + WHILE( C.EQ.A )LOOP */
|
|
kusano |
7d535a |
L20:
|
|
kusano |
7d535a |
if (c == a) {
|
|
kusano |
7d535a |
b *= 2;
|
|
kusano |
7d535a |
c = dlamc3_(&a, &b);
|
|
kusano |
7d535a |
goto L20;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* + END WHILE
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Now compute the base. a and c are neighbouring floating po
|
|
kusano |
7d535a |
int
|
|
kusano |
7d535a |
numbers in the interval ( beta**t, beta**( t + 1 ) ) and
|
|
kusano |
7d535a |
so
|
|
kusano |
7d535a |
their difference is beta. Adding 0.25 to c is to ensure that
|
|
kusano |
7d535a |
it
|
|
kusano |
7d535a |
is truncated to beta and not ( beta - 1 ). */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
qtr = one / 4;
|
|
kusano |
7d535a |
savec = c;
|
|
kusano |
7d535a |
d__1 = -a;
|
|
kusano |
7d535a |
c = dlamc3_(&c, &d__1);
|
|
kusano |
7d535a |
lbeta = (int) (c + qtr);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Now determine whether rounding or chopping occurs, by addin
|
|
kusano |
7d535a |
g a
|
|
kusano |
7d535a |
bit less than beta/2 and a bit more than beta/2 to
|
|
kusano |
7d535a |
a. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
b = (double) lbeta;
|
|
kusano |
7d535a |
d__1 = b / 2;
|
|
kusano |
7d535a |
d__2 = -b / 100;
|
|
kusano |
7d535a |
f = dlamc3_(&d__1, &d__2);
|
|
kusano |
7d535a |
c = dlamc3_(&f, &a);
|
|
kusano |
7d535a |
if (c == a) {
|
|
kusano |
7d535a |
lrnd = TRUE_;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
lrnd = FALSE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
d__1 = b / 2;
|
|
kusano |
7d535a |
d__2 = b / 100;
|
|
kusano |
7d535a |
f = dlamc3_(&d__1, &d__2);
|
|
kusano |
7d535a |
c = dlamc3_(&f, &a);
|
|
kusano |
7d535a |
if (lrnd && c == a) {
|
|
kusano |
7d535a |
lrnd = FALSE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Try and decide whether rounding is done in the IEEE 'round
|
|
kusano |
7d535a |
to
|
|
kusano |
7d535a |
nearest' style. B/2 is half a unit in the last place of the
|
|
kusano |
7d535a |
two
|
|
kusano |
7d535a |
numbers A and SAVEC. Furthermore, A is even, i.e. has last
|
|
kusano |
7d535a |
bit
|
|
kusano |
7d535a |
zero, and SAVEC is odd. Thus adding B/2 to A should not cha
|
|
kusano |
7d535a |
nge
|
|
kusano |
7d535a |
A, but adding B/2 to SAVEC should change SAVEC. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
d__1 = b / 2;
|
|
kusano |
7d535a |
t1 = dlamc3_(&d__1, &a);
|
|
kusano |
7d535a |
d__1 = b / 2;
|
|
kusano |
7d535a |
t2 = dlamc3_(&d__1, &savec);
|
|
kusano |
7d535a |
lieee1 = t1 == a && t2 > savec && lrnd;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Now find the mantissa, t. It should be the integer part
|
|
kusano |
7d535a |
of
|
|
kusano |
7d535a |
log to the base beta of a, however it is safer to determine
|
|
kusano |
7d535a |
t
|
|
kusano |
7d535a |
by powering. So we find t as the smallest positive integer
|
|
kusano |
7d535a |
for
|
|
kusano |
7d535a |
which
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
fl( beta**t + 1.0 ) = 1.0. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
lt = 0;
|
|
kusano |
7d535a |
a = 1.;
|
|
kusano |
7d535a |
c = 1.;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* + WHILE( C.EQ.ONE )LOOP */
|
|
kusano |
7d535a |
L30:
|
|
kusano |
7d535a |
if (c == one) {
|
|
kusano |
7d535a |
++lt;
|
|
kusano |
7d535a |
a *= lbeta;
|
|
kusano |
7d535a |
c = dlamc3_(&a, &one);
|
|
kusano |
7d535a |
d__1 = -a;
|
|
kusano |
7d535a |
c = dlamc3_(&c, &d__1);
|
|
kusano |
7d535a |
goto L30;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* + END WHILE */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
*beta = lbeta;
|
|
kusano |
7d535a |
*t = lt;
|
|
kusano |
7d535a |
*rnd = lrnd;
|
|
kusano |
7d535a |
*ieee1 = lieee1;
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* End of DLAMC1 */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} /* dlamc1_ */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Subroutine */ int dlamc2_(int *beta, int *t, int *rnd,
|
|
kusano |
7d535a |
double *eps, int *emin, double *rmin, int *emax,
|
|
kusano |
7d535a |
double *rmax)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
/* -- LAPACK auxiliary routine (version 2.0) --
|
|
kusano |
7d535a |
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
kusano |
7d535a |
Courant Institute, Argonne National Lab, and Rice University
|
|
kusano |
7d535a |
October 31, 1992
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Purpose
|
|
kusano |
7d535a |
=======
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DLAMC2 determines the machine parameters specified in its argument
|
|
kusano |
7d535a |
list.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Arguments
|
|
kusano |
7d535a |
=========
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
BETA (output) INT
|
|
kusano |
7d535a |
The base of the machine.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
T (output) INT
|
|
kusano |
7d535a |
The number of ( BETA ) digits in the mantissa.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
RND (output) INT
|
|
kusano |
7d535a |
Specifies whether proper rounding ( RND = .TRUE. ) or
|
|
kusano |
7d535a |
chopping ( RND = .FALSE. ) occurs in addition. This may not
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
be a reliable guide to the way in which the machine performs
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
its arithmetic.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
EPS (output) DOUBLE PRECISION
|
|
kusano |
7d535a |
The smallest positive number such that
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
fl( 1.0 - EPS ) .LT. 1.0,
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
where fl denotes the computed value.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
EMIN (output) INT
|
|
kusano |
7d535a |
The minimum exponent before (gradual) underflow occurs.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
RMIN (output) DOUBLE PRECISION
|
|
kusano |
7d535a |
The smallest normalized number for the machine, given by
|
|
kusano |
7d535a |
BASE**( EMIN - 1 ), where BASE is the floating point value
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
of BETA.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
EMAX (output) INT
|
|
kusano |
7d535a |
The maximum exponent before overflow occurs.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
RMAX (output) DOUBLE PRECISION
|
|
kusano |
7d535a |
The largest positive number for the machine, given by
|
|
kusano |
7d535a |
BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
value of BETA.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Further Details
|
|
kusano |
7d535a |
===============
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
The computation of EPS is based on a routine PARANOIA by
|
|
kusano |
7d535a |
W. Kahan of the University of California at Berkeley.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
=====================================================================
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
/* Table of constant values */
|
|
kusano |
7d535a |
static int c__1 = 1;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Initialized data */
|
|
kusano |
7d535a |
static int first = TRUE_;
|
|
kusano |
7d535a |
static int iwarn = FALSE_;
|
|
kusano |
7d535a |
/* System generated locals */
|
|
kusano |
7d535a |
int i__1;
|
|
kusano |
7d535a |
double d__1, d__2, d__3, d__4, d__5;
|
|
kusano |
7d535a |
/* Builtin functions */
|
|
kusano |
7d535a |
double pow_di(double *, int *);
|
|
kusano |
7d535a |
/* Local variables */
|
|
kusano |
7d535a |
static int ieee;
|
|
kusano |
7d535a |
static double half;
|
|
kusano |
7d535a |
static int lrnd;
|
|
kusano |
7d535a |
static double leps, zero, a, b, c;
|
|
kusano |
7d535a |
static int i, lbeta;
|
|
kusano |
7d535a |
static double rbase;
|
|
kusano |
7d535a |
static int lemin, lemax, gnmin;
|
|
kusano |
7d535a |
static double small;
|
|
kusano |
7d535a |
static int gpmin;
|
|
kusano |
7d535a |
static double third, lrmin, lrmax, sixth;
|
|
kusano |
7d535a |
extern /* Subroutine */ int dlamc1_(int *, int *, int *,
|
|
kusano |
7d535a |
int *);
|
|
kusano |
7d535a |
extern double dlamc3_(double *, double *);
|
|
kusano |
7d535a |
static int lieee1;
|
|
kusano |
7d535a |
extern /* Subroutine */ int dlamc4_(int *, double *, int *),
|
|
kusano |
7d535a |
dlamc5_(int *, int *, int *, int *, int *,
|
|
kusano |
7d535a |
double *);
|
|
kusano |
7d535a |
static int lt, ngnmin, ngpmin;
|
|
kusano |
7d535a |
static double one, two;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (first) {
|
|
kusano |
7d535a |
first = FALSE_;
|
|
kusano |
7d535a |
zero = 0.;
|
|
kusano |
7d535a |
one = 1.;
|
|
kusano |
7d535a |
two = 2.;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values
|
|
kusano |
7d535a |
of
|
|
kusano |
7d535a |
BETA, T, RND, EPS, EMIN and RMIN.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Throughout this routine we use the function DLAMC3 to ens
|
|
kusano |
7d535a |
ure
|
|
kusano |
7d535a |
that relevant values are stored and not held in registers,
|
|
kusano |
7d535a |
or
|
|
kusano |
7d535a |
are not affected by optimizers.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
dlamc1_(&lbeta, <, &lrnd, &lieee1);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Start to find EPS. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
b = (double) lbeta;
|
|
kusano |
7d535a |
i__1 = -lt;
|
|
kusano |
7d535a |
a = pow_di(&b, &i__1);
|
|
kusano |
7d535a |
leps = a;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Try some tricks to see whether or not this is the correct E
|
|
kusano |
7d535a |
PS. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
b = two / 3;
|
|
kusano |
7d535a |
half = one / 2;
|
|
kusano |
7d535a |
d__1 = -half;
|
|
kusano |
7d535a |
sixth = dlamc3_(&b, &d__1);
|
|
kusano |
7d535a |
third = dlamc3_(&sixth, &sixth);
|
|
kusano |
7d535a |
d__1 = -half;
|
|
kusano |
7d535a |
b = dlamc3_(&third, &d__1);
|
|
kusano |
7d535a |
b = dlamc3_(&b, &sixth);
|
|
kusano |
7d535a |
b = abs(b);
|
|
kusano |
7d535a |
if (b < leps) {
|
|
kusano |
7d535a |
b = leps;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
leps = 1.;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
|
|
kusano |
7d535a |
L10:
|
|
kusano |
7d535a |
if (leps > b && b > zero) {
|
|
kusano |
7d535a |
leps = b;
|
|
kusano |
7d535a |
d__1 = half * leps;
|
|
kusano |
7d535a |
/* Computing 5th power */
|
|
kusano |
7d535a |
d__3 = two, d__4 = d__3, d__3 *= d__3;
|
|
kusano |
7d535a |
/* Computing 2nd power */
|
|
kusano |
7d535a |
d__5 = leps;
|
|
kusano |
7d535a |
d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
|
|
kusano |
7d535a |
c = dlamc3_(&d__1, &d__2);
|
|
kusano |
7d535a |
d__1 = -c;
|
|
kusano |
7d535a |
c = dlamc3_(&half, &d__1);
|
|
kusano |
7d535a |
b = dlamc3_(&half, &c);
|
|
kusano |
7d535a |
d__1 = -b;
|
|
kusano |
7d535a |
c = dlamc3_(&half, &d__1);
|
|
kusano |
7d535a |
b = dlamc3_(&half, &c);
|
|
kusano |
7d535a |
goto L10;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* + END WHILE */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (a < leps) {
|
|
kusano |
7d535a |
leps = a;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Computation of EPS complete.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3
|
|
kusano |
7d535a |
)).
|
|
kusano |
7d535a |
Keep dividing A by BETA until (gradual) underflow occurs. T
|
|
kusano |
7d535a |
his
|
|
kusano |
7d535a |
is detected when we cannot recover the previous A. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
rbase = one / lbeta;
|
|
kusano |
7d535a |
small = one;
|
|
kusano |
7d535a |
for (i = 1; i <= 3; ++i) {
|
|
kusano |
7d535a |
d__1 = small * rbase;
|
|
kusano |
7d535a |
small = dlamc3_(&d__1, &zero);
|
|
kusano |
7d535a |
/* L20: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
a = dlamc3_(&one, &small);
|
|
kusano |
7d535a |
dlamc4_(&ngpmin, &one, &lbeta);
|
|
kusano |
7d535a |
d__1 = -one;
|
|
kusano |
7d535a |
dlamc4_(&ngnmin, &d__1, &lbeta);
|
|
kusano |
7d535a |
dlamc4_(&gpmin, &a, &lbeta);
|
|
kusano |
7d535a |
d__1 = -a;
|
|
kusano |
7d535a |
dlamc4_(&gnmin, &d__1, &lbeta);
|
|
kusano |
7d535a |
ieee = FALSE_;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (ngpmin == ngnmin && gpmin == gnmin) {
|
|
kusano |
7d535a |
if (ngpmin == gpmin) {
|
|
kusano |
7d535a |
lemin = ngpmin;
|
|
kusano |
7d535a |
/* ( Non twos-complement machines, no gradual under
|
|
kusano |
7d535a |
flow;
|
|
kusano |
7d535a |
e.g., VAX ) */
|
|
kusano |
7d535a |
} else if (gpmin - ngpmin == 3) {
|
|
kusano |
7d535a |
lemin = ngpmin - 1 + lt;
|
|
kusano |
7d535a |
ieee = TRUE_;
|
|
kusano |
7d535a |
/* ( Non twos-complement machines, with gradual und
|
|
kusano |
7d535a |
erflow;
|
|
kusano |
7d535a |
e.g., IEEE standard followers ) */
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
lemin = min(ngpmin,gpmin);
|
|
kusano |
7d535a |
/* ( A guess; no known machine ) */
|
|
kusano |
7d535a |
iwarn = TRUE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} else if (ngpmin == gpmin && ngnmin == gnmin) {
|
|
kusano |
7d535a |
if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
|
|
kusano |
7d535a |
lemin = max(ngpmin,ngnmin);
|
|
kusano |
7d535a |
/* ( Twos-complement machines, no gradual underflow
|
|
kusano |
7d535a |
;
|
|
kusano |
7d535a |
e.g., CYBER 205 ) */
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
lemin = min(ngpmin,ngnmin);
|
|
kusano |
7d535a |
/* ( A guess; no known machine ) */
|
|
kusano |
7d535a |
iwarn = TRUE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
if (gpmin - min(ngpmin,ngnmin) == 3) {
|
|
kusano |
7d535a |
lemin = max(ngpmin,ngnmin) - 1 + lt;
|
|
kusano |
7d535a |
/* ( Twos-complement machines with gradual underflo
|
|
kusano |
7d535a |
w;
|
|
kusano |
7d535a |
no known machine ) */
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
lemin = min(ngpmin,ngnmin);
|
|
kusano |
7d535a |
/* ( A guess; no known machine ) */
|
|
kusano |
7d535a |
iwarn = TRUE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
/* Computing MIN */
|
|
kusano |
7d535a |
i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin);
|
|
kusano |
7d535a |
lemin = min(i__1,gnmin);
|
|
kusano |
7d535a |
/* ( A guess; no known machine ) */
|
|
kusano |
7d535a |
iwarn = TRUE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* **
|
|
kusano |
7d535a |
Comment out this if block if EMIN is ok */
|
|
kusano |
7d535a |
if (iwarn) {
|
|
kusano |
7d535a |
first = TRUE_;
|
|
kusano |
7d535a |
printf("\n\n WARNING. The value EMIN may be incorrect:- ");
|
|
kusano |
7d535a |
printf("EMIN = %8i\n",lemin);
|
|
kusano |
7d535a |
printf("If, after inspection, the value EMIN looks acceptable");
|
|
kusano |
7d535a |
printf("please comment out \n the IF block as marked within the");
|
|
kusano |
7d535a |
printf("code of routine DLAMC2, \n otherwise supply EMIN");
|
|
kusano |
7d535a |
printf("explicitly.\n");
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* **
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Assume IEEE arithmetic if we found denormalised numbers abo
|
|
kusano |
7d535a |
ve,
|
|
kusano |
7d535a |
or if arithmetic seems to round in the IEEE style, determi
|
|
kusano |
7d535a |
ned
|
|
kusano |
7d535a |
in routine DLAMC1. A true IEEE machine should have both thi
|
|
kusano |
7d535a |
ngs
|
|
kusano |
7d535a |
true; however, faulty machines may have one or the other. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
ieee = ieee || lieee1;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Compute RMIN by successive division by BETA. We could comp
|
|
kusano |
7d535a |
ute
|
|
kusano |
7d535a |
RMIN as BASE**( EMIN - 1 ), but some machines underflow dur
|
|
kusano |
7d535a |
ing
|
|
kusano |
7d535a |
this computation. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
lrmin = 1.;
|
|
kusano |
7d535a |
i__1 = 1 - lemin;
|
|
kusano |
7d535a |
for (i = 1; i <= 1-lemin; ++i) {
|
|
kusano |
7d535a |
d__1 = lrmin * rbase;
|
|
kusano |
7d535a |
lrmin = dlamc3_(&d__1, &zero);
|
|
kusano |
7d535a |
/* L30: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Finally, call DLAMC5 to compute EMAX and RMAX. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
dlamc5_(&lbeta, <, &lemin, &ieee, &lemax, &lrmax);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
*beta = lbeta;
|
|
kusano |
7d535a |
*t = lt;
|
|
kusano |
7d535a |
*rnd = lrnd;
|
|
kusano |
7d535a |
*eps = leps;
|
|
kusano |
7d535a |
*emin = lemin;
|
|
kusano |
7d535a |
*rmin = lrmin;
|
|
kusano |
7d535a |
*emax = lemax;
|
|
kusano |
7d535a |
*rmax = lrmax;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* End of DLAMC2 */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} /* dlamc2_ */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
double dlamc3_(double *a, double *b)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
/* -- LAPACK auxiliary routine (version 2.0) --
|
|
kusano |
7d535a |
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
kusano |
7d535a |
Courant Institute, Argonne National Lab, and Rice University
|
|
kusano |
7d535a |
October 31, 1992
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Purpose
|
|
kusano |
7d535a |
=======
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DLAMC3 is intended to force A and B to be stored prior to doing
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
the addition of A and B , for use in situations where optimizers
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
might hold one of these in a register.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Arguments
|
|
kusano |
7d535a |
=========
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
A, B (input) DOUBLE PRECISION
|
|
kusano |
7d535a |
The values A and B.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
=====================================================================
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
/* >>Start of File<<
|
|
kusano |
7d535a |
System generated locals */
|
|
kusano |
7d535a |
double ret_val;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
ret_val = *a + *b;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
return ret_val;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* End of DLAMC3 */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} /* dlamc3_ */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Subroutine */ int dlamc4_(int *emin, double *start, int *base)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
/* -- LAPACK auxiliary routine (version 2.0) --
|
|
kusano |
7d535a |
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
kusano |
7d535a |
Courant Institute, Argonne National Lab, and Rice University
|
|
kusano |
7d535a |
October 31, 1992
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Purpose
|
|
kusano |
7d535a |
=======
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DLAMC4 is a service routine for DLAMC2.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Arguments
|
|
kusano |
7d535a |
=========
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
EMIN (output) EMIN
|
|
kusano |
7d535a |
The minimum exponent before (gradual) underflow, computed by
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
setting A = START and dividing by BASE until the previous A
|
|
kusano |
7d535a |
can not be recovered.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
START (input) DOUBLE PRECISION
|
|
kusano |
7d535a |
The starting point for determining EMIN.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
BASE (input) INT
|
|
kusano |
7d535a |
The base of the machine.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
=====================================================================
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
/* System generated locals */
|
|
kusano |
7d535a |
int i__1;
|
|
kusano |
7d535a |
double d__1;
|
|
kusano |
7d535a |
/* Local variables */
|
|
kusano |
7d535a |
static double zero, a;
|
|
kusano |
7d535a |
static int i;
|
|
kusano |
7d535a |
static double rbase, b1, b2, c1, c2, d1, d2;
|
|
kusano |
7d535a |
extern double dlamc3_(double *, double *);
|
|
kusano |
7d535a |
static double one;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
a = *start;
|
|
kusano |
7d535a |
one = 1.;
|
|
kusano |
7d535a |
rbase = one / *base;
|
|
kusano |
7d535a |
zero = 0.;
|
|
kusano |
7d535a |
*emin = 1;
|
|
kusano |
7d535a |
d__1 = a * rbase;
|
|
kusano |
7d535a |
b1 = dlamc3_(&d__1, &zero);
|
|
kusano |
7d535a |
c1 = a;
|
|
kusano |
7d535a |
c2 = a;
|
|
kusano |
7d535a |
d1 = a;
|
|
kusano |
7d535a |
d2 = a;
|
|
kusano |
7d535a |
/* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
|
|
kusano |
7d535a |
$ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
|
|
kusano |
7d535a |
L10:
|
|
kusano |
7d535a |
if (c1 == a && c2 == a && d1 == a && d2 == a) {
|
|
kusano |
7d535a |
--(*emin);
|
|
kusano |
7d535a |
a = b1;
|
|
kusano |
7d535a |
d__1 = a / *base;
|
|
kusano |
7d535a |
b1 = dlamc3_(&d__1, &zero);
|
|
kusano |
7d535a |
d__1 = b1 * *base;
|
|
kusano |
7d535a |
c1 = dlamc3_(&d__1, &zero);
|
|
kusano |
7d535a |
d1 = zero;
|
|
kusano |
7d535a |
i__1 = *base;
|
|
kusano |
7d535a |
for (i = 1; i <= *base; ++i) {
|
|
kusano |
7d535a |
d1 += b1;
|
|
kusano |
7d535a |
/* L20: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
d__1 = a * rbase;
|
|
kusano |
7d535a |
b2 = dlamc3_(&d__1, &zero);
|
|
kusano |
7d535a |
d__1 = b2 / rbase;
|
|
kusano |
7d535a |
c2 = dlamc3_(&d__1, &zero);
|
|
kusano |
7d535a |
d2 = zero;
|
|
kusano |
7d535a |
i__1 = *base;
|
|
kusano |
7d535a |
for (i = 1; i <= *base; ++i) {
|
|
kusano |
7d535a |
d2 += b2;
|
|
kusano |
7d535a |
/* L30: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
goto L10;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* + END WHILE */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* End of DLAMC4 */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} /* dlamc4_ */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Subroutine */ int dlamc5_(int *beta, int *p, int *emin,
|
|
kusano |
7d535a |
int *ieee, int *emax, double *rmax)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
/* -- LAPACK auxiliary routine (version 2.0) --
|
|
kusano |
7d535a |
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
kusano |
7d535a |
Courant Institute, Argonne National Lab, and Rice University
|
|
kusano |
7d535a |
October 31, 1992
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Purpose
|
|
kusano |
7d535a |
=======
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DLAMC5 attempts to compute RMAX, the largest machine floating-point
|
|
kusano |
7d535a |
number, without overflow. It assumes that EMAX + abs(EMIN) sum
|
|
kusano |
7d535a |
approximately to a power of 2. It will fail on machines where this
|
|
kusano |
7d535a |
assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
EMAX = 28718). It will also fail if the value supplied for EMIN is
|
|
kusano |
7d535a |
too large (i.e. too close to zero), probably with overflow.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Arguments
|
|
kusano |
7d535a |
=========
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
BETA (input) INT
|
|
kusano |
7d535a |
The base of floating-point arithmetic.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
P (input) INT
|
|
kusano |
7d535a |
The number of base BETA digits in the mantissa of a
|
|
kusano |
7d535a |
floating-point value.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
EMIN (input) INT
|
|
kusano |
7d535a |
The minimum exponent before (gradual) underflow.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
IEEE (input) INT
|
|
kusano |
7d535a |
A int flag specifying whether or not the arithmetic
|
|
kusano |
7d535a |
system is thought to comply with the IEEE standard.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
EMAX (output) INT
|
|
kusano |
7d535a |
The largest exponent before overflow
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
RMAX (output) DOUBLE PRECISION
|
|
kusano |
7d535a |
The largest machine floating-point number.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
=====================================================================
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
First compute LEXP and UEXP, two powers of 2 that bound
|
|
kusano |
7d535a |
abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
|
|
kusano |
7d535a |
approximately to the bound that is closest to abs(EMIN).
|
|
kusano |
7d535a |
(EMAX is the exponent of the required number RMAX). */
|
|
kusano |
7d535a |
/* Table of constant values */
|
|
kusano |
7d535a |
static double c_b5 = 0.;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* System generated locals */
|
|
kusano |
7d535a |
int i__1;
|
|
kusano |
7d535a |
double d__1;
|
|
kusano |
7d535a |
/* Local variables */
|
|
kusano |
7d535a |
static int lexp;
|
|
kusano |
7d535a |
static double oldy;
|
|
kusano |
7d535a |
static int uexp, i;
|
|
kusano |
7d535a |
static double y, z;
|
|
kusano |
7d535a |
static int nbits;
|
|
kusano |
7d535a |
extern double dlamc3_(double *, double *);
|
|
kusano |
7d535a |
static double recbas;
|
|
kusano |
7d535a |
static int exbits, expsum, try__;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
lexp = 1;
|
|
kusano |
7d535a |
exbits = 1;
|
|
kusano |
7d535a |
L10:
|
|
kusano |
7d535a |
try__ = lexp << 1;
|
|
kusano |
7d535a |
if (try__ <= -(*emin)) {
|
|
kusano |
7d535a |
lexp = try__;
|
|
kusano |
7d535a |
++exbits;
|
|
kusano |
7d535a |
goto L10;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
if (lexp == -(*emin)) {
|
|
kusano |
7d535a |
uexp = lexp;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
uexp = try__;
|
|
kusano |
7d535a |
++exbits;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
|
|
kusano |
7d535a |
than or equal to EMIN. EXBITS is the number of bits needed to
|
|
kusano |
7d535a |
store the exponent. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (uexp + *emin > -lexp - *emin) {
|
|
kusano |
7d535a |
expsum = lexp << 1;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
expsum = uexp << 1;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* EXPSUM is the exponent range, approximately equal to
|
|
kusano |
7d535a |
EMAX - EMIN + 1 . */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
*emax = expsum + *emin - 1;
|
|
kusano |
7d535a |
nbits = exbits + 1 + *p;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* NBITS is the total number of bits needed to store a
|
|
kusano |
7d535a |
floating-point number. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (nbits % 2 == 1 && *beta == 2) {
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Either there are an odd number of bits used to store a
|
|
kusano |
7d535a |
floating-point number, which is unlikely, or some bits are
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
not used in the representation of numbers, which is possible
|
|
kusano |
7d535a |
,
|
|
kusano |
7d535a |
(e.g. Cray machines) or the mantissa has an implicit bit,
|
|
kusano |
7d535a |
(e.g. IEEE machines, Dec Vax machines), which is perhaps the
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
most likely. We have to assume the last alternative.
|
|
kusano |
7d535a |
If this is true, then we need to reduce EMAX by one because
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
there must be some way of representing zero in an implicit-b
|
|
kusano |
7d535a |
it
|
|
kusano |
7d535a |
system. On machines like Cray, we are reducing EMAX by one
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
unnecessarily. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
--(*emax);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (*ieee) {
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Assume we are on an IEEE machine which reserves one exponent
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
for infinity and NaN. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
--(*emax);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Now create RMAX, the largest machine number, which should
|
|
kusano |
7d535a |
be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
First compute 1.0 - BETA**(-P), being careful that the
|
|
kusano |
7d535a |
result is less than 1.0 . */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
recbas = 1. / *beta;
|
|
kusano |
7d535a |
z = *beta - 1.;
|
|
kusano |
7d535a |
y = 0.;
|
|
kusano |
7d535a |
i__1 = *p;
|
|
kusano |
7d535a |
for (i = 1; i <= *p; ++i) {
|
|
kusano |
7d535a |
z *= recbas;
|
|
kusano |
7d535a |
if (y < 1.) {
|
|
kusano |
7d535a |
oldy = y;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
y = dlamc3_(&y, &z);
|
|
kusano |
7d535a |
/* L20: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
if (y >= 1.) {
|
|
kusano |
7d535a |
y = oldy;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Now multiply by BETA**EMAX to get RMAX. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
i__1 = *emax;
|
|
kusano |
7d535a |
for (i = 1; i <= *emax; ++i) {
|
|
kusano |
7d535a |
d__1 = y * *beta;
|
|
kusano |
7d535a |
y = dlamc3_(&d__1, &c_b5);
|
|
kusano |
7d535a |
/* L30: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
*rmax = y;
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* End of DLAMC5 */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} /* dlamc5_ */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
double pow_di(double *ap, int *bp)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
double pow, x;
|
|
kusano |
7d535a |
int n;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
pow = 1;
|
|
kusano |
7d535a |
x = *ap;
|
|
kusano |
7d535a |
n = *bp;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if(n != 0){
|
|
kusano |
7d535a |
if(n < 0) {
|
|
kusano |
7d535a |
n = -n;
|
|
kusano |
7d535a |
x = 1/x;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
for( ; ; ) {
|
|
kusano |
7d535a |
if(n & 01) pow *= x;
|
|
kusano |
7d535a |
if(n >>= 1) x *= x;
|
|
kusano |
7d535a |
else break;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
return(pow);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|