Mercurial > hg > octave-max
changeset 2425:679068a18eee
[project @ 1996-10-25 01:24:59 by jwe]
| author | jwe |
|---|---|
| date | Fri, 25 Oct 1996 01:26:47 +0000 |
| parents | b5c3b08f1bab |
| children | 1b5536a0bbb4 |
| files | libcruft/ChangeLog libcruft/Makefile.in libcruft/misc/Makefile.in libcruft/misc/d1mach.f libcruft/misc/gen-d1mach.c libcruft/misc/machar.c liboctave/Quad.h |
| diffstat | 7 files changed, 414 insertions(+), 429 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,3 +1,10 @@ +Thu Oct 24 20:22:47 1996 John W. Eaton <jwe@bevo.che.wisc.edu> + + * Makefile.in (CRUFT_OBJ): No special treatment for d1mach.o. + + * misc/machar.c, misc/d1mach.f: New files + * misc/Makefile.in: Fix to not generate d1mach.f. + Mon Oct 14 11:07:25 1996 John W. Eaton <jwe@bevo.che.wisc.edu> * Makefile.in (distclean): Remove stamp-shared too.
--- a/libcruft/Makefile.in +++ b/libcruft/Makefile.in @@ -40,12 +40,11 @@ cd $@; $(MAKE) all .PHONY: $(SUBDIRS) -MISC_OBJ := misc/d1mach.o misc/dostop.o misc/f77-extern.o misc/lo-error.o +MISC_OBJ := misc/machar.o misc/dostop.o misc/f77-extern.o misc/lo-error.o CRUFT_FSRC := $(foreach dir, $(SUBDIRS), $(wildcard $(srcdir)/$(dir)/*.f)) -CRUFT_OBJ3 := $(patsubst $(srcdir)/%, %, $(CRUFT_FSRC)) -CRUFT_OBJ2 := $(patsubst %.f, %.o, $(CRUFT_OBJ3)) -CRUFT_OBJ1 := $(subst misc/d1mach.o, , $(CRUFT_OBJ2)) +CRUFT_OBJ2 := $(patsubst $(srcdir)/%, %, $(CRUFT_FSRC)) +CRUFT_OBJ1 := $(patsubst %.f, %.o, $(CRUFT_OBJ2)) CRUFT_OBJ := $(CRUFT_OBJ1) $(MISC_OBJ) ifeq ($(SHARED_LIBS), true)
--- a/libcruft/misc/Makefile.in +++ b/libcruft/misc/Makefile.in @@ -12,11 +12,9 @@ top_srcdir = @top_srcdir@ VPATH = @srcdir@ -SPECIAL := gen-d1mach.c d1mach-tst.for dostop.c f77-extern.cc lo-error.cc +SPECIAL := machar.c d1mach-tst.for dostop.c f77-extern.cc lo-error.cc -SPECIAL_DEPEND := d1mach.o dostop.o f77-extern.o lo-error.o - -DISTFILES := $(subst d1mach.f, , $(DISTFILES)) +SPECIAL_DEPEND := d1mach.o machar.o dostop.o f77-extern.o lo-error.o EXTERNAL_DISTFILES = $(DISTFILES) @@ -34,15 +32,13 @@ include ../Makerules -d1mach.f: gen-d1mach - ./gen-d1mach > d1mach.f - # Don't optimize. XCC = $(patsubst -O%, , $(CC)) +XALL_CFLAGS = $(patsubst -O%, , $(ALL_CFLAGS)) -gen-d1mach: $(srcdir)/gen-d1mach.c - $(XCC) -DDP -o gen-d1mach $(srcdir)/gen-d1mach.c -lm +machar.o: $(srcdir)/machar.c + $(XCC) -c $(CPP_FLAGS) $(XALL_CFLAGS) -DDP $(srcdir)/machar.c distclean:: rm -f d1mach.f gen-d1mach
new file mode 100644 --- /dev/null +++ b/libcruft/misc/d1mach.f @@ -0,0 +1,18 @@ + double precision function d1mach (i) + integer i + logical init + double precision dmach(5) + data init /.false./ + save init, dmach + if (.not. init) then + call machar (dmach(1), dmach(2), dmach(3), dmach(4), dmach(5)) + init = .true. + endif + if (i .lt. 1 .or. i .gt. 5) goto 999 + d1mach = dmach(i) + return + 999 write(*,1999) i + 1999 format(' d1mach - i out of bounds', i10) + call xstopx (' ') + d1mach = 0 + end
deleted file mode 100644 --- a/libcruft/misc/gen-d1mach.c +++ /dev/null @@ -1,415 +0,0 @@ -/* - -This file combines the single and double precision versions of machar, -selected by cc -DSP or cc -DDP. This feature provided by D. G. Hough, -August 3, 1988. - -*/ - -#ifdef SP -#define REAL float -#define ZERO 0.0 -#define ONE 1.0 -#define PREC "Single " -#define REALSIZE 1 -#endif - -#ifdef DP -#define REAL double -#define ZERO 0.0e0 -#define ONE 1.0e0 -#define PREC "Double " -#define REALSIZE 2 -#endif - -#include <math.h> -#include <stdio.h> - -#define ABS(xxx) ((xxx>ZERO)?(xxx):(-xxx)) - -void -rmachar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp, - maxexp,eps,epsneg,xmin,xmax) - - int *ibeta,*iexp,*irnd,*it,*machep,*maxexp,*minexp,*negep,*ngrd; - REAL *eps,*epsneg,*xmax,*xmin; - -/* - - This subroutine is intended to determine the parameters of the - floating-point arithmetic system specified below. The - determination of the first three uses an extension of an algorithm - due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some, - but not all, of the improvements suggested by M. Gentleman and S. - Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this - program was published in the book Software Manual for the - Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall, - Englewood Cliffs, NJ, 1980. The present program is a - translation of the Fortran 77 program in W. J. Cody, "MACHAR: - A subroutine to dynamically determine machine parameters". - TOMS (14), 1988. - - Parameter values reported are as follows: - - ibeta - the radix for the floating-point representation - it - the number of base ibeta digits in the floating-point - significand - irnd - 0 if floating-point addition chops - 1 if floating-point addition rounds, but not in the - IEEE style - 2 if floating-point addition rounds in the IEEE style - 3 if floating-point addition chops, and there is - partial underflow - 4 if floating-point addition rounds, but not in the - IEEE style, and there is partial underflow - 5 if floating-point addition rounds in the IEEE style, - and there is partial underflow - ngrd - the number of guard digits for multiplication with - truncating arithmetic. It is - 0 if floating-point arithmetic rounds, or if it - truncates and only it base ibeta digits - participate in the post-normalization shift of the - floating-point significand in multiplication; - 1 if floating-point arithmetic truncates and more - than it base ibeta digits participate in the - post-normalization shift of the floating-point - significand in multiplication. - machep - the largest negative integer such that - 1.0+FLOAT(ibeta)**machep .NE. 1.0, except that - machep is bounded below by -(it+3) - negeps - the largest negative integer such that - 1.0-FLOAT(ibeta)**negeps .NE. 1.0, except that - negeps is bounded below by -(it+3) - iexp - the number of bits (decimal places if ibeta = 10) - reserved for the representation of the exponent - (including the bias or sign) of a floating-point - number - minexp - the largest in magnitude negative integer such that - FLOAT(ibeta)**minexp is positive and normalized - maxexp - the smallest positive power of BETA that overflows - eps - the smallest positive floating-point number such - that 1.0+eps .NE. 1.0. In particular, if either - ibeta = 2 or IRND = 0, eps = FLOAT(ibeta)**machep. - Otherwise, eps = (FLOAT(ibeta)**machep)/2 - epsneg - A small positive floating-point number such that - 1.0-epsneg .NE. 1.0. In particular, if ibeta = 2 - or IRND = 0, epsneg = FLOAT(ibeta)**negeps. - Otherwise, epsneg = (ibeta**negeps)/2. Because - negeps is bounded below by -(it+3), epsneg may not - be the smallest number that can alter 1.0 by - subtraction. - xmin - the smallest non-vanishing normalized floating-point - power of the radix, i.e., xmin = FLOAT(ibeta)**minexp - xmax - the largest finite floating-point number. In - particular xmax = (1.0-epsneg)*FLOAT(ibeta)**maxexp - Note - on some machines xmax will be only the - second, or perhaps third, largest number, being - too small by 1 or 2 units in the last digit of - the significand. - - Latest revision - August 4, 1988 - - Author - W. J. Cody - Argonne National Laboratory - -*/ - -{ - int i,iz,j,k; - int mx,itmp,nxres; - REAL a,b,beta,betain,one,y,z,zero; - REAL betah,t,tmp,tmpa,tmp1,two; - - (*irnd) = 1; - one = (REAL)(*irnd); - two = one + one; - a = two; - b = a; - zero = 0.0e0; - -/* - determine ibeta,beta ala malcolm -*/ - - tmp = ((a+one)-a)-one; - - while (tmp == zero) { - a = a+a; - tmp = a+one; - tmp1 = tmp-a; - tmp = tmp1-one; - } - - tmp = a+b; - itmp = (int)(tmp-a); - while (itmp == 0) { - b = b+b; - tmp = a+b; - itmp = (int)(tmp-a); - } - - *ibeta = itmp; - beta = (REAL)(*ibeta); - -/* - determine irnd, it -*/ - - (*it) = 0; - b = one; - tmp = ((b+one)-b)-one; - - while (tmp == zero) { - *it = *it+1; - b = b*beta; - tmp = b+one; - tmp1 = tmp-b; - tmp = tmp1-one; - } - - *irnd = 0; - betah = beta/two; - tmp = a+betah; - tmp1 = tmp-a; - if (tmp1 != zero) *irnd = 1; - tmpa = a+beta; - tmp = tmpa+betah; - if ((*irnd == 0) && (tmp-tmpa != zero)) *irnd = 2; - -/* - determine negep, epsneg -*/ - - (*negep) = (*it) + 3; - betain = one / beta; - a = one; - - for (i = 1; i<=(*negep); i++) { - a = a * betain; - } - - b = a; - tmp = (one-a); - tmp = tmp-one; - - while (tmp == zero) { - a = a*beta; - *negep = *negep-1; - tmp1 = one-a; - tmp = tmp1-one; - } - - (*negep) = -(*negep); - (*epsneg) = a; - -/* - determine machep, eps -*/ - - (*machep) = -(*it) - 3; - a = b; - tmp = one+a; - - while (tmp-one == zero) { - a = a*beta; - *machep = *machep+1; - tmp = one+a; - } - - *eps = a; - -/* - determine ngrd -*/ - - (*ngrd) = 0; - tmp = one+*eps; - tmp = tmp*one; - if (((*irnd) == 0) && (tmp-one) != zero) (*ngrd) = 1; - -/* - determine iexp, minexp, xmin - - loop to determine largest i such that - (1/beta) ** (2**(i)) - does not underflow. - exit from loop is signaled by an underflow. -*/ - - i = 0; - k = 1; - z = betain; - t = one+*eps; - nxres = 0; - - for (;;) { - y = z; - z = y * y; - -/* - check for underflow -*/ - - a = z * one; - tmp = z*t; - if ((a+a == zero) || (ABS(z) > y)) break; - tmp1 = tmp*betain; - if (tmp1*beta == z) break; - i = i + 1; - k = k+k; - } - -/* - determine k such that (1/beta)**k does not underflow - first set k = 2 ** i -*/ - - (*iexp) = i + 1; - mx = k + k; - if (*ibeta == 10) { - -/* - for decimal machines only -*/ - - (*iexp) = 2; - iz = *ibeta; - while (k >= iz) { - iz = iz * (*ibeta); - (*iexp) = (*iexp) + 1; - } - mx = iz + iz - 1; - } - -/* - loop to determine minexp, xmin. - exit from loop is signaled by an underflow. -*/ - - for (;;) { - (*xmin) = y; - y = y * betain; - a = y * one; - tmp = y*t; - tmp1 = a+a; - if ((tmp1 == zero) || (ABS(y) >= (*xmin))) break; - k = k + 1; - tmp1 = tmp*betain; - tmp1 = tmp1*beta; - - if ((tmp1 == y) && (tmp != y)) { - nxres = 3; - *xmin = y; - break; - } - - } - - (*minexp) = -k; - -/* - determine maxexp, xmax -*/ - - if ((mx <= k+k-3) && ((*ibeta) != 10)) { - mx = mx + mx; - (*iexp) = (*iexp) + 1; - } - - (*maxexp) = mx + (*minexp); - -/* - Adjust *irnd to reflect partial underflow. -*/ - - (*irnd) = (*irnd)+nxres; - -/* - Adjust for IEEE style machines. -*/ - - if ((*irnd) >= 2) (*maxexp) = (*maxexp)-2; - -/* - adjust for machines with implicit leading bit in binary - significand and machines with radix point at extreme - right of significand. -*/ - - i = (*maxexp) + (*minexp); - if (((*ibeta) == 2) && (i == 0)) (*maxexp) = (*maxexp) - 1; - if (i > 20) (*maxexp) = (*maxexp) - 1; - if (a != y) (*maxexp) = (*maxexp) - 2; - (*xmax) = one - (*epsneg); - tmp = (*xmax)*one; - if (tmp != (*xmax)) (*xmax) = one - beta * (*epsneg); - (*xmax) = (*xmax) / (beta * beta * beta * (*xmin)); - i = (*maxexp) + (*minexp) + 3; - if (i > 0) { - - for (j = 1; j<=i; j++ ) { - if ((*ibeta) == 2) (*xmax) = (*xmax) + (*xmax); - if ((*ibeta) != 2) (*xmax) = (*xmax) * beta; - } - - } - - return; - -} - -typedef union -{ - double d; - int i[2]; -} equiv; - -int -main (void) -{ - /* Works for 32 bit machines with 32 bit ints and 64 bit doubles */ - - int ibeta, iexp, irnd, it, machep, maxexp, minexp, negep, ngrd; - REAL eps, epsneg, xmax, xmin; - int i; - equiv flt_params[6]; - - rmachar (&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp, - &maxexp, &eps, &epsneg, &xmin, &xmax); - - flt_params[1].d = xmin; - flt_params[2].d = xmax; - flt_params[3].d = epsneg; - flt_params[4].d = eps; - flt_params[5].d = log10 ((double) ibeta); - - printf ("* d1mach.f Do not edit. Generated automatically by gen-d1mach.c\n\ - double precision function d1mach(i)\n\ - integer i\n\ - integer i1var (2)\n\ - integer i2var (2)\n\ - integer i3var (2)\n\ - integer i4var (2)\n\ - integer i5var (2)\n\ - double precision dmach(5)\n\ - equivalence (dmach(1), i1var(1))\n\ - equivalence (dmach(2), i2var(1))\n\ - equivalence (dmach(3), i3var(1))\n\ - equivalence (dmach(4), i4var(1))\n\ - equivalence (dmach(5), i5var(1))\n"); - - for (i = 1; i < 6; i++) - printf (" data i%dvar(1), i%dvar(2) / %ld , %ld /\n", - i, i, flt_params[i].i[0], flt_params[i].i[1]); - - printf (" if (i .lt. 1 .or. i .gt. 5) goto 999\n\ - d1mach = dmach(i)\n\ - return\n\ - 999 write(*,1999) i\n\ - 1999 format(' d1mach - i out of bounds', i10)\n\ - call xstopx (' ')\n\ - d1mach = 0\n\ - end\n"); - - return 0; -}
new file mode 100644 --- /dev/null +++ b/libcruft/misc/machar.c @@ -0,0 +1,380 @@ +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "f77-fcn.h" + +/* + +This file combines the single and double precision versions of machar, +selected by cc -DSP or cc -DDP. This feature provided by D. G. Hough, +August 3, 1988. + +*/ + +#ifdef SP +#define REAL float +#define ZERO 0.0 +#define ONE 1.0 +#define PREC "Single " +#define REALSIZE 1 +#endif + +#ifdef DP +#define REAL double +#define ZERO 0.0e0 +#define ONE 1.0e0 +#define PREC "Double " +#define REALSIZE 2 +#endif + +#include <math.h> +#include <stdio.h> + +#define ABS(xxx) ((xxx>ZERO)?(xxx):(-xxx)) + +void +rmachar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp, + maxexp,eps,epsneg,xmin,xmax) + + int *ibeta,*iexp,*irnd,*it,*machep,*maxexp,*minexp,*negep,*ngrd; + REAL *eps,*epsneg,*xmax,*xmin; + +/* + + This subroutine is intended to determine the parameters of the + floating-point arithmetic system specified below. The + determination of the first three uses an extension of an algorithm + due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some, + but not all, of the improvements suggested by M. Gentleman and S. + Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this + program was published in the book Software Manual for the + Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall, + Englewood Cliffs, NJ, 1980. The present program is a + translation of the Fortran 77 program in W. J. Cody, "MACHAR: + A subroutine to dynamically determine machine parameters". + TOMS (14), 1988. + + Parameter values reported are as follows: + + ibeta - the radix for the floating-point representation + it - the number of base ibeta digits in the floating-point + significand + irnd - 0 if floating-point addition chops + 1 if floating-point addition rounds, but not in the + IEEE style + 2 if floating-point addition rounds in the IEEE style + 3 if floating-point addition chops, and there is + partial underflow + 4 if floating-point addition rounds, but not in the + IEEE style, and there is partial underflow + 5 if floating-point addition rounds in the IEEE style, + and there is partial underflow + ngrd - the number of guard digits for multiplication with + truncating arithmetic. It is + 0 if floating-point arithmetic rounds, or if it + truncates and only it base ibeta digits + participate in the post-normalization shift of the + floating-point significand in multiplication; + 1 if floating-point arithmetic truncates and more + than it base ibeta digits participate in the + post-normalization shift of the floating-point + significand in multiplication. + machep - the largest negative integer such that + 1.0+FLOAT(ibeta)**machep .NE. 1.0, except that + machep is bounded below by -(it+3) + negeps - the largest negative integer such that + 1.0-FLOAT(ibeta)**negeps .NE. 1.0, except that + negeps is bounded below by -(it+3) + iexp - the number of bits (decimal places if ibeta = 10) + reserved for the representation of the exponent + (including the bias or sign) of a floating-point + number + minexp - the largest in magnitude negative integer such that + FLOAT(ibeta)**minexp is positive and normalized + maxexp - the smallest positive power of BETA that overflows + eps - the smallest positive floating-point number such + that 1.0+eps .NE. 1.0. In particular, if either + ibeta = 2 or IRND = 0, eps = FLOAT(ibeta)**machep. + Otherwise, eps = (FLOAT(ibeta)**machep)/2 + epsneg - A small positive floating-point number such that + 1.0-epsneg .NE. 1.0. In particular, if ibeta = 2 + or IRND = 0, epsneg = FLOAT(ibeta)**negeps. + Otherwise, epsneg = (ibeta**negeps)/2. Because + negeps is bounded below by -(it+3), epsneg may not + be the smallest number that can alter 1.0 by + subtraction. + xmin - the smallest non-vanishing normalized floating-point + power of the radix, i.e., xmin = FLOAT(ibeta)**minexp + xmax - the largest finite floating-point number. In + particular xmax = (1.0-epsneg)*FLOAT(ibeta)**maxexp + Note - on some machines xmax will be only the + second, or perhaps third, largest number, being + too small by 1 or 2 units in the last digit of + the significand. + + Latest revision - August 4, 1988 + + Author - W. J. Cody + Argonne National Laboratory + +*/ + +{ + int i,iz,j,k; + int mx,itmp,nxres; + REAL a,b,beta,betain,one,y,z,zero; + REAL betah,t,tmp,tmpa,tmp1,two; + + (*irnd) = 1; + one = (REAL)(*irnd); + two = one + one; + a = two; + b = a; + zero = 0.0e0; + +/* + determine ibeta,beta ala malcolm +*/ + + tmp = ((a+one)-a)-one; + + while (tmp == zero) { + a = a+a; + tmp = a+one; + tmp1 = tmp-a; + tmp = tmp1-one; + } + + tmp = a+b; + itmp = (int)(tmp-a); + while (itmp == 0) { + b = b+b; + tmp = a+b; + itmp = (int)(tmp-a); + } + + *ibeta = itmp; + beta = (REAL)(*ibeta); + +/* + determine irnd, it +*/ + + (*it) = 0; + b = one; + tmp = ((b+one)-b)-one; + + while (tmp == zero) { + *it = *it+1; + b = b*beta; + tmp = b+one; + tmp1 = tmp-b; + tmp = tmp1-one; + } + + *irnd = 0; + betah = beta/two; + tmp = a+betah; + tmp1 = tmp-a; + if (tmp1 != zero) *irnd = 1; + tmpa = a+beta; + tmp = tmpa+betah; + if ((*irnd == 0) && (tmp-tmpa != zero)) *irnd = 2; + +/* + determine negep, epsneg +*/ + + (*negep) = (*it) + 3; + betain = one / beta; + a = one; + + for (i = 1; i<=(*negep); i++) { + a = a * betain; + } + + b = a; + tmp = (one-a); + tmp = tmp-one; + + while (tmp == zero) { + a = a*beta; + *negep = *negep-1; + tmp1 = one-a; + tmp = tmp1-one; + } + + (*negep) = -(*negep); + (*epsneg) = a; + +/* + determine machep, eps +*/ + + (*machep) = -(*it) - 3; + a = b; + tmp = one+a; + + while (tmp-one == zero) { + a = a*beta; + *machep = *machep+1; + tmp = one+a; + } + + *eps = a; + +/* + determine ngrd +*/ + + (*ngrd) = 0; + tmp = one+*eps; + tmp = tmp*one; + if (((*irnd) == 0) && (tmp-one) != zero) (*ngrd) = 1; + +/* + determine iexp, minexp, xmin + + loop to determine largest i such that + (1/beta) ** (2**(i)) + does not underflow. + exit from loop is signaled by an underflow. +*/ + + i = 0; + k = 1; + z = betain; + t = one+*eps; + nxres = 0; + + for (;;) { + y = z; + z = y * y; + +/* + check for underflow +*/ + + a = z * one; + tmp = z*t; + if ((a+a == zero) || (ABS(z) > y)) break; + tmp1 = tmp*betain; + if (tmp1*beta == z) break; + i = i + 1; + k = k+k; + } + +/* + determine k such that (1/beta)**k does not underflow + first set k = 2 ** i +*/ + + (*iexp) = i + 1; + mx = k + k; + if (*ibeta == 10) { + +/* + for decimal machines only +*/ + + (*iexp) = 2; + iz = *ibeta; + while (k >= iz) { + iz = iz * (*ibeta); + (*iexp) = (*iexp) + 1; + } + mx = iz + iz - 1; + } + +/* + loop to determine minexp, xmin. + exit from loop is signaled by an underflow. +*/ + + for (;;) { + (*xmin) = y; + y = y * betain; + a = y * one; + tmp = y*t; + tmp1 = a+a; + if ((tmp1 == zero) || (ABS(y) >= (*xmin))) break; + k = k + 1; + tmp1 = tmp*betain; + tmp1 = tmp1*beta; + + if ((tmp1 == y) && (tmp != y)) { + nxres = 3; + *xmin = y; + break; + } + + } + + (*minexp) = -k; + +/* + determine maxexp, xmax +*/ + + if ((mx <= k+k-3) && ((*ibeta) != 10)) { + mx = mx + mx; + (*iexp) = (*iexp) + 1; + } + + (*maxexp) = mx + (*minexp); + +/* + Adjust *irnd to reflect partial underflow. +*/ + + (*irnd) = (*irnd)+nxres; + +/* + Adjust for IEEE style machines. +*/ + + if ((*irnd) >= 2) (*maxexp) = (*maxexp)-2; + +/* + adjust for machines with implicit leading bit in binary + significand and machines with radix point at extreme + right of significand. +*/ + + i = (*maxexp) + (*minexp); + if (((*ibeta) == 2) && (i == 0)) (*maxexp) = (*maxexp) - 1; + if (i > 20) (*maxexp) = (*maxexp) - 1; + if (a != y) (*maxexp) = (*maxexp) - 2; + (*xmax) = one - (*epsneg); + tmp = (*xmax)*one; + if (tmp != (*xmax)) (*xmax) = one - beta * (*epsneg); + (*xmax) = (*xmax) / (beta * beta * beta * (*xmin)); + i = (*maxexp) + (*minexp) + 3; + if (i > 0) { + + for (j = 1; j<=i; j++ ) { + if ((*ibeta) == 2) (*xmax) = (*xmax) + (*xmax); + if ((*ibeta) != 2) (*xmax) = (*xmax) * beta; + } + + } + + return; + +} + +void +#if defined (F77_APPEND_UNDERSCORE) +machar_ (REAL *xmin, REAL *xmax, REAL *epsneg, REAL *eps, REAL *log10_ibeta) +#else +machar (REAL *xmin, REAL *xmax, REAL *epsneg, REAL *eps, REAL *log10_ibeta) +#endif +{ + int ibeta, iexp, irnd, it, machep, maxexp, minexp, negep, ngrd; + + rmachar (&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp, + &maxexp, eps, epsneg, xmin, xmax); + + *log10_ibeta = log10 ((REAL) ibeta); +}