Re: R-alpha: Machine Constants

George White (gwhite@cabot.bio.dfo.ca)
Fri, 6 Dec 96 8:06:53 AST


Message-Id: <199612061203.NAA24909@hypatia.math.ethz.ch>
From: George White <gwhite@cabot.bio.dfo.ca>
Subject: Re: R-alpha: Machine Constants
To: r-testers@stat.math.ethz.ch
Date: Fri, 6 Dec 96 8:06:53 AST

My two bits to the discussion:

First comment: I don't recall seeing anything for the IEEE rounding
mode in the standard .h files.  This is the setting that has caused
me the greatest problems with portability in the past, so I would like
to see it included in .Machine, and since it can change for the same
binary, it needs to determined at runtime.

2nd comment:

In Ross's machar.c, which appears to have been converted from the 
Fortran, it says:

 *     latest revision - april 20, 1987

I just checked the current machar.c on netlib, and I see:
      Latest revision - August 4, 1988

Perhaps someone who is seeing incorrect results will try the netlib 
version with -DDP added to the compiler flags (included below), and
with the addition of "volatile" if it appears necessary.  On 
the SGI Irix 5.3 system here I get (in agreement with .Machine):

$ uname -a
IRIX caligo 5.3 11091812 IP22 mips
$ hinv
Iris Audio Processor: version A2 revision 4.1.0
1 132 MHZ IP22 Processor
FPU: MIPS R4610 Floating Point Chip Revision: 2.0
CPU: MIPS R4600 Processor Chip Revision: 2.0
[...]
$ cc -DDP -mips2 -O machar.c && ./a.out
Double  precision MACHAR constants
ibeta  = 2
it     = 53
irnd   = 5
ngrd   = 0
machep = -52
negep  = -53
iexp   = 11
minexp = -1022
maxexp = 1024
eps      2.2204460492503131e-16   3CB00000          0
epsneg   1.1102230246251565e-16   3CA00000          0
xmin    2.2250738585072014e-308     100000          0
xmax    1.7976931348623157e+308   7FEFFFFF   FFFFFFFF

The source for the newer version is:

/* obtained from netlib, 5 Dec., 1996 */
/*

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))


main()

{

/*
      This program prints hardware-determined double-precision
      machine constants obtained from rmachar.  Dmachar is a C
      translation of the Fortran routine MACHAR from W. J. Cody,
      "MACHAR: A subroutine to dynamically determine machine
      parameters".  TOMS (14), 1988.

      Descriptions of the machine constants are given in the
      prologue comments in rmachar.

      Subprograms called

        rmachar

      Original driver: Richard Bartels, October 16, 1985

      Modified by: W. J. Cody
                   July 26, 1988

      **********
*/
      int ibeta,iexp,irnd,it,machep,maxexp,minexp,negep,ngrd;
        int i ;
      REAL eps,epsneg,xmax,xmin;
      union wjc{
              long int jj[REALSIZE];
              REAL xbig;
         } uval;

      rmachar(&ibeta,&it,&irnd,&ngrd,&machep,&negep,
              &iexp,&minexp,&maxexp,&eps,&epsneg,&xmin,&xmax);
      printf(PREC);
        printf(" precision MACHAR 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);

#define DISPLAY(s,x) { \
        uval.xbig = x ; \
        printf(s); \
        printf(" %24.16e ",(double) x) ; \
        for(i=0;i<REALSIZE;i++) printf(" %9X ",uval.jj[i]) ; \
        printf("\n"); \
        }
                        
      DISPLAY("eps   ",eps);
      DISPLAY("epsneg",epsneg);
      DISPLAY("xmin  ",xmin);
      DISPLAY("xmax  ",xmax);

}


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;

}

--
  George N. White III                | Internet: <aa056@chebucto.ns.ca>
  Department of Fisheries and Oceans | DFO email: GWhite@BIOnet.BIO.DFO.ca
  Ocean Sciences Division            | Time Zone: Atlantic (AST4ADT)
  Bedford Institute of Oceanography  | phone: (08:00-16:00)
  P.O. Box 1006                      |     902.426.8509 (office)
  Dartmouth, Nova Scotia             |     902.426.9388 (FAX)
  CANADA               B2Y 4A2       |     902.426.3793 (emergencies)
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
r-testers mailing list -- For info or help, send "info" or "help",
To [un]subscribe, send "[un]subscribe"
(in the "body", not the subject !)  To: r-testers-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-