changeset 981:73cc20a6976b

[project @ 1994-12-14 03:46:24 by jwe] Initial revision
author jwe
date Wed, 14 Dec 1994 03:48:48 +0000
parents 4793e60ad17c
children db38fe433efd
files libcruft/misc/d1mach-tst.for libcruft/misc/gen-d1mach.c libcruft/misc/i1mach.f libcruft/misc/xstopx.f libcruft/ranlib/ignuin.f
diffstat 5 files changed, 1060 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
new file mode 100644
--- /dev/null
+++ b/libcruft/misc/d1mach-tst.for
@@ -0,0 +1,7 @@
+      program main
+      integer i
+      double precision d1mach
+      do 10 i = 1, 5
+        print *, d1mach (i)
+   10 continue
+      end
new file mode 100644
--- /dev/null
+++ b/libcruft/misc/gen-d1mach.c
@@ -0,0 +1,414 @@
+/*
+
+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\
+      end\n");
+
+  return 0;
+}
new file mode 100644
--- /dev/null
+++ b/libcruft/misc/i1mach.f
@@ -0,0 +1,527 @@
+      INTEGER FUNCTION I1MACH(I)
+C
+C  I/O UNIT NUMBERS.
+C
+C    I1MACH( 1) = THE STANDARD INPUT UNIT.
+C
+C    I1MACH( 2) = THE STANDARD OUTPUT UNIT.
+C
+C    I1MACH( 3) = THE STANDARD PUNCH UNIT.
+C
+C    I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT.
+C
+C  WORDS.
+C
+C    I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT.
+C
+C    I1MACH( 6) = THE NUMBER OF CHARACTERS PER CHARACTER STORAGE UNIT.
+C                 FOR FORTRAN 77, THIS IS ALWAYS 1.  FOR FORTRAN 66,
+C                 CHARACTER STORAGE UNIT = INTEGER STORAGE UNIT.
+C
+C  INTEGERS.
+C
+C    ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM
+C
+C               SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
+C
+C               WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1.
+C
+C    I1MACH( 7) = A, THE BASE.
+C
+C    I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS.
+C
+C    I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE.
+C
+C  FLOATING-POINT NUMBERS.
+C
+C    ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT,
+C    BASE-B FORM
+C
+C               SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
+C
+C               WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T,
+C               0 .LT. X(1), AND EMIN .LE. E .LE. EMAX.
+C
+C    I1MACH(10) = B, THE BASE.
+C
+C  SINGLE-PRECISION
+C
+C    I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
+C
+C    I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
+C
+C    I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
+C
+C  DOUBLE-PRECISION
+C
+C    I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
+C
+C    I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
+C
+C    I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
+C
+C  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
+C  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
+C  REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
+C  I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
+C  WITH THE LOCAL OPERATING SYSTEM.  FOR FORTRAN 77, YOU MAY WISH
+C  TO ADJUST THE DATA STATEMENT SO IMACH(6) IS SET TO 1, AND
+C  THEN TO COMMENT OUT THE EXECUTABLE TEST ON I .EQ. 6 BELOW.
+C  ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED.
+C  (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.)
+C
+C  FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), THE FIRST
+C  SET OF CONSTANTS BELOW SHOULD BE APPROPRIATE, EXCEPT PERHAPS
+C  FOR IMACH(1) - IMACH(4).
+C
+      INTEGER IMACH(16),OUTPUT,SANITY
+C
+      EQUIVALENCE (IMACH(4),OUTPUT)
+C
+C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
+C     3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
+C     PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300).
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   32 /
+C      DATA IMACH( 6) /    4 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   31 /
+C      DATA IMACH( 9) / 2147483647 /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   24 /
+C      DATA IMACH(12) / -125 /
+C      DATA IMACH(13) /  128 /
+C      DATA IMACH(14) /   53 /
+C      DATA IMACH(15) / -1021 /
+C      DATA IMACH(16) /  1024 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
+C
+C      DATA IMACH( 1) /   5 /
+C      DATA IMACH( 2) /   6 /
+C      DATA IMACH( 3) /   7 /
+C      DATA IMACH( 4) /   6 /
+C      DATA IMACH( 5) /  32 /
+C      DATA IMACH( 6) /   4 /
+C      DATA IMACH( 7) /   2 /
+C      DATA IMACH( 8) /  31 /
+C      DATA IMACH( 9) / 2147483647 /
+C      DATA IMACH(10) /  16 /
+C      DATA IMACH(11) /   6 /
+C      DATA IMACH(12) / -64 /
+C      DATA IMACH(13) /  63 /
+C      DATA IMACH(14) /  14 /
+C      DATA IMACH(15) / -64 /
+C      DATA IMACH(16) /  63 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
+C
+C      DATA IMACH( 1) /    7 /
+C      DATA IMACH( 2) /    2 /
+C      DATA IMACH( 3) /    2 /
+C      DATA IMACH( 4) /    2 /
+C      DATA IMACH( 5) /   36 /
+C      DATA IMACH( 6) /    4 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   33 /
+C      DATA IMACH( 9) / Z1FFFFFFFF /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   24 /
+C      DATA IMACH(12) / -256 /
+C      DATA IMACH(13) /  255 /
+C      DATA IMACH(14) /   60 /
+C      DATA IMACH(15) / -256 /
+C      DATA IMACH(16) /  255 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
+C
+C      DATA IMACH( 1) /   5 /
+C      DATA IMACH( 2) /   6 /
+C      DATA IMACH( 3) /   7 /
+C      DATA IMACH( 4) /   6 /
+C      DATA IMACH( 5) /  48 /
+C      DATA IMACH( 6) /   6 /
+C      DATA IMACH( 7) /   2 /
+C      DATA IMACH( 8) /  39 /
+C      DATA IMACH( 9) / O0007777777777777 /
+C      DATA IMACH(10) /   8 /
+C      DATA IMACH(11) /  13 /
+C      DATA IMACH(12) / -50 /
+C      DATA IMACH(13) /  76 /
+C      DATA IMACH(14) /  26 /
+C      DATA IMACH(15) / -50 /
+C      DATA IMACH(16) /  76 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
+C
+C      DATA IMACH( 1) /   5 /
+C      DATA IMACH( 2) /   6 /
+C      DATA IMACH( 3) /   7 /
+C      DATA IMACH( 4) /   6 /
+C      DATA IMACH( 5) /  48 /
+C      DATA IMACH( 6) /   6 /
+C      DATA IMACH( 7) /   2 /
+C      DATA IMACH( 8) /  39 /
+C      DATA IMACH( 9) / O0007777777777777 /
+C      DATA IMACH(10) /   8 /
+C      DATA IMACH(11) /  13 /
+C      DATA IMACH(12) / -50 /
+C      DATA IMACH(13) /  76 /
+C      DATA IMACH(14) /  26 /
+C      DATA IMACH(15) / -32754 /
+C      DATA IMACH(16) /  32780 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES.
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   60 /
+C      DATA IMACH( 6) /   10 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   48 /
+C      DATA IMACH( 9) / 00007777777777777777B /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   47 /
+C      DATA IMACH(12) / -929 /
+C      DATA IMACH(13) / 1070 /
+C      DATA IMACH(14) /   94 /
+C      DATA IMACH(15) / -929 /
+C      DATA IMACH(16) / 1069 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES.
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   60 /
+C      DATA IMACH( 6) /   10 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   48 /
+C      DATA IMACH( 9) / O"00007777777777777777" /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   47 /
+C      DATA IMACH(12) / -929 /
+C      DATA IMACH(13) / 1070 /
+C      DATA IMACH(14) /   94 /
+C      DATA IMACH(15) / -929 /
+C      DATA IMACH(16) / 1069 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR CONVEX C-1.
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   32 /
+C      DATA IMACH( 6) /    4 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   31 /
+C      DATA IMACH( 9) / 2147483647 /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   24 /
+C      DATA IMACH(12) / -128 /
+C      DATA IMACH(13) /  127 /
+C      DATA IMACH(14) /   53 /
+C      DATA IMACH(15) /-1024 /
+C      DATA IMACH(16) / 1023 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
+C
+C      DATA IMACH( 1) /     5 /
+C      DATA IMACH( 2) /     6 /
+C      DATA IMACH( 3) /   102 /
+C      DATA IMACH( 4) /     6 /
+C      DATA IMACH( 5) /    64 /
+C      DATA IMACH( 6) /     8 /
+C      DATA IMACH( 7) /     2 /
+C      DATA IMACH( 8) /    63 /
+C      DATA IMACH( 9) /  777777777777777777777B /
+C      DATA IMACH(10) /     2 /
+C      DATA IMACH(11) /    47 /
+C      DATA IMACH(12) / -8189 /
+C      DATA IMACH(13) /  8190 /
+C      DATA IMACH(14) /    94 /
+C      DATA IMACH(15) / -8099 /
+C      DATA IMACH(16) /  8190 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200.
+C
+C      DATA IMACH( 1) /   11 /
+C      DATA IMACH( 2) /   12 /
+C      DATA IMACH( 3) /    8 /
+C      DATA IMACH( 4) /   10 /
+C      DATA IMACH( 5) /   16 /
+C      DATA IMACH( 6) /    2 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   15 /
+C      DATA IMACH( 9) /32767 /
+C      DATA IMACH(10) /   16 /
+C      DATA IMACH(11) /    6 /
+C      DATA IMACH(12) /  -64 /
+C      DATA IMACH(13) /   63 /
+C      DATA IMACH(14) /   14 /
+C      DATA IMACH(15) /  -64 /
+C      DATA IMACH(16) /   63 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7.
+C
+C      DATA IMACH( 1) /       5 /
+C      DATA IMACH( 2) /       6 /
+C      DATA IMACH( 3) /       0 /
+C      DATA IMACH( 4) /       6 /
+C      DATA IMACH( 5) /      24 /
+C      DATA IMACH( 6) /       3 /
+C      DATA IMACH( 7) /       2 /
+C      DATA IMACH( 8) /      23 /
+C      DATA IMACH( 9) / 8388607 /
+C      DATA IMACH(10) /       2 /
+C      DATA IMACH(11) /      23 /
+C      DATA IMACH(12) /    -127 /
+C      DATA IMACH(13) /     127 /
+C      DATA IMACH(14) /      38 /
+C      DATA IMACH(15) /    -127 /
+C      DATA IMACH(16) /     127 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /   43 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   36 /
+C      DATA IMACH( 6) /    4 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   35 /
+C      DATA IMACH( 9) / O377777777777 /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   27 /
+C      DATA IMACH(12) / -127 /
+C      DATA IMACH(13) /  127 /
+C      DATA IMACH(14) /   63 /
+C      DATA IMACH(15) / -127 /
+C      DATA IMACH(16) /  127 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
+C     THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86.
+C
+C      DATA IMACH( 1) /   5 /
+C      DATA IMACH( 2) /   6 /
+C      DATA IMACH( 3) /   7 /
+C      DATA IMACH( 4) /   6 /
+C      DATA IMACH( 5) /  32 /
+C      DATA IMACH( 6) /   4 /
+C      DATA IMACH( 7) /   2 /
+C      DATA IMACH( 8) /  31 /
+C      DATA IMACH( 9) / Z7FFFFFFF /
+C      DATA IMACH(10) /  16 /
+C      DATA IMACH(11) /   6 /
+C      DATA IMACH(12) / -64 /
+C      DATA IMACH(13) /  63 /
+C      DATA IMACH(14) /  14 /
+C      DATA IMACH(15) / -64 /
+C      DATA IMACH(16) /  63 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE INTERDATA 8/32
+C     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
+C
+C     FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
+C     THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.
+C
+C      DATA IMACH( 1) /   5 /
+C      DATA IMACH( 2) /   6 /
+C      DATA IMACH( 3) /   6 /
+C      DATA IMACH( 4) /   6 /
+C      DATA IMACH( 5) /  32 /
+C      DATA IMACH( 6) /   4 /
+C      DATA IMACH( 7) /   2 /
+C      DATA IMACH( 8) /  31 /
+C      DATA IMACH( 9) / Z'7FFFFFFF' /
+C      DATA IMACH(10) /  16 /
+C      DATA IMACH(11) /   6 /
+C      DATA IMACH(12) / -64 /
+C      DATA IMACH(13) /  62 /
+C      DATA IMACH(14) /  14 /
+C      DATA IMACH(15) / -64 /
+C      DATA IMACH(16) /  62 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   36 /
+C      DATA IMACH( 6) /    5 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   35 /
+C      DATA IMACH( 9) / "377777777777 /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   27 /
+C      DATA IMACH(12) / -128 /
+C      DATA IMACH(13) /  127 /
+C      DATA IMACH(14) /   54 /
+C      DATA IMACH(15) / -101 /
+C      DATA IMACH(16) /  127 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   36 /
+C      DATA IMACH( 6) /    5 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   35 /
+C      DATA IMACH( 9) / "377777777777 /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   27 /
+C      DATA IMACH(12) / -128 /
+C      DATA IMACH(13) /  127 /
+C      DATA IMACH(14) /   62 /
+C      DATA IMACH(15) / -128 /
+C      DATA IMACH(16) /  127 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
+C     32-BIT INTEGER ARITHMETIC.
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   32 /
+C      DATA IMACH( 6) /    4 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   31 /
+C      DATA IMACH( 9) / 2147483647 /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   24 /
+C      DATA IMACH(12) / -127 /
+C      DATA IMACH(13) /  127 /
+C      DATA IMACH(14) /   56 /
+C      DATA IMACH(15) / -127 /
+C      DATA IMACH(16) /  127 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
+C     16-BIT INTEGER ARITHMETIC.
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   16 /
+C      DATA IMACH( 6) /    2 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   15 /
+C      DATA IMACH( 9) / 32767 /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   24 /
+C      DATA IMACH(12) / -127 /
+C      DATA IMACH(13) /  127 /
+C      DATA IMACH(14) /   56 /
+C      DATA IMACH(15) / -127 /
+C      DATA IMACH(16) /  127 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS
+C     WTIH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS,
+C     SUPPLIED BY IGOR BRAY.
+C
+C      DATA IMACH( 1) /            1 /
+C      DATA IMACH( 2) /            1 /
+C      DATA IMACH( 3) /            2 /
+C      DATA IMACH( 4) /            1 /
+C      DATA IMACH( 5) /           32 /
+C      DATA IMACH( 6) /            4 /
+C      DATA IMACH( 7) /            2 /
+C      DATA IMACH( 8) /           31 /
+C      DATA IMACH( 9) / :17777777777 /
+C      DATA IMACH(10) /            2 /
+C      DATA IMACH(11) /           23 /
+C      DATA IMACH(12) /         -127 /
+C      DATA IMACH(13) /         +127 /
+C      DATA IMACH(14) /           47 /
+C      DATA IMACH(15) /       -32895 /
+C      DATA IMACH(16) /       +32637 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
+C
+C      DATA IMACH( 1) /     0 /
+C      DATA IMACH( 2) /     0 /
+C      DATA IMACH( 3) /     7 /
+C      DATA IMACH( 4) /     0 /
+C      DATA IMACH( 5) /    32 /
+C      DATA IMACH( 6) /     1 /
+C      DATA IMACH( 7) /     2 /
+C      DATA IMACH( 8) /    31 /
+C      DATA IMACH( 9) /  2147483647 /
+C      DATA IMACH(10) /     2 /
+C      DATA IMACH(11) /    24 /
+C      DATA IMACH(12) /  -125 /
+C      DATA IMACH(13) /   128 /
+C      DATA IMACH(14) /    53 /
+C      DATA IMACH(15) / -1021 /
+C      DATA IMACH(16) /  1024 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
+C
+C     NOTE THAT THE PUNCH UNIT, I1MACH(3), HAS BEEN SET TO 7
+C     WHICH IS APPROPRIATE FOR THE UNIVAC-FOR SYSTEM.
+C     IF YOU HAVE THE UNIVAC-FTN SYSTEM, SET IT TO 1.
+C
+C      DATA IMACH( 1) /    5 /
+C      DATA IMACH( 2) /    6 /
+C      DATA IMACH( 3) /    7 /
+C      DATA IMACH( 4) /    6 /
+C      DATA IMACH( 5) /   36 /
+C      DATA IMACH( 6) /    6 /
+C      DATA IMACH( 7) /    2 /
+C      DATA IMACH( 8) /   35 /
+C      DATA IMACH( 9) / O377777777777 /
+C      DATA IMACH(10) /    2 /
+C      DATA IMACH(11) /   27 /
+C      DATA IMACH(12) / -128 /
+C      DATA IMACH(13) /  127 /
+C      DATA IMACH(14) /   60 /
+C      DATA IMACH(15) /-1024 /
+C      DATA IMACH(16) / 1023 /, SANITY/987/
+C
+C     MACHINE CONSTANTS FOR VAX.
+C
+      DATA IMACH( 1) /    5 /
+      DATA IMACH( 2) /    6 /
+      DATA IMACH( 3) /    7 /
+      DATA IMACH( 4) /    6 /
+      DATA IMACH( 5) /   32 /
+      DATA IMACH( 6) /    4 /
+      DATA IMACH( 7) /    2 /
+      DATA IMACH( 8) /   31 /
+      DATA IMACH( 9) / 2147483647 /
+      DATA IMACH(10) /    2 /
+      DATA IMACH(11) /   24 /
+      DATA IMACH(12) / -127 /
+      DATA IMACH(13) /  127 /
+      DATA IMACH(14) /   56 /
+      DATA IMACH(15) / -127 /
+      DATA IMACH(16) /  127 /, SANITY/987/
+C
+C  ***  ISSUE STOP 777 IF ALL DATA STATEMENTS ARE COMMENTED...
+      IF (SANITY .NE. 987) STOP 777
+      IF (I .LT. 1  .OR.  I .GT. 16) GO TO 10
+C
+      I1MACH = IMACH(I)
+C/6S
+C/7S
+      IF(I.EQ.6) I1MACH=1
+C/
+      RETURN
+   10 WRITE(OUTPUT,1999) I
+ 1999 FORMAT(' I1MACH - I OUT OF BOUNDS',I10)
+      CALL XSTOPX (' ')
+      END
new file mode 100644
--- /dev/null
+++ b/libcruft/misc/xstopx.f
@@ -0,0 +1,9 @@
+      subroutine xstopx (string)
+      character *(*) string
+      integer slen
+      slen = len (string)
+      if (slen .eq. 0) goto 9999
+      if (slen .eq. 1 .and. string(1:1) .eq. ' ') goto 9999
+      write (*, *) string
+ 9999 call dostop ()
+      end
new file mode 100644
--- /dev/null
+++ b/libcruft/ranlib/ignuin.f
@@ -0,0 +1,103 @@
+      INTEGER FUNCTION ignuin(low,high)
+C**********************************************************************
+C
+C     INTEGER FUNCTION IGNUIN( LOW, HIGH )
+C
+C               GeNerate Uniform INteger
+C
+C
+C                              Function
+C
+C
+C     Generates an integer uniformly distributed between LOW and HIGH.
+C
+C
+C                              Arguments
+C
+C
+C     LOW --> Low bound (inclusive) on integer value to be generated
+C                         INTEGER LOW
+C
+C     HIGH --> High bound (inclusive) on integer value to be generated
+C                         INTEGER HIGH
+C
+C
+C                              Note
+C
+C
+C     If (HIGH-LOW) > 2,147,483,561 prints error message on * unit and
+C     stops the program.
+C
+C**********************************************************************
+
+C     IGNLGI generates integers between 1 and 2147483562
+C     MAXNUM is 1 less than maximum generable value
+C     .. Parameters ..
+      INTEGER maxnum
+      PARAMETER (maxnum=2147483561)
+      CHARACTER*(*) err1,err2
+      PARAMETER (err1='LOW > HIGH in IGNUIN',
+     +          err2=' ( HIGH - LOW ) > 2,147,483,561 in IGNUIN')
+C     ..
+C     .. Scalar Arguments ..
+      INTEGER high,low
+C     ..
+C     .. Local Scalars ..
+      INTEGER err,ign,maxnow,range,ranp1
+C     ..
+C     .. External Functions ..
+      INTEGER ignlgi
+      EXTERNAL ignlgi
+C     ..
+C     .. Intrinsic Functions ..
+      INTRINSIC mod
+C     ..
+C     .. Executable Statements ..
+      IF (.NOT. (low.GT.high)) GO TO 10
+      err = 1
+C      ABORT-PROGRAM
+      GO TO 80
+
+   10 range = high - low
+      IF (.NOT. (range.GT.maxnum)) GO TO 20
+      err = 2
+C      ABORT-PROGRAM
+      GO TO 80
+
+   20 IF (.NOT. (low.EQ.high)) GO TO 30
+      ignuin = low
+      RETURN
+
+      GO TO 70
+
+C     Number to be generated should be in range 0..RANGE
+C     Set MAXNOW so that the number of integers in 0..MAXNOW is an
+C     integral multiple of the number in 0..RANGE
+
+   30 ranp1 = range + 1
+      maxnow = (maxnum/ranp1)*ranp1
+   40 ign = ignlgi() - 1
+      IF (.NOT. (ign.LE.maxnow)) GO TO 50
+      ignuin = low + mod(ign,ranp1)
+      RETURN
+
+   50 GO TO 40
+
+   60 CONTINUE
+   70 CONTINUE
+   80 IF (.NOT. (err.EQ.1)) GO TO 90
+      WRITE (*,*) err1
+      GO TO 100
+
+C     TO ABORT-PROGRAM
+   90 WRITE (*,*) err2
+  100 WRITE (*,*) ' LOW: ',low,' HIGH: ',high
+      WRITE (*,*) ' Abort on Fatal ERROR'
+      IF (.NOT. (err.EQ.1)) GO TO 110
+      CALL XSTOPX ('LOW > HIGH in IGNUIN')
+
+      GO TO 120
+
+  110 STOP ' ( HIGH - LOW ) > 2,147,483,561 in IGNUIN'
+
+  120 END