From GOLUB@SU-SCORE.ARPA Wed Oct 16 21:27:10 1985 Received: from SU-SCORE.ARPA (su-score.arpa.ARPA) by anl-mcs.ARPA ; Wed, 16 Oct 85 21:26:37 cdt Return-Path: Received: from CSNET-RELAY.ARPA by SU-SCORE.ARPA with TCP; Wed 16 Oct 85 13:52:49-PDT Received: from waterloo by csnet-relay.csnet id ad25408; 16 Oct 85 16:48 EDT Date: Wed, 16 Oct 85 12:38:29 edt From: Richard Bartels Message-Id: <8510161638.AA01992@watcgl.UUCP> To: na.dis@su-score.ARPA Subject: UNIX/f77 constants I Resent-Date: Wed 16 Oct 85 18:48:23-PDT Resent-From: Gene Golub (415/497-3124) Resent-To: :; Resent-Message-Id: <12151704851.14.GOLUB@SU-SCORE.ARPA> Status: RO The following, gratuitously submitted, is the SMCHAR routine from MINPACK hacked into C. To change to double, change all ocurrences of "float" in the code into "double". This program works equivalently to the FORTRAN version distributed with MINPACK and run under UNIX/f77 except that the octal numbers are really printed out in octal. Suggestions, bug fixes, and improvements solicited. ------------------- cut here ----------------------------------------------- #include #include #define ZERO 0.0 #define ABS(xxx) ((xxx>ZERO)?(xxx):(-xxx)) main() { /* c ********** c c this program prints hardware-determined c machine constants obtained by smchar, a subroutine due to c w. j. cody. c c descriptions of the machine constants are c given in the prologue comments of smchar. c c subprograms called c c smchar c c minpack. version of january 1979. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c modified: october 16, 1985 c richard bartels c c version (sort of): october 16 1985 c richard bartels c c ********** */ int ibeta,iexp,irnd,it,machep,maxexp,minexp,negep,ngrd; float dwarf,eps,epsmch,epsneg,giant,xmax,xmin; smchar(&ibeta,&it,&irnd,&ngrd,&machep,&negep, &iexp,&minexp,&maxexp,&eps,&epsneg,&xmin,&xmax); printf("smchar constants\n"); printf("ibeta = %d\n",ibeta); printf("it = %d\n",it); printf("irnd = %d\n",irnd); printf("ngrd = %d\n",ngrd); printf("machep = %d\n",machep); printf("negep = %d\n",negep); printf("iexp = %d\n",iexp); printf("minexp = %d\n",minexp); printf("maxexp = %d\n",maxexp); printf("eps = %o\n",eps); printf("epsneg = %o\n",epsneg); printf("xmin = %o\n",xmin); printf("xmax = %o\n",xmax); } smchar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp, maxexp,eps,epsneg,xmin,xmax) int *ibeta,*iexp,*irnd,*it,*machep,*maxexp,*minexp,*negep,*ngrd; float *eps,*epsneg,*xmax,*xmin; { int i,iz,j,k; int mx; float a,b,beta,betain,betam1,one,y,z,zero; /* c c this subroutine is intended to determine the characteristics c of the floating-point arithmetic system that are specified c below. the first three are determined according to an c algorithm due to m. malcolm, cacm 15 (1972), pp. 949-951, c incorporating some, but not all, of the improvements c suggested by m. gentleman and s. marovich, cacm 17 (1974), c pp. 276-277. the version given here is for single precision. c c latest revision - june 25, 1979 c c author - w. j. cody c argonne national laboratory c c c version (sort of) - october 16, 1985 c c hacker - richard bartels c university of waterloo c c ibeta is the radix of the floating-point representation c it is the number of base ibeta digits in the floating-point c significand c irnd = 0 if the arithmetic chops, c 1 if the arithmetic rounds c ngrd = 0 if irnd=1, or if irnd=0 and only it base ibeta c digits participate in the post normalization shift c of the floating-point significand in multiplication c 1 if irnd=0 and more than it base ibeta digits c participate in the post normalization shift of the c floating-point significand in multiplication c machep is the largest negative integer such that c 1.0+float(ibeta)**machep .ne. 1.0, except that c machep is bounded below by it-3 c negeps is the largest negative integer such that c 1.0-float(ibeta)**negeps .ne. 1.0, except that c negeps is bounded below by it-3 c iexp is the number of bits (decimal places if ibeta = 10) c reserved for the representation of the exponent c (including the bias or sign) of a floating-point c number c minexp is the largest in magnitude negative integer such that c float(ibeta)**minexp is a positive floating-point c number c maxexp is the largest positive integer exponent for a finite c floating-point number c eps is the smallest positive floating-point number such c that 1.0+eps .ne. 1.0. in particular, if either c ibeta = 2 or irnd = 0, eps = float(ibeta)**machep. c otherwise, eps = (float(ibeta)**machep)/2 c epsneg is a small positive floating-point number such that c 1.0-epsneg .ne. 1.0. in particular, if ibeta = 2 c or irnd = 0, epsneg = float(ibeta)**negeps. c otherwise, epsneg = (ibeta**negeps)/2. because c negeps is bounded below by it-3, epsneg may not c be the smallest number which can alter 1.0 by c subtraction. c xmin is the smallest positive floating-point power of the c radix. in particular, xmin = float(ibeta)**minexp c xmax is the largest finite floating-point number. in c particular xmax = (1.0-epsneg)*float(ibeta)**maxexp c note - on some machines xmax will be only the c second, or perhaps third, largest number, being c too small by 1 or 2 units in the last digit of c the significand. c c----------------------------------------------------------------- */ (*irnd) = 1; one = (float)(*irnd); a = one + one; b = a; zero = 0.0e0; /* c----------------------------------------------------------------- c determine ibeta,beta ala malcolm c----------------------------------------------------------------- */ s10: if (((a+one)-a)-one != zero) goto s20; a = a + a; goto s10; s20: if ((a+b)-a != zero) goto s30; b = b + b; goto s20; s30: *ibeta = (int)((a+b)-a); beta = (float)(*ibeta); betam1 = beta - one; /* c----------------------------------------------------------------- c determine irnd, it c----------------------------------------------------------------- */ if ((a+betam1)-a == zero) (*irnd) = 0; (*it) = 0; a = one; s100: (*it) = (*it) + 1; a = a * beta; if (((a+one)-a)-one == zero) goto s100; /* c----------------------------------------------------------------- c determine negep, epsneg c----------------------------------------------------------------- */ (*negep) = (*it) + 3; a = one; for (i = 1; i<=(*negep); i++) { a = a / beta; } b = a; s210: if ((one-a)-one != zero) goto s220; a = a * beta; (*negep) = (*negep) - 1; goto s210; s220: (*negep) = -(*negep); (*epsneg) = a; if ((*ibeta == 2) || ((*irnd) == 0)) goto s300; a = (a*(one+a)) / (one+one); if ((one-a)-one != zero) (*epsneg) = a; /* c----------------------------------------------------------------- c determine machep, eps c----------------------------------------------------------------- */ s300: (*machep) = -(*it) - 3; a = b; s310: if((one+a)-one != zero) goto s320; a = a * beta; (*machep) = (*machep) + 1; goto s310; s320: (*eps) = a; if ((*ibeta == 2) || ((*irnd) == 0)) goto s350; a = (a*(one+a)) / (one+one); if ((one+a)-one != zero) (*eps) = a; /* c----------------------------------------------------------------- c determine ngrd c----------------------------------------------------------------- */ s350: (*ngrd) = 0; if (((*irnd) == 0) && ((one+(*eps))*one-one) != zero) (*ngrd) = 1; /* c----------------------------------------------------------------- c determine iexp, minexp, xmin c c loop to determine largest i such that c (1/beta) ** (2**(i)) c does not underflow c exit from loop is signaled by an underflow. c----------------------------------------------------------------- */ i = 0; betain = one / beta; z = betain; s400: y = z; z = y * y; /* c----------------------------------------------------------------- c check for underflow c----------------------------------------------------------------- */ a = z * one; if ((a+a == zero) || (ABS(z) > y)) goto s410; i = i + 1; goto s400; /* c----------------------------------------------------------------- c determine k such that (1/beta)**k does not underflow c c first set k = 2 ** i c----------------------------------------------------------------- */ s410: k = 1; for (j = 1; j<=i; j++ ) { k = k + k; } (*iexp) = i + 1; mx = k + k; if (*ibeta != 10) goto s450; /* c----------------------------------------------------------------- c for decimal machines only c----------------------------------------------------------------- */ (*iexp) = 2; iz = *ibeta; s430: if (k < iz) goto s440; iz = iz * (*ibeta); (*iexp) = (*iexp) + 1; goto s430; s440: mx = iz + iz - 1; /* c----------------------------------------------------------------- c loop to construct xmin. c exit from loop is signaled by an underflow. c----------------------------------------------------------------- */ s450: (*xmin) = y; y = y * betain; a = y * one; if (((a+a) == zero) || (ABS(y) > (*xmin))) goto s460; k = k + 1; goto s450; s460: (*minexp) = -k; /* c----------------------------------------------------------------- c determine maxexp, xmax c----------------------------------------------------------------- */ if ((mx > k+k-3) || ((*ibeta) == 10)) goto s500; mx = mx + mx; (*iexp) = (*iexp) + 1; s500: (*maxexp) = mx + (*minexp); /* c----------------------------------------------------------------- c adjust for machines with implicit leading c bit in binary significand and machines with c radix point at extreme right of significand c----------------------------------------------------------------- */ 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); if ((*xmax)*one != (*xmax)) (*xmax) = one - beta * (*epsneg); (*xmax) = (*xmax) / (beta * beta * beta * (*xmin)); i = (*maxexp) + (*minexp) + 3; if (i <= 0) goto s520; for (j = 1; j<=i; j++ ) { if ((*ibeta) == 2) (*xmax) = (*xmax) + (*xmax); if ((*ibeta) != 2) (*xmax) = (*xmax) * beta; } s520: return; } From GOLUB@SU-SCORE.ARPA Wed Oct 16 22:37:48 1985 Received: from SU-SCORE.ARPA (su-score.arpa.ARPA) by anl-mcs.ARPA ; Wed, 16 Oct 85 22:37:11 cdt Return-Path: Received: from CSNET-RELAY.ARPA by SU-SCORE.ARPA with TCP; Wed 16 Oct 85 13:52:49-PDT Received: from waterloo by csnet-relay.csnet id ad25408; 16 Oct 85 16:48 EDT Date: Wed, 16 Oct 85 12:38:29 edt From: Richard Bartels Message-Id: <8510161638.AA01992@watcgl.UUCP> To: na.dis@su-score.ARPA Subject: UNIX/f77 constants I Resent-Date: Wed 16 Oct 85 19:01:21-PDT Resent-From: Gene Golub (415/497-3124) Resent-To: :; Resent-Message-Id: <12151707211.14.GOLUB@SU-SCORE.ARPA> Status: RO The following, gratuitously submitted, is the SMCHAR routine from MINPACK hacked into C. To change to double, change all ocurrences of "float" in the code into "double". This program works equivalently to the FORTRAN version distributed with MINPACK and run under UNIX/f77 except that the octal numbers are really printed out in octal. Suggestions, bug fixes, and improvements solicited. ------------------- cut here ----------------------------------------------- #include #include #define ZERO 0.0 #define ABS(xxx) ((xxx>ZERO)?(xxx):(-xxx)) main() { /* c ********** c c this program prints hardware-determined c machine constants obtained by smchar, a subroutine due to c w. j. cody. c c descriptions of the machine constants are c given in the prologue comments of smchar. c c subprograms called c c smchar c c minpack. version of january 1979. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c modified: october 16, 1985 c richard bartels c c version (sort of): october 16 1985 c richard bartels c c ********** */ int ibeta,iexp,irnd,it,machep,maxexp,minexp,negep,ngrd; float dwarf,eps,epsmch,epsneg,giant,xmax,xmin; smchar(&ibeta,&it,&irnd,&ngrd,&machep,&negep, &iexp,&minexp,&maxexp,&eps,&epsneg,&xmin,&xmax); printf("smchar constants\n"); printf("ibeta = %d\n",ibeta); printf("it = %d\n",it); printf("irnd = %d\n",irnd); printf("ngrd = %d\n",ngrd); printf("machep = %d\n",machep); printf("negep = %d\n",negep); printf("iexp = %d\n",iexp); printf("minexp = %d\n",minexp); printf("maxexp = %d\n",maxexp); printf("eps = %o\n",eps); printf("epsneg = %o\n",epsneg); printf("xmin = %o\n",xmin); printf("xmax = %o\n",xmax); } smchar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp, maxexp,eps,epsneg,xmin,xmax) int *ibeta,*iexp,*irnd,*it,*machep,*maxexp,*minexp,*negep,*ngrd; float *eps,*epsneg,*xmax,*xmin; { int i,iz,j,k; int mx; float a,b,beta,betain,betam1,one,y,z,zero; /* c c this subroutine is intended to determine the characteristics c of the floating-point arithmetic system that are specified c below. the first three are determined according to an c algorithm due to m. malcolm, cacm 15 (1972), pp. 949-951, c incorporating some, but not all, of the improvements c suggested by m. gentleman and s. marovich, cacm 17 (1974), c pp. 276-277. the version given here is for single precision. c c latest revision - june 25, 1979 c c author - w. j. cody c argonne national laboratory c c c version (sort of) - october 16, 1985 c c hacker - richard bartels c university of waterloo c c ibeta is the radix of the floating-point representation c it is the number of base ibeta digits in the floating-point c significand c irnd = 0 if the arithmetic chops, c 1 if the arithmetic rounds c ngrd = 0 if irnd=1, or if irnd=0 and only it base ibeta c digits participate in the post normalization shift c of the floating-point significand in multiplication c 1 if irnd=0 and more than it base ibeta digits c participate in the post normalization shift of the c floating-point significand in multiplication c machep is the largest negative integer such that c 1.0+float(ibeta)**machep .ne. 1.0, except that c machep is bounded below by it-3 c negeps is the largest negative integer such that c 1.0-float(ibeta)**negeps .ne. 1.0, except that c negeps is bounded below by it-3 c iexp is the number of bits (decimal places if ibeta = 10) c reserved for the representation of the exponent c (including the bias or sign) of a floating-point c number c minexp is the largest in magnitude negative integer such that c float(ibeta)**minexp is a positive floating-point c number c maxexp is the largest positive integer exponent for a finite c floating-point number c eps is the smallest positive floating-point number such c that 1.0+eps .ne. 1.0. in particular, if either c ibeta = 2 or irnd = 0, eps = float(ibeta)**machep. c otherwise, eps = (float(ibeta)**machep)/2 c epsneg is a small positive floating-point number such that c 1.0-epsneg .ne. 1.0. in particular, if ibeta = 2 c or irnd = 0, epsneg = float(ibeta)**negeps. c otherwise, epsneg = (ibeta**negeps)/2. because c negeps is bounded below by it-3, epsneg may not c be the smallest number which can alter 1.0 by c subtraction. c xmin is the smallest positive floating-point power of the c radix. in particular, xmin = float(ibeta)**minexp c xmax is the largest finite floating-point number. in c particular xmax = (1.0-epsneg)*float(ibeta)**maxexp c note - on some machines xmax will be only the c second, or perhaps third, largest number, being c too small by 1 or 2 units in the last digit of c the significand. c c----------------------------------------------------------------- */ (*irnd) = 1; one = (float)(*irnd); a = one + one; b = a; zero = 0.0e0; /* c----------------------------------------------------------------- c determine ibeta,beta ala malcolm c----------------------------------------------------------------- */ s10: if (((a+one)-a)-one != zero) goto s20; a = a + a; goto s10; s20: if ((a+b)-a != zero) goto s30; b = b + b; goto s20; s30: *ibeta = (int)((a+b)-a); beta = (float)(*ibeta); betam1 = beta - one; /* c----------------------------------------------------------------- c determine irnd, it c----------------------------------------------------------------- */ if ((a+betam1)-a == zero) (*irnd) = 0; (*it) = 0; a = one; s100: (*it) = (*it) + 1; a = a * beta; if (((a+one)-a)-one == zero) goto s100; /* c----------------------------------------------------------------- c determine negep, epsneg c----------------------------------------------------------------- */ (*negep) = (*it) + 3; a = one; for (i = 1; i<=(*negep); i++) { a = a / beta; } b = a; s210: if ((one-a)-one != zero) goto s220; a = a * beta; (*negep) = (*negep) - 1; goto s210; s220: (*negep) = -(*negep); (*epsneg) = a; if ((*ibeta == 2) || ((*irnd) == 0)) goto s300; a = (a*(one+a)) / (one+one); if ((one-a)-one != zero) (*epsneg) = a; /* c----------------------------------------------------------------- c determine machep, eps c----------------------------------------------------------------- */ s300: (*machep) = -(*it) - 3; a = b; s310: if((one+a)-one != zero) goto s320; a = a * beta; (*machep) = (*machep) + 1; goto s310; s320: (*eps) = a; if ((*ibeta == 2) || ((*irnd) == 0)) goto s350; a = (a*(one+a)) / (one+one); if ((one+a)-one != zero) (*eps) = a; /* c----------------------------------------------------------------- c determine ngrd c----------------------------------------------------------------- */ s350: (*ngrd) = 0; if (((*irnd) == 0) && ((one+(*eps))*one-one) != zero) (*ngrd) = 1; /* c----------------------------------------------------------------- c determine iexp, minexp, xmin c c loop to determine largest i such that c (1/beta) ** (2**(i)) c does not underflow c exit from loop is signaled by an underflow. c----------------------------------------------------------------- */ i = 0; betain = one / beta; z = betain; s400: y = z; z = y * y; /* c----------------------------------------------------------------- c check for underflow c----------------------------------------------------------------- */ a = z * one; if ((a+a == zero) || (ABS(z) > y)) goto s410; i = i + 1; goto s400; /* c----------------------------------------------------------------- c determine k such that (1/beta)**k does not underflow c c first set k = 2 ** i c----------------------------------------------------------------- */ s410: k = 1; for (j = 1; j<=i; j++ ) { k = k + k; } (*iexp) = i + 1; mx = k + k; if (*ibeta != 10) goto s450; /* c----------------------------------------------------------------- c for decimal machines only c----------------------------------------------------------------- */ (*iexp) = 2; iz = *ibeta; s430: if (k < iz) goto s440; iz = iz * (*ibeta); (*iexp) = (*iexp) + 1; goto s430; s440: mx = iz + iz - 1; /* c----------------------------------------------------------------- c loop to construct xmin. c exit from loop is signaled by an underflow. c----------------------------------------------------------------- */ s450: (*xmin) = y; y = y * betain; a = y * one; if (((a+a) == zero) || (ABS(y) > (*xmin))) goto s460; k = k + 1; goto s450; s460: (*minexp) = -k; /* c----------------------------------------------------------------- c determine maxexp, xmax c----------------------------------------------------------------- */ if ((mx > k+k-3) || ((*ibeta) == 10)) goto s500; mx = mx + mx; (*iexp) = (*iexp) + 1; s500: (*maxexp) = mx + (*minexp); /* c----------------------------------------------------------------- c adjust for machines with implicit leading c bit in binary significand and machines with c radix point at extreme right of significand c----------------------------------------------------------------- */ 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); if ((*xmax)*one != (*xmax)) (*xmax) = one - beta * (*epsneg); (*xmax) = (*xmax) / (beta * beta * beta * (*xmin)); i = (*maxexp) + (*minexp) + 3; if (i <= 0) goto s520; for (j = 1; j<=i; j++ ) { if ((*ibeta) == 2) (*xmax) = (*xmax) + (*xmax); if ((*ibeta) != 2) (*xmax) = (*xmax) * beta; } s520: return; } .