The R Project SVN R

/* Formerly src/appl/machar.c:
 * void machar()  -- computes ALL `machine constants' at once.
 * -------------  -- compare with ../nmath/i1mach.c & ../nmath/d1mach.c
 *           which use the C  <float.h> constants !
 *      algorithm 665, collected algorithms from acm.
 *      this work published in transactions on mathematical software,
 *      vol. 14, no. 4, pp. 303-311.
 *
 *  this fortran 77 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 program as given here must be modified before compiling.  if
 *   a single (double) precision version is desired, change all
 *   occurrences of cs (  ) in columns 1 and 2 to blanks.
 *
 *  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 - april 20, 1987
 *
 *     author - w. j. cody
 *              argonne national laboratory
 *
 */

/* #include  from ./platform.c , now as
 a "template" to be used with    DTYPE in {double, long double, ..}
                             MACH_NAME in {machar, machar_LD,   ..}
                     ABS() in {fabs(), fabsl()  ,   ..}
*/
static void
MACH_NAME(int *ibeta, int *it, int *irnd, int *ngrd, int *machep, int *negep,
      int *iexp, int *minexp, int *maxexp,
      DTYPE *eps, DTYPE *epsneg, DTYPE *xmin, DTYPE *xmax)
{
    volatile DTYPE a, b, beta, betain, betah, one,
        t, temp, tempa, temp1, two, y, z, zero;
    int i, iz, j, k, mx, nxres;

    one = 1;
    two = one+one;
    zero = one-one;

        /* determine ibeta, beta ala malcolm. */
    a = one; // a = <large> = 9.0072e+15 for 'double' is used later
    do {
        a = a + a;
        temp = a + one;
        temp1 = temp - a;
    }
    while(temp1 - one == zero);
#ifdef _no_longer___did_overflow_ // on IBM PowerPPC ('Power 8')
    int itemp;
    b = one;
    do {
        b = b + b;
        temp = a + b;
        itemp = (int)(temp - a);
    }
    while (itemp == 0);
    *ibeta = itemp;
#else
    *ibeta = (int) FLT_RADIX;
#endif
    beta = *ibeta;

        /* determine it, irnd */

    *it = 0;
    b = one;
    do {
        *it = *it + 1;
        b = b * beta;
        temp = b + one;
        temp1 = temp - b;
    }
    while(temp1 - one == zero);
    *irnd = 0;
    betah = beta / two;
    temp = a + betah;
    if (temp - a != zero)
        *irnd = 1;
    tempa = a + beta;
    temp = tempa + betah;
    if (*irnd == 0 && temp - tempa != 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;
    for(;;) {
        temp = one - a;
        if (temp - one != zero)
            break;
        a = a * beta;
        *negep = *negep - 1;
    }
    *negep = -*negep;
    *epsneg = a;
    if (*ibeta != 2 && *irnd != 0) {
        a = (a * (one + a)) / two;
        temp = one - a;
        if (temp - one != zero)
            *epsneg = a;
    }

        /* determine machep, eps */

    *machep = -*it - 3;
    a = b;
    for(;;) {
        temp = one + a;
        if (temp - one != zero)
            break;
        a = a * beta;
        *machep = *machep + 1;
    }
    *eps = a;
    temp = tempa + beta * (one + *eps);
    if (*ibeta != 2 && *irnd != 0) {
        a = (a * (one + a)) / two;
        temp = one + a;
        if (temp - one != zero)
            *eps = a;
    }

        /* determine ngrd */

    *ngrd = 0;
    temp = one + *eps;
    if (*irnd == 0 && temp * one - one != zero)
        *ngrd = 1;

    /* determine iexp, minexp, xmin */

    /* loop to determine largest i and k = 2**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 here */

        a = z * one;
        temp = z * t;
        if (a+a == zero || ABS(z) >= y)
            break;
        temp1 = temp * betain;
        if (temp1 * beta == z)
            break;
        i = i+1;
        k = k+k;
    }
    if (*ibeta != 10) {
        *iexp = i + 1;
        mx = k + k;
    }
    else {
        /* this segment is for decimal machines only */

        *iexp = 2;
        iz = *ibeta;
        while (k >= iz) {
            iz = iz * *ibeta;
            iexp = iexp + 1;
        }
        mx = iz + iz - 1;
    }
    do {
        /* loop to determine minexp, xmin */
        /* exit from loop is signaled by an underflow */

        *xmin = y;
        y = y * betain;

        /* check for underflow here */

        a = y * one;
        temp = y * t;
        if (a+a == zero || ABS(y) >= *xmin)
            goto L10;
        k = k + 1;
        temp1 = temp * betain;
    }
    while(temp1 * beta != y);
    nxres = 3;
    *xmin = y;
L10:    *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 || *irnd == 5)
        *maxexp = *maxexp - 2;

    /* adjust for non-ieee machines with partial underflow */

    if (*irnd == 3 || *irnd == 4)
        *maxexp = *maxexp - *it;

    /* 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;
    if (*xmax * one != *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;
        }
}