home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / misc / gen-d1mach.c < prev    next >
C/C++ Source or Header  |  1996-09-28  |  11KB  |  416 lines
  1. /*
  2.  
  3. This file combines the single and double precision versions of machar,
  4. selected by cc -DSP or cc -DDP.  This feature provided by D. G. Hough,
  5. August 3, 1988.
  6.  
  7. */
  8.  
  9. #ifdef SP
  10. #define REAL float
  11. #define ZERO 0.0
  12. #define ONE 1.0
  13. #define PREC "Single "
  14. #define REALSIZE 1
  15. #endif
  16.  
  17. #ifdef DP
  18. #define REAL double
  19. #define ZERO 0.0e0
  20. #define ONE 1.0e0
  21. #define PREC "Double "
  22. #define REALSIZE 2
  23. #endif
  24.  
  25. #include <math.h>
  26. #include <stdio.h>
  27.  
  28. #define ABS(xxx) ((xxx>ZERO)?(xxx):(-xxx))
  29.  
  30. void
  31. rmachar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,
  32.         maxexp,eps,epsneg,xmin,xmax)
  33.  
  34.       int *ibeta,*iexp,*irnd,*it,*machep,*maxexp,*minexp,*negep,*ngrd;
  35.       REAL *eps,*epsneg,*xmax,*xmin;
  36.  
  37. /*
  38.  
  39.    This subroutine is intended to determine the parameters of the
  40.     floating-point arithmetic system specified below.  The
  41.     determination of the first three uses an extension of an algorithm
  42.     due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
  43.     but not all, of the improvements suggested by M. Gentleman and S.
  44.     Marovich, CACM 17 (1974), pp. 276-277.  An earlier version of this
  45.     program was published in the book Software Manual for the
  46.     Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
  47.     Englewood Cliffs, NJ, 1980.  The present program is a
  48.     translation of the Fortran 77 program in W. J. Cody, "MACHAR:
  49.     A subroutine to dynamically determine machine parameters".
  50.     TOMS (14), 1988.
  51.  
  52.    Parameter values reported are as follows:
  53.  
  54.         ibeta   - the radix for the floating-point representation
  55.         it      - the number of base ibeta digits in the floating-point
  56.                   significand
  57.         irnd    - 0 if floating-point addition chops
  58.                   1 if floating-point addition rounds, but not in the
  59.                     IEEE style
  60.                   2 if floating-point addition rounds in the IEEE style
  61.                   3 if floating-point addition chops, and there is
  62.                     partial underflow
  63.                   4 if floating-point addition rounds, but not in the
  64.                     IEEE style, and there is partial underflow
  65.                   5 if floating-point addition rounds in the IEEE style,
  66.                     and there is partial underflow
  67.         ngrd    - the number of guard digits for multiplication with
  68.                   truncating arithmetic.  It is
  69.                   0 if floating-point arithmetic rounds, or if it
  70.                     truncates and only  it  base  ibeta digits
  71.                     participate in the post-normalization shift of the
  72.                     floating-point significand in multiplication;
  73.                   1 if floating-point arithmetic truncates and more
  74.                     than  it  base  ibeta  digits participate in the
  75.                     post-normalization shift of the floating-point
  76.                     significand in multiplication.
  77.         machep  - the largest negative integer such that
  78.                   1.0+FLOAT(ibeta)**machep .NE. 1.0, except that
  79.                   machep is bounded below by  -(it+3)
  80.         negeps  - the largest negative integer such that
  81.                   1.0-FLOAT(ibeta)**negeps .NE. 1.0, except that
  82.                   negeps is bounded below by  -(it+3)
  83.         iexp    - the number of bits (decimal places if ibeta = 10)
  84.                   reserved for the representation of the exponent
  85.                   (including the bias or sign) of a floating-point
  86.                   number
  87.         minexp  - the largest in magnitude negative integer such that
  88.                   FLOAT(ibeta)**minexp is positive and normalized
  89.         maxexp  - the smallest positive power of  BETA  that overflows
  90.         eps     - the smallest positive floating-point number such
  91.                   that  1.0+eps .NE. 1.0. In particular, if either
  92.                   ibeta = 2  or  IRND = 0, eps = FLOAT(ibeta)**machep.
  93.                   Otherwise,  eps = (FLOAT(ibeta)**machep)/2
  94.         epsneg  - A small positive floating-point number such that
  95.                   1.0-epsneg .NE. 1.0. In particular, if ibeta = 2
  96.                   or  IRND = 0, epsneg = FLOAT(ibeta)**negeps.
  97.                   Otherwise,  epsneg = (ibeta**negeps)/2.  Because
  98.                   negeps is bounded below by -(it+3), epsneg may not
  99.                   be the smallest number that can alter 1.0 by
  100.                   subtraction.
  101.         xmin    - the smallest non-vanishing normalized floating-point
  102.                   power of the radix, i.e.,  xmin = FLOAT(ibeta)**minexp
  103.         xmax    - the largest finite floating-point number.  In
  104.                   particular  xmax = (1.0-epsneg)*FLOAT(ibeta)**maxexp
  105.                   Note - on some machines  xmax  will be only the
  106.                   second, or perhaps third, largest number, being
  107.                   too small by 1 or 2 units in the last digit of
  108.                   the significand.
  109.  
  110.       Latest revision - August 4, 1988
  111.  
  112.       Author - W. J. Cody
  113.                Argonne National Laboratory
  114.  
  115. */
  116.  
  117. {
  118.       int i,iz,j,k;
  119.       int mx,itmp,nxres;
  120.       REAL a,b,beta,betain,one,y,z,zero;
  121.       REAL betah,t,tmp,tmpa,tmp1,two;
  122.  
  123.       (*irnd) = 1;
  124.       one = (REAL)(*irnd);
  125.       two = one + one;
  126.       a = two;
  127.       b = a;
  128.       zero = 0.0e0;
  129.  
  130. /*
  131.   determine ibeta,beta ala malcolm
  132. */
  133.  
  134.       tmp = ((a+one)-a)-one;
  135.  
  136.       while (tmp == zero) {
  137.          a = a+a;
  138.          tmp = a+one;
  139.          tmp1 = tmp-a;
  140.          tmp = tmp1-one;
  141.       }
  142.  
  143.       tmp = a+b;
  144.       itmp = (int)(tmp-a);
  145.       while (itmp == 0) {
  146.          b = b+b;
  147.          tmp = a+b;
  148.          itmp = (int)(tmp-a);
  149.       }
  150.  
  151.       *ibeta = itmp;
  152.       beta = (REAL)(*ibeta);
  153.  
  154. /*
  155.   determine irnd, it
  156. */
  157.  
  158.       (*it) = 0;
  159.       b = one;
  160.       tmp = ((b+one)-b)-one;
  161.  
  162.       while (tmp == zero) {
  163.          *it = *it+1;
  164.          b = b*beta;
  165.          tmp = b+one;
  166.          tmp1 = tmp-b;
  167.          tmp = tmp1-one;
  168.       }
  169.  
  170.       *irnd = 0;
  171.       betah = beta/two;
  172.       tmp = a+betah;
  173.       tmp1 = tmp-a;
  174.       if (tmp1 != zero) *irnd = 1;
  175.       tmpa = a+beta;
  176.       tmp = tmpa+betah;
  177.       if ((*irnd == 0) && (tmp-tmpa != zero)) *irnd = 2;
  178.  
  179. /*
  180.   determine negep, epsneg
  181. */
  182.  
  183.       (*negep) = (*it) + 3;
  184.       betain = one / beta;
  185.       a = one;
  186.  
  187.       for (i = 1; i<=(*negep); i++) {
  188.          a = a * betain;
  189.       }
  190.  
  191.       b = a;
  192.       tmp = (one-a);
  193.       tmp = tmp-one;
  194.  
  195.       while (tmp == zero) {
  196.          a = a*beta;
  197.          *negep = *negep-1;
  198.          tmp1 = one-a;
  199.          tmp = tmp1-one;
  200.       }
  201.  
  202.       (*negep) = -(*negep);
  203.       (*epsneg) = a;
  204.  
  205. /*
  206.   determine machep, eps
  207. */
  208.  
  209.       (*machep) = -(*it) - 3;
  210.       a = b;
  211.       tmp = one+a;
  212.  
  213.       while (tmp-one == zero) {
  214.          a = a*beta;
  215.          *machep = *machep+1;
  216.          tmp = one+a;
  217.       }
  218.  
  219.       *eps = a;
  220.       
  221. /*
  222.   determine ngrd
  223. */
  224.  
  225.       (*ngrd) = 0;
  226.       tmp = one+*eps;
  227.       tmp = tmp*one;
  228.       if (((*irnd) == 0) && (tmp-one) != zero) (*ngrd) = 1;
  229.  
  230. /*
  231.   determine iexp, minexp, xmin
  232.  
  233.   loop to determine largest i such that
  234.          (1/beta) ** (2**(i))
  235.     does not underflow.
  236.     exit from loop is signaled by an underflow.
  237. */
  238.  
  239.       i = 0;
  240.       k = 1;
  241.       z = betain;
  242.       t = one+*eps;
  243.       nxres = 0;
  244.  
  245.       for (;;) {
  246.          y = z;
  247.          z = y * y;
  248.  
  249. /*
  250.   check for underflow
  251. */
  252.  
  253.          a = z * one;
  254.          tmp = z*t;
  255.          if ((a+a == zero) || (ABS(z) > y)) break;
  256.          tmp1 = tmp*betain;
  257.          if (tmp1*beta == z) break;
  258.          i = i + 1;
  259.          k = k+k;
  260.       }
  261.  
  262. /*
  263.   determine k such that (1/beta)**k does not underflow
  264.     first set  k = 2 ** i
  265. */
  266.  
  267.       (*iexp) = i + 1;
  268.       mx = k + k;
  269.       if (*ibeta == 10) {
  270.  
  271. /*
  272.   for decimal machines only
  273. */
  274.  
  275.          (*iexp) = 2;
  276.          iz = *ibeta;
  277.          while (k >= iz) {
  278.             iz = iz * (*ibeta);
  279.             (*iexp) = (*iexp) + 1;
  280.          }
  281.          mx = iz + iz - 1;
  282.       }
  283.  
  284. /*
  285.   loop to determine minexp, xmin.
  286.     exit from loop is signaled by an underflow.
  287. */
  288.  
  289.       for (;;) {
  290.          (*xmin) = y;
  291.          y = y * betain;
  292.          a = y * one;
  293.          tmp = y*t;
  294.          tmp1 = a+a;
  295.          if ((tmp1 == zero) || (ABS(y) >= (*xmin))) break;
  296.          k = k + 1;
  297.          tmp1 = tmp*betain;
  298.          tmp1 = tmp1*beta;
  299.  
  300.          if ((tmp1 == y) && (tmp != y)) {
  301.             nxres = 3;
  302.             *xmin = y;
  303.             break;
  304.          }
  305.  
  306.       }
  307.  
  308.       (*minexp) = -k;
  309.  
  310. /*
  311.   determine maxexp, xmax
  312. */
  313.  
  314.       if ((mx <= k+k-3) && ((*ibeta) != 10)) {
  315.          mx = mx + mx;
  316.          (*iexp) = (*iexp) + 1;
  317.       }
  318.  
  319.       (*maxexp) = mx + (*minexp);
  320.  
  321. /*
  322.   Adjust *irnd to reflect partial underflow.
  323. */
  324.  
  325.       (*irnd) = (*irnd)+nxres;
  326.  
  327. /*
  328.   Adjust for IEEE style machines.
  329. */
  330.  
  331.       if ((*irnd) >= 2) (*maxexp) = (*maxexp)-2;
  332.  
  333. /*
  334.   adjust for machines with implicit leading bit in binary
  335.     significand and machines with radix point at extreme
  336.     right of significand.
  337. */
  338.  
  339.       i = (*maxexp) + (*minexp);
  340.       if (((*ibeta) == 2) && (i == 0)) (*maxexp) = (*maxexp) - 1;
  341.       if (i > 20) (*maxexp) = (*maxexp) - 1;
  342.       if (a != y) (*maxexp) = (*maxexp) - 2;
  343.       (*xmax) = one - (*epsneg);
  344.       tmp = (*xmax)*one;
  345.       if (tmp != (*xmax)) (*xmax) = one - beta * (*epsneg);
  346.       (*xmax) = (*xmax) / (beta * beta * beta * (*xmin));
  347.       i = (*maxexp) + (*minexp) + 3;
  348.       if (i > 0) {
  349.  
  350.          for (j = 1; j<=i; j++ ) {
  351.              if ((*ibeta) == 2) (*xmax) = (*xmax) + (*xmax);
  352.              if ((*ibeta) != 2) (*xmax) = (*xmax) * beta;
  353.          }
  354.  
  355.       }
  356.  
  357.     return;
  358.  
  359. }
  360.  
  361. typedef union
  362. {
  363.   double d;
  364.   int i[2];
  365. } equiv;
  366.  
  367. int
  368. main (void)
  369. {
  370.   /* Works for 32 bit machines with 32 bit ints and 64 bit doubles */
  371.  
  372.   int ibeta, iexp, irnd, it, machep, maxexp, minexp, negep, ngrd;
  373.   REAL eps, epsneg, xmax, xmin;
  374.   int i;
  375.   equiv flt_params[6];
  376.  
  377.   rmachar (&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp,
  378.        &maxexp, &eps, &epsneg, &xmin, &xmax);
  379.  
  380.   flt_params[1].d = xmin;
  381.   flt_params[2].d = xmax;
  382.   flt_params[3].d = epsneg;
  383.   flt_params[4].d = eps;
  384.   flt_params[5].d = log10 ((double) ibeta);
  385.  
  386.   printf ("* d1mach.f  Do not edit.  Generated automatically by gen-d1mach.c\n\
  387.       double precision function d1mach(i)\n\
  388.       integer i\n\
  389.       integer i1var (2)\n\
  390.       integer i2var (2)\n\
  391.       integer i3var (2)\n\
  392.       integer i4var (2)\n\
  393.       integer i5var (2)\n\
  394.       double precision dmach(5)\n\
  395.       equivalence (dmach(1), i1var(1))\n\
  396.       equivalence (dmach(2), i2var(1))\n\
  397.       equivalence (dmach(3), i3var(1))\n\
  398.       equivalence (dmach(4), i4var(1))\n\
  399.       equivalence (dmach(5), i5var(1))\n");
  400.  
  401.   for (i = 1; i < 6; i++)
  402.     printf ("      data i%dvar(1), i%dvar(2) / %ld , %ld /\n",
  403.         i, i, flt_params[i].i[0], flt_params[i].i[1]);
  404.  
  405.   printf ("      if (i .lt. 1  .or.  i .gt. 5) goto 999\n\
  406.       d1mach = dmach(i)\n\
  407.       return\n\
  408.   999 write(*,1999) i\n\
  409.  1999 format(' d1mach - i out of bounds', i10)\n\
  410.       call xstopx (' ')\n\
  411.       d1mach = 0\n\
  412.       end\n");
  413.  
  414.   return 0;
  415. }
  416.