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