diff --git a/math/ieee/README b/math/ieee/README
deleted file mode 100644
index 80b30bfd5..000000000
--- a/math/ieee/README
+++ /dev/null
@@ -1,8 +0,0 @@
-     This directory contians programs and subroutines from the
-text "Programs for Digital Signal Processing" by the IEEE Press.
-Directories of the name "chap*" correspond to each chapter of
-the text.
-     The files "*mach*" are the machine dependent routines which
-should be edited for the hardware which you are going to run the
-routines on.  The file "uni.c" is a uniform random number generator,
-not the one from the text, but one which uses the UNIX generator.
diff --git a/math/ieee/chap1/README b/math/ieee/chap1/README
deleted file mode 100644
index bf86b45f1..000000000
--- a/math/ieee/chap1/README
+++ /dev/null
@@ -1,14 +0,0 @@
-     This directory contains the IEEE routines from Chapter 1 -
-Discrete Fourier Transform Programs.  Some of the routines are
-provided here as separete files (fourea.f, wfta.f, ...) but are
-purposely not put into the library, as these routines are not for
-general use (fourea.f is an inefficient, demonstration version;
-wfta.f is the Winograd DFT which is slower than the regular DFT
-and uses far too much memory ona  16-bit machine to be of
-practical use; ...).
-     The directory "test" contains the test routines from the chapter.
-The directory "time" contains some routines to time various of the
-routines.
-     The file "compall" compiles the appropiate routines and then
-"mklib" will make the library, which one will probably want to
-move to "/usr/lib/libieee.a".
diff --git a/math/ieee/chap1/const.f b/math/ieee/chap1/const.f
deleted file mode 100644
index 54eb4e410..000000000
--- a/math/ieee/chap1/const.f
+++ /dev/null
@@ -1,114 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: const
-c   computes the multipliers for the various modules
-c-----------------------------------------------------------------------
-c
-      subroutine const(co3,co8,co16,co9,cdc,cdd)
-      double precision dtheta,dtwopi,dsq32,dsq2
-      double precision dcos1,dcos2,dcos3,dcos4
-      double precision dsin1,dsin2,dsin3,dsin4
-      dimension co3(3),co8(8),co16(18),co9(11),cdc(9),cdd(6)
-      dtwopi=8.0d0*datan(1.0d0)
-      dsq32=dsqrt(0.75d0)
-      dsq2=dsqrt(0.5d0)
-c
-c   multipliers for the three point module
-c
-      co3(1)=1.0
-      co3(2)=-1.5
-      co3(3)=-dsq32
-c
-c   multipliers for the five point module
-c
-      dtheta=dtwopi/5.0d0
-      dcos1=dcos(dtheta)
-      dcos2=dcos(2.0d0*dtheta)
-      dsin1=dsin(dtheta)
-      dsin2=dsin(2.0d0*dtheta)
-      cdd(1)=1.0
-      cdd(2)=-1.25
-      cdd(3)=-dsin1-dsin2
-      cdd(4)=0.5*(dcos1-dcos2)
-      cdd(5)=dsin1-dsin2
-      cdd(6)=dsin2
-c
-c
-c   multipliers for the seven point module
-c
-      dtheta=dtwopi/7.0d0
-      dcos1=dcos(dtheta)
-      dcos2=dcos(2.0d0*dtheta)
-      dcos3=dcos(3.0d0*dtheta)
-      dsin1=dsin(dtheta)
-      dsin2=dsin(2.0d0*dtheta)
-      dsin3=dsin(3.0d0*dtheta)
-      cdc(1)=1.0
-      cdc(2)=-7.0d0/6.0d0
-      cdc(3)=-(dsin1+dsin2-dsin3)/3.0d0
-      cdc(4)=(dcos1+dcos2-2.0d0*dcos3)/3.0d0
-      cdc(5)=(2.0d0*dcos1-dcos2-dcos3)/3.0d0
-      cdc(6)=-(2.0d0*dsin1-dsin2+dsin3)/3.0d0
-      cdc(7)=-(dsin1+dsin2+2.0d0*dsin3)/3.0d0
-      cdc(8)=(dcos1-2.0d0*dcos2+dcos3)/3.0d0
-      cdc(9)=-(dsin1-2.0d0*dsin2-dsin3)/3.0d0
-c
-c   multipliers for the eight point module
-c
-      co8(1)=1.0
-      co8(2)=1.0
-      co8(3)=1.0
-      co8(4)=-1.0
-      co8(5)=1.0
-      co8(6)=-dsq2
-      co8(7)=-1.0
-      co8(8)=dsq2
-c
-c   multipliers for the nine point module
-c
-      dtheta=dtwopi/9.0d0
-      dcos1=dcos(dtheta)
-      dcos2=dcos(2.0d0*dtheta)
-      dcos4=dcos(4.0d0*dtheta)
-      dsin1=dsin(dtheta)
-      dsin2=dsin(2.0d0*dtheta)
-      dsin4=dsin(4.0d0*dtheta)
-      co9(1)=1.0
-      co9(2)=-1.5
-      co9(3)=-dsq32
-      co9(4)=0.5
-      co9(5)=(2.0d0*dcos1-dcos2-dcos4)/3.0d0
-      co9(6)=(dcos1-2.0d0*dcos2+dcos4)/3.0d0
-      co9(7)=(dcos1+dcos2-2.0d0*dcos4)/3.0d0
-      co9(8)=-(2.0d0*dsin1+dsin2-dsin4)/3.0d0
-      co9(9)=-(dsin1+2.0d0*dsin2+dsin4)/3.0d0
-      co9(10)=-(dsin1-dsin2-2.0d0*dsin4)/3.0d0
-      co9(11)=-dsq32
-c
-c   multipliers for the sixteen point module
-c
-      dtheta=dtwopi/16.0d0
-      dcos1=dcos(dtheta)
-      dcos3=dcos(3.0d0*dtheta)
-      dsin1=dsin(dtheta)
-      dsin3=dsin(3.0d0*dtheta)
-      co16(1)=1.0
-      co16(2)=1.0
-      co16(3)=1.0
-      co16(4)=-1.0
-      co16(5)=1.0
-      co16(6)=-dsq2
-      co16(7)=-1.0
-      co16(8)=dsq2
-      co16(9)=1.0
-      co16(10)=-(dsin1-dsin3)
-      co16(11)=-dsq2
-      co16(12)=-co16(10)
-      co16(13)=-1.0
-      co16(14)=-(dsin1+dsin3)
-      co16(15)=dsq2
-      co16(16)=-co16(14)
-      co16(17)=-dsin3
-      co16(18)=dcos3
-      return
-      end
diff --git a/math/ieee/chap1/fast.f b/math/ieee/chap1/fast.f
deleted file mode 100644
index 9a0ad713f..000000000
--- a/math/ieee/chap1/fast.f
+++ /dev/null
@@ -1,73 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  fast
-c replaces the real vector b(k), for k=1,2,...,n,
-c with its finite discrete fourier transform
-c-----------------------------------------------------------------------
-c
-      subroutine fast(b, n)
-c
-c the dc term is returned in location b(1) with b(2) set to 0.
-c thereafter the jth harmonic is returned as a complex
-c number stored as  b(2*j+1) + i b(2*j+2).
-c the n/2 harmonic is returned in b(n+1) with b(n+2) set to 0.
-c hence, b must be dimensioned to size n+2.
-c the subroutine is called as  fast(b,n) where n=2**m and
-c b is the real array described above.
-c
-      dimension b(2)
-      common /cons/ pii, p7, p7two, c22, s22, pi2
-c
-c iw is a machine dependent write device number
-c
-      iw = i1mach(2)
-c
-      pii = 4.*atan(1.)
-      pi8 = pii/8.
-      p7 = 1./sqrt(2.)
-      p7two = 2.*p7
-      c22 = cos(pi8)
-      s22 = sin(pi8)
-      pi2 = 2.*pii
-      do 10 i=1,15
-        m = i
-        nt = 2**i
-        if (n.eq.nt) go to 20
-  10  continue
-      write (iw,9999)
-9999  format (33h n is not a power of two for fast)
-      stop
-  20  n4pow = m/2
-c
-c do a radix 2 iteration first if one is required.
-c
-      if (m-n4pow*2) 40, 40, 30
-  30  nn = 2
-      int = n/nn
-      call fr2tr(int, b(1), b(int+1))
-      go to 50
-  40  nn = 1
-c
-c perform radix 4 iterations.
-c
-  50  if (n4pow.eq.0) go to 70
-      do 60 it=1,n4pow
-        nn = nn*4
-        int = n/nn
-        call fr4tr(int, nn, b(1), b(int+1), b(2*int+1), b(3*int+1),
-     *      b(1), b(int+1), b(2*int+1), b(3*int+1))
-  60  continue
-c
-c perform in-place reordering.
-c
-  70  call ford1(m, b)
-      call ford2(m, b)
-      t = b(2)
-      b(2) = 0.
-      b(n+1) = t
-      b(n+2) = 0.
-      do 80 it=4,n,2
-        b(it) = -b(it)
-  80  continue
-      return
-      end
diff --git a/math/ieee/chap1/ffa.f b/math/ieee/chap1/ffa.f
deleted file mode 100644
index dd62ea01e..000000000
--- a/math/ieee/chap1/ffa.f
+++ /dev/null
@@ -1,84 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  ffa
-c fast fourier analysis subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine ffa(b, nfft)
-c
-c this subroutine replaces the real vector b(k),  (k=1,2,...,n),
-c with its finite discrete fourier transform.  the dc term is
-c returned in location b(1) with b(2) set to 0.  thereafter, the
-c jth harmonic is returned as a complex number stored as
-c b(2*j+1) + i b(2*j+2).  note that the n/2 harmonic is returned
-c in b(n+1) with b(n+2) set to 0.  hence, b must be dimensioned
-c to size n+2.
-c subroutine is called as ffa (b,n) where n=2**m and b is an
-c n term real array.  a real-valued, radix 8  algorithm is used
-c with in-place reordering and the trig functions are computed as
-c needed.
-c
-      dimension b(2)
-      common /con/ pii, p7, p7two, c22, s22, pi2
-c
-c iw is a machine dependent write device number
-c
-      iw = i1mach(2)
-c
-      pii = 4.*atan(1.)
-      pi8 = pii/8.
-      p7 = 1./sqrt(2.)
-      p7two = 2.*p7
-      c22 = cos(pi8)
-      s22 = sin(pi8)
-      pi2 = 2.*pii
-      n = 1
-      do 10 i=1,15
-        m = i
-        n = n*2
-        if (n.eq.nfft) go to 20
-  10  continue
-      write (iw,9999)
-9999  format (30h nfft not a power of 2 for ffa)
-      stop
-  20  continue
-      n8pow = m/3
-c
-c do a radix 2 or radix 4 iteration first if one is required
-c
-      if (m-n8pow*3-1) 50, 40, 30
-  30  nn = 4
-      int = n/nn
-      call r4tr(int, b(1), b(int+1), b(2*int+1), b(3*int+1))
-      go to 60
-  40  nn = 2
-      int = n/nn
-      call r2tr(int, b(1), b(int+1))
-      go to 60
-  50  nn = 1
-c
-c perform radix 8 iterations
-c
-  60  if (n8pow) 90, 90, 70
-  70  do 80 it=1,n8pow
-        nn = nn*8
-        int = n/nn
-        call r8tr(int, nn, b(1), b(int+1), b(2*int+1), b(3*int+1),
-     *      b(4*int+1), b(5*int+1), b(6*int+1), b(7*int+1), b(1),
-     *      b(int+1), b(2*int+1), b(3*int+1), b(4*int+1), b(5*int+1),
-     *      b(6*int+1), b(7*int+1))
-  80  continue
-c
-c perform in-place reordering
-c
-  90  call ord1(m, b)
-      call ord2(m, b)
-      t = b(2)
-      b(2) = 0.
-      b(nfft+1) = t
-      b(nfft+2) = 0.
-      do 100 i=4,nfft,2
-        b(i) = -b(i)
- 100  continue
-      return
-      end
diff --git a/math/ieee/chap1/ffs.f b/math/ieee/chap1/ffs.f
deleted file mode 100644
index 2858fdb07..000000000
--- a/math/ieee/chap1/ffs.f
+++ /dev/null
@@ -1,80 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  ffs
-c fast fourier synthesis subroutine
-c radix 8-4-2
-c-----------------------------------------------------------------------
-c
-      subroutine ffs(b, nfft)
-c
-c this subroutine synthesizes the real vector b(k), where
-c k=1,2,...,n. the initial fourier coefficients are placed in
-c the b array of size n+2.  the dc term is in b(1) with
-c b(2) equal to 0.
-c the jth harmonic is stored as b(2*j+1) + i b(2*j+2).
-c the n/2 harmonic is in b(n+1) with b(n+2) equal to 0.
-c the subroutine is called as ffs(b,n) where n=2**m and
-c b is the n term real array discussed above.
-c
-      dimension b(2)
-      common /con1/ pii, p7, p7two, c22, s22, pi2
-c
-c iw is a machine dependent write device number
-c
-      iw = i1mach(2)
-c
-      pii = 4.*atan(1.)
-      pi8 = pii/8.
-      p7 = 1./sqrt(2.)
-      p7two = 2.*p7
-      c22 = cos(pi8)
-      s22 = sin(pi8)
-      pi2 = 2.*pii
-      n = 1
-      do 10 i=1,15
-        m = i
-        n = n*2
-        if (n.eq.nfft) go to 20
-  10  continue
-      write (iw,9999)
-9999  format (30h nfft not a power of 2 for ffs)
-      stop
-  20  continue
-      b(2) = b(nfft+1)
-      do 30 i=1,nfft
-        b(i) = b(i)/float(nfft)
-  30  continue
-      do 40 i=4,nfft,2
-        b(i) = -b(i)
-  40  continue
-      n8pow = m/3
-c
-c reorder the input fourier coefficients
-c
-      call ord2(m, b)
-      call ord1(m, b)
-c
-      if (n8pow.eq.0) go to 60
-c
-c perform the radix 8 iterations
-c
-      nn = n
-      do 50 it=1,n8pow
-        int = n/nn
-        call r8syn(int, nn, b, b(int+1), b(2*int+1), b(3*int+1),
-     *      b(4*int+1), b(5*int+1), b(6*int+1), b(7*int+1), b(1),
-     *      b(int+1), b(2*int+1), b(3*int+1), b(4*int+1), b(5*int+1),
-     *      b(6*int+1), b(7*int+1))
-        nn = nn/8
-  50  continue
-c
-c do a radix 2 or radix 4 iteration if one is required
-c
-  60  if (m-n8pow*3-1) 90, 80, 70
-  70  int = n/4
-      call r4syn(int, b(1), b(int+1), b(2*int+1), b(3*int+1))
-      go to 90
-  80  int = n/2
-      call r2tr(int, b(1), b(int+1))
-  90  return
-      end
diff --git a/math/ieee/chap1/fft842.f b/math/ieee/chap1/fft842.f
deleted file mode 100644
index 37a3a651c..000000000
--- a/math/ieee/chap1/fft842.f
+++ /dev/null
@@ -1,116 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  fft842
-c fast fourier transform for n=2**m
-c complex input
-c-----------------------------------------------------------------------
-c
-      subroutine fft842(in, n, x, y)
-c
-c this program replaces the vector z=x+iy by its  finite
-c discrete, complex fourier transform if in=0.  the inverse transform
-c is calculated for in=1.  it performs as many base
-c 8 iterations as possible and then finishes with a base 4 iteration
-c or a base 2 iteration if needed.
-c
-c the subroutine is called as subroutine fft842 (in,n,x,y).
-c the integer n (a power of 2), the n real location array x, and
-c the n real location array y must be supplied to the subroutine.
-c
-      dimension x(2), y(2), l(15)
-      common /con2/ pi2, p7
-      equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
-     *    (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
-     *    (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
-     *    (l1,l(15))
-c
-c
-c iw is a machine dependent write device number
-c
-      iw = i1mach(2)
-c
-      pi2 = 8.*atan(1.)
-      p7 = 1./sqrt(2.)
-      do 10 i=1,15
-        m = i
-        nt = 2**i
-        if (n.eq.nt) go to 20
-  10  continue
-      write (iw,9999)
-9999  format (35h n is not a power of two for fft842)
-      stop
-  20  n2pow = m
-      nthpo = n
-      fn = nthpo
-      if (in.eq.1) go to 40
-      do 30 i=1,nthpo
-        y(i) = -y(i)
-  30  continue
-  40  n8pow = n2pow/3
-      if (n8pow.eq.0) go to 60
-c
-c radix 8 passes,if any.
-c
-      do 50 ipass=1,n8pow
-        nxtlt = 2**(n2pow-3*ipass)
-        lengt = 8*nxtlt
-        call r8tx(nxtlt, nthpo, lengt, x(1), x(nxtlt+1), x(2*nxtlt+1),
-     *      x(3*nxtlt+1), x(4*nxtlt+1), x(5*nxtlt+1), x(6*nxtlt+1),
-     *      x(7*nxtlt+1), y(1), y(nxtlt+1), y(2*nxtlt+1), y(3*nxtlt+1),
-     *      y(4*nxtlt+1), y(5*nxtlt+1), y(6*nxtlt+1), y(7*nxtlt+1))
-  50  continue
-c
-c is there a four factor left
-c
-  60  if (n2pow-3*n8pow-1) 90, 70, 80
-c
-c go through the base 2 iteration
-c
-c
-  70  call r2tx(nthpo, x(1), x(2), y(1), y(2))
-      go to 90
-c
-c go through the base 4 iteration
-c
-  80  call r4tx(nthpo, x(1), x(2), x(3), x(4), y(1), y(2), y(3), y(4))
-c
-  90  do 110 j=1,15
-        l(j) = 1
-        if (j-n2pow) 100, 100, 110
- 100    l(j) = 2**(n2pow+1-j)
- 110  continue
-      ij = 1
-      do 130 j1=1,l1
-      do 130 j2=j1,l2,l1
-      do 130 j3=j2,l3,l2
-      do 130 j4=j3,l4,l3
-      do 130 j5=j4,l5,l4
-      do 130 j6=j5,l6,l5
-      do 130 j7=j6,l7,l6
-      do 130 j8=j7,l8,l7
-      do 130 j9=j8,l9,l8
-      do 130 j10=j9,l10,l9
-      do 130 j11=j10,l11,l10
-      do 130 j12=j11,l12,l11
-      do 130 j13=j12,l13,l12
-      do 130 j14=j13,l14,l13
-      do 130 ji=j14,l15,l14
-        if (ij-ji) 120, 130, 130
- 120    r = x(ij)
-        x(ij) = x(ji)
-        x(ji) = r
-        fi = y(ij)
-        y(ij) = y(ji)
-        y(ji) = fi
- 130    ij = ij + 1
-      if (in.eq.1) go to 150
-      do 140 i=1,nthpo
-        y(i) = -y(i)
- 140  continue
-      go to 170
- 150  do 160 i=1,nthpo
-        x(i) = x(i)/fn
-        y(i) = y(i)/fn
- 160  continue
- 170  return
-      end
diff --git a/math/ieee/chap1/fftaoh.f b/math/ieee/chap1/fftaoh.f
deleted file mode 100644
index f703c98bb..000000000
--- a/math/ieee/chap1/fftaoh.f
+++ /dev/null
@@ -1,82 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: fftaoh
-c compute dft for real, antisymmetric, odd harmonic, n-point sequence
-c using n/4-point fft
-c antisymmetric sequence means x(m)=-x(n-m), m=1,...,n/2-1
-c odd harmonic means x(2*k)=0, all k, where x(k) is the dft of x(m)
-c x(m) has the property x(m)=x(n/2-m), m=0,1,...,n/4-1,  x(0)=0
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine fftaoh(x, n, y)
-      dimension x(1), y(1)
-c
-c  x = real array which on input contains the (n/4+1) points of the
-c      input sequence (antisymmetrical)
-c      on output x contains the n/4 imaginary points of the odd
-c      harmonics of the transform of the input--i.e. the zero
-c      valued real parts are not given nor are the zero-valued
-c      even harmonics
-c  n = true size of input
-c  y = scratch array of size n/4+2
-c
-c
-c handle n = 2 and n = 4 cases separately
-c
-      if (n.gt.4) go to 20
-      if (n.eq.4) go to 10
-c
-c for n=2, assume x(1)=0, x(2)=0, compute dft directly
-c
-      x(1) = 0.
-      return
-c
-c n = 4 case, assume x(1)=x(3)=0, x(2)=-x(4)=x0, compute dft directly
-c
-  10  x(1) = -2.*x(2)
-      return
-  20  twopi = 8.*atan(1.0)
-c
-c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1))
-c
-      no2 = n/2
-      no4 = n/4
-      no8 = n/8
-      if (no8.eq.1) go to 40
-      do 30 i=2,no8
-        ind = 2*i
-        t1 = x(ind) - x(ind-2)
-        y(i) = x(ind-1) + t1
-        ind1 = n/4 + 2 - i
-        y(ind1) = x(ind-1) - t1
-  30  continue
-  40  y(1) = 2.*x(2)
-      y(no8+1) = x(no4+1)
-c
-c the sequence y (n/4 points) has only odd harmonics
-c call subroutine fftohm to exploit odd harmonics
-c
-      call fftohm(y, no2)
-c
-c form original dft from complex odd harmonics of y(k)
-c by unscrambling y(k)
-c
-      tpn = twopi/float(n)
-      cosi = 2.*cos(tpn)
-      sini = 2.*sin(tpn)
-      cosd = cos(tpn*2.)
-      sind = sin(tpn*2.)
-      do 50 i=1,no8
-        ind = 2*i
-        bk = y(ind-1)/sini
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-        ak = y(ind)
-        x(i) = ak - bk
-        ind1 = n/4 + 1 - i
-        x(ind1) = -ak - bk
-  50  continue
-      return
-      end
diff --git a/math/ieee/chap1/fftasm.f b/math/ieee/chap1/fftasm.f
deleted file mode 100644
index 30c84b794..000000000
--- a/math/ieee/chap1/fftasm.f
+++ /dev/null
@@ -1,67 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: fftasm
-c compute dft for real, antisymmetric, n-point sequence x(m) using
-c n/2-point fft
-c antisymmetric sequence means x(m)=-x(n-m), m=1,...,n/2-1
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine fftasm(x, n, y)
-      dimension x(1), y(1)
-c
-c x = real array which on input contains the n/2 points of the
-c     input sequence (asymmetrical)
-c     on output x contains the n/2+1 imaginary points of the transform
-c     of the input--i.e. the zero valued real parts are not returned
-c n = true size of input
-c y = scratch array of size n/2+2
-c
-c
-c for n = 2, assume x(1)=0, x(2)=0, compute dft directly
-c
-      if (n.eq.2) go to 30
-      twopi = 8.*atan(1.0)
-c
-c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1))
-c
-      no2 = n/2
-      no4 = n/4
-      do 10 i=2,no4
-        ind = 2*i
-        t1 = x(ind) - x(ind-2)
-        y(i) = x(ind-1) + t1
-        ind1 = no2 + 2 - i
-        y(ind1) = -x(ind-1) + t1
-  10  continue
-      y(1) = 2.*x(2)
-      y(no4+1) = -2.*x(no2)
-c
-c take n/2 point (real) fft of y
-c
-      call fast(y, no2)
-c
-c form original dft by unscrambling y(k)
-c use recursion relation to generate sin(tpn*i) multiplier
-c
-      tpn = twopi/float(n)
-      cosi = 2.*cos(tpn)
-      sini = 2.*sin(tpn)
-      cosd = cosi/2.
-      sind = sini/2.
-      nind = no4 + 1
-      do 20 i=2,nind
-        ind = 2*i
-        bk = y(ind-1)/sini
-        ak = y(ind)
-        x(i) = ak - bk
-        ind1 = no2 + 2 - i
-        x(ind1) = -ak - bk
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-  20  continue
-  30  x(1) = 0.
-      x(no2+1) = 0.
-      return
-      end
diff --git a/math/ieee/chap1/fftohm.f b/math/ieee/chap1/fftohm.f
deleted file mode 100644
index 6e3ecc9d7..000000000
--- a/math/ieee/chap1/fftohm.f
+++ /dev/null
@@ -1,101 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: fftohm
-c compute dft for real, n-point, odd harmonic sequences using an
-c n/2 point fft
-c odd harmonic means x(2*k)=0, all k where x(k) is the dft of x(m)
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine fftohm(x, n)
-      dimension x(1)
-c
-c x = real array which on input contains the first n/2 points of the
-c     input
-c     on output x contains the n/4 complex values of the odd
-c     harmonics of the input--stored in the sequence re(x(1)),im(x(1)),
-c     re(x(2)),im(x(2)),...
-c ****note: x must be dimensioned to size n/2+2 for fft routine
-c n = true size of x sequence
-c
-c first compute real(x(1)) and real(x(n/2-1)) separately
-c also simultaneously multiply original sequence by sin(twopi*(m-1)/n)
-c sin and cos are computed recursively
-c
-c
-c for n = 2, assume x(1)=x0, x(2)=-x0, compute dft directly
-c
-      if (n.gt.2) go to 10
-      x(1) = 2.*x(1)
-      x(2) = 0.
-      return
-  10  twopi = 8.*atan(1.0)
-      tpn = twopi/float(n)
-c
-c compute x1=real(x(1)) and x2=imaginary(x(n/2-1))
-c x(n) = x(n)*4.*sin(twopi*(i-1)/n)
-c
-      t1 = 0.
-c
-c cosd and sind are multipliers for recursion for sin and cos
-c cosi and sini are initial conditions for recursion for sin and cos
-c
-      cosd = cos(tpn*2.)
-      sind = sin(tpn*2.)
-      cosi = 1.
-      sini = 0.
-      no2 = n/2
-      do 20 i=1,no2,2
-        t = x(i)*cosi
-        x(i) = x(i)*4.*sini
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-        t1 = t1 + t
-  20  continue
-c
-c reset initial conditions (cosi,sini) for new recursion
-c
-      cosi = cos(tpn)
-      sini = sin(tpn)
-      t2 = 0.
-      do 30 i=2,no2,2
-        t = x(i)*cosi
-        x(i) = x(i)*4.*sini
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-        t2 = t2 + t
-  30  continue
-      x1 = 2.*(t1+t2)
-      x2 = 2.*(t1-t2)
-c
-c take n/2 point (real) fft of preprocessed sequence x
-c
-      call fast(x, no2)
-c
-c for n = 4--skip recursion and initial conditions
-c
-      if (n.eq.4) go to 50
-c
-c initial conditions for recursion
-c
-      x(2) = -x(1)/2.
-      x(1) = x1
-c
-c for n = 8, skip recursion
-c
-      if (n.eq.8) go to 50
-c
-c unscramble y(k) using recursion formula
-c
-      nind = no2 - 2
-      do 40 i=3,nind,2
-        t = x(i)
-        x(i) = x(i-2) + x(i+1)
-        x(i+1) = x(i-1) - t
-  40  continue
-  50  x(no2) = x(no2+1)/2.
-      x(no2-1) = x2
-      return
-      end
diff --git a/math/ieee/chap1/fftsoh.f b/math/ieee/chap1/fftsoh.f
deleted file mode 100644
index afc18c179..000000000
--- a/math/ieee/chap1/fftsoh.f
+++ /dev/null
@@ -1,81 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: fftsoh
-c compute dft for real, symmetric, odd harmonic, n-point sequence
-c using n/4-point fft
-c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1
-c odd harmonic means x(2*k)=0, all k, where x(k) is the dft of x(m)
-c x(m) has the property x(m)=-x(n/2-m), m=0,1,...,n/4-1,  x(n/4)=0
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine fftsoh(x, n, y)
-      dimension x(1), y(1)
-c
-c  x = real array which on input contains the n/4 points of the
-c      input sequence (symmetrical)
-c      on output x contains  the n/4 real points of the odd harmonics
-c      of the transform of the input--i.e. the zero valued imaginary
-c      parts are not given nor are the zero-valued even harmonics
-c  n = true size of input
-c  y = scratch array of size n/4+2
-c
-c
-c handle n = 2 and n = 4 cases separately
-c
-      if (n.gt.4) go to 20
-      if (n.eq.4) go to 10
-c
-c for n=2, assume x(1)=x0, x(2)=-x0, compute dft directly
-c
-      x(1) = 2.*x(1)
-      return
-c
-c n = 4 case, compute dft directly
-c
-  10  x(1) = 2.*x(1)
-      return
-  20  twopi = 8.*atan(1.0)
-c
-c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1))
-c
-      no2 = n/2
-      no4 = n/4
-      no8 = n/8
-      if (no8.eq.1) go to 40
-      do 30 i=2,no8
-        ind = 2*i
-        t1 = x(ind) - x(ind-2)
-        y(i) = x(ind-1) + t1
-        ind1 = n/4 + 2 - i
-        y(ind1) = -x(ind-1) + t1
-  30  continue
-  40  y(1) = x(1)
-      y(no8+1) = -2.*x(no4)
-c
-c the sequence y (n/4 points) has only odd harmonics
-c call subroutine fftohm to exploit odd harmonics
-c
-      call fftohm(y, no2)
-c
-c form original dft from complex odd harmonics of y(k)
-c by unscrambling y(k)
-c
-      tpn = twopi/float(n)
-      cosi = 2.*cos(tpn)
-      sini = 2.*sin(tpn)
-      cosd = cos(tpn*2.)
-      sind = sin(tpn*2.)
-      do 50 i=1,no8
-        ind = 2*i
-        bk = y(ind)/sini
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-        ak = y(ind-1)
-        x(i) = ak + bk
-        ind1 = n/4 + 1 - i
-        x(ind1) = ak - bk
-  50  continue
-      return
-      end
diff --git a/math/ieee/chap1/fftsym.f b/math/ieee/chap1/fftsym.f
deleted file mode 100644
index b1e795bb6..000000000
--- a/math/ieee/chap1/fftsym.f
+++ /dev/null
@@ -1,84 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: fftsym
-c compute dft for real, symmetric, n-point sequence x(m) using
-c n/2-point fft
-c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine fftsym(x, n, y)
-      dimension x(1), y(1)
-c
-c x = real array which on input contains the n/2+1 points of the
-c     input sequence (symmetrical)
-c     on output x contains the n/2+1 real points of the transform of
-c     the input--i.e. the zero valued imaginary parts are not returned
-c n = true size of input
-c y = scratch array of size n/2+2
-c
-c
-c for n = 2, compute dft directly
-c
-      if (n.gt.2) go to 10
-      t = x(1) + x(2)
-      x(2) = x(1) - x(2)
-      x(1) = t
-      return
-  10  twopi = 8.*atan(1.0)
-c
-c first compute b0 term, where b0=sum of odd values of x(m)
-c
-      no2 = n/2
-      no4 = n/4
-      nind = no2 + 1
-      b0 = 0.
-      do 20 i=2,nind,2
-        b0 = b0 + x(i)
-  20  continue
-      b0 = b0*2.
-c
-c for n = 4 skip recursion loop
-c
-      if (n.eq.4) go to 40
-c
-c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1))
-c
-      do 30 i=2,no4
-        ind = 2*i
-        t1 = x(ind) - x(ind-2)
-        y(i) = x(ind-1) + t1
-        ind1 = no2 + 2 - i
-        y(ind1) = x(ind-1) - t1
-  30  continue
-  40  y(1) = x(1)
-      y(no4+1) = x(no2+1)
-c
-c take n/2 point (real) fft of y
-c
-      call fast(y, no2)
-c
-c form original dft by unscrambling y(k)
-c use recursion to give sin(tpn*i) multiplier
-c
-      tpn = twopi/float(n)
-      cosi = 2.*cos(tpn)
-      sini = 2.*sin(tpn)
-      cosd = cosi/2.
-      sind = sini/2.
-      nind = no4 + 1
-      do 50 i=2,nind
-        ind = 2*i
-        bk = y(ind)/sini
-        ak = y(ind-1)
-        x(i) = ak + bk
-        nind1 = n/2 + 2 - i
-        x(nind1) = ak - bk
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-  50  continue
-      x(1) = b0 + y(1)
-      x(no2+1) = y(1) - b0
-      return
-      end
diff --git a/math/ieee/chap1/ford1.f b/math/ieee/chap1/ford1.f
deleted file mode 100644
index 4982246f3..000000000
--- a/math/ieee/chap1/ford1.f
+++ /dev/null
@@ -1,24 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  ford1
-c in-place reordering subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine ford1(m, b)
-      dimension b(2)
-c
-      k = 4
-      kl = 2
-      n = 2**m
-      do 40 j=4,n,2
-        if (k-j) 20, 20, 10
-  10    t = b(j)
-        b(j) = b(k)
-        b(k) = t
-  20    k = k - 2
-        if (k-kl) 30, 30, 40
-  30    k = 2*j
-        kl = j
-  40  continue
-      return
-      end
diff --git a/math/ieee/chap1/ford2.f b/math/ieee/chap1/ford2.f
deleted file mode 100644
index 8ee26d697..000000000
--- a/math/ieee/chap1/ford2.f
+++ /dev/null
@@ -1,46 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  ford2
-c in-place reordering subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine ford2(m, b)
-      dimension l(15), b(2)
-      equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
-     *    (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
-     *    (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
-     *    (l1,l(15))
-      n = 2**m
-      l(1) = n
-      do 10 k=2,m
-        l(k) = l(k-1)/2
-  10  continue
-      do 20 k=m,14
-        l(k+1) = 2
-  20  continue
-      ij = 2
-      do 40 j1=2,l1,2
-      do 40 j2=j1,l2,l1
-      do 40 j3=j2,l3,l2
-      do 40 j4=j3,l4,l3
-      do 40 j5=j4,l5,l4
-      do 40 j6=j5,l6,l5
-      do 40 j7=j6,l7,l6
-      do 40 j8=j7,l8,l7
-      do 40 j9=j8,l9,l8
-      do 40 j10=j9,l10,l9
-      do 40 j11=j10,l11,l10
-      do 40 j12=j11,l12,l11
-      do 40 j13=j12,l13,l12
-      do 40 j14=j13,l14,l13
-      do 40 ji=j14,l15,l14
-        if (ij-ji) 30, 40, 40
-  30    t = b(ij-1)
-        b(ij-1) = b(ji-1)
-        b(ji-1) = t
-        t = b(ij)
-        b(ij) = b(ji)
-        b(ji) = t
-  40    ij = ij + 2
-      return
-      end
diff --git a/math/ieee/chap1/fourea.f b/math/ieee/chap1/fourea.f
deleted file mode 100644
index d412c0e26..000000000
--- a/math/ieee/chap1/fourea.f
+++ /dev/null
@@ -1,98 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: fourea
-c performs cooley-tukey fast fourier transform
-c-----------------------------------------------------------------------
-c
-      subroutine fourea(data, n, isi)
-c
-c the cooley-tukey fast fourier transform in ansi fortran
-c
-c data is a one-dimensional complex array whose length, n, is a
-c power of two.  isi is +1 for an inverse transform and -1 for a
-c forward transform.  transform values are returned in the input
-c array, replacing the input.
-c transform(j)=sum(data(i)*w**((i-1)*(j-1))), where i and j run
-c from 1 to n and w = exp (isi*2*pi*sqrt(-1)/n).  program also
-c computes inverse transform, for which the defining expression
-c is invtr(j)=(1/n)*sum(data(i)*w**((i-1)*(j-1))).
-c running time is proportional to n*log2(n), rather than to the
-c classical n**2.
-c after program by brenner, june 1967. this is a very short version
-c of the fft and is intended mainly for demonstration. programs
-c are available in this collection which run faster and are not
-c restricted to powers of 2 or to one-dimensional arrays.
-c see -- ieee trans audio (june 1967), special issue on fft.
-c
-      complex data(1)
-      complex temp, w
-      ioutd = i1mach(2)
-c
-c check for power of two up to 15
-c
-      nn = 1
-      do 10 i=1,15
-        m = i
-        nn = nn*2
-        if (nn.eq.n) go to 20
-  10  continue
-      write (ioutd,9999)
-9999  format (30h n not a power of 2 for fourea)
-      stop
-  20  continue
-c
-      pi = 4.*atan(1.)
-      fn = n
-c
-c this section puts data in bit-reversed order
-c
-      j = 1
-      do 80 i=1,n
-c
-c at this point, i and j are a bit reversed pair (except for the
-c displacement of +1)
-c
-        if (i-j) 30, 40, 40
-c
-c exchange data(i) with data(j) if i.lt.j.
-c
-  30    temp = data(j)
-        data(j) = data(i)
-        data(i) = temp
-c
-c implement j=j+1, bit-reversed counter
-c
-  40    m = n/2
-  50    if (j-m) 70, 70, 60
-  60    j = j - m
-        m = (m+1)/2
-        go to 50
-  70    j = j + m
-  80  continue
-c
-c now compute the butterflies
-c
-      mmax = 1
-  90  if (mmax-n) 100, 130, 130
- 100  istep = 2*mmax
-      do 120 m=1,mmax
-        theta = pi*float(isi*(m-1))/float(mmax)
-        w = cmplx(cos(theta),sin(theta))
-        do 110 i=m,n,istep
-          j = i + mmax
-          temp = w*data(j)
-          data(j) = data(i) - temp
-          data(i) = data(i) + temp
- 110    continue
- 120  continue
-      mmax = istep
-      go to 90
- 130  if (isi) 160, 140, 140
-c
-c for inv trans -- isi=1 -- multiply output by 1/n
-c
- 140  do 150 i=1,n
-        data(i) = data(i)/fn
- 150  continue
- 160  return
-      end
diff --git a/math/ieee/chap1/fr2tr.f b/math/ieee/chap1/fr2tr.f
deleted file mode 100644
index 24a103ffc..000000000
--- a/math/ieee/chap1/fr2tr.f
+++ /dev/null
@@ -1,15 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  fr2tr
-c radix 2 iteration subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine fr2tr(int, b0, b1)
-      dimension b0(2), b1(2)
-      do 10 k=1,int
-        t = b0(k) + b1(k)
-        b1(k) = b0(k) - b1(k)
-        b0(k) = t
-  10  continue
-      return
-      end
diff --git a/math/ieee/chap1/fr4syn.f b/math/ieee/chap1/fr4syn.f
deleted file mode 100644
index 5ce978f00..000000000
--- a/math/ieee/chap1/fr4syn.f
+++ /dev/null
@@ -1,109 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  fr4syn
-c radix 4 synthesis
-c-----------------------------------------------------------------------
-c
-c
-      subroutine fr4syn(int, nn, b0, b1, b2, b3, b4, b5, b6, b7)
-      dimension l(15), b0(2), b1(2), b2(2), b3(2), b4(2), b5(2), b6(2),
-     *    b7(2)
-      common /const/ pii, p7, p7two, c22, s22, pi2
-      equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
-     *    (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
-     *    (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
-     *    (l1,l(15))
-c
-      l(1) = nn/4
-      do 40 k=2,15
-        if (l(k-1)-2) 10, 20, 30
-  10    l(k-1) = 2
-  20    l(k) = 2
-        go to 40
-  30    l(k) = l(k-1)/2
-  40  continue
-c
-      piovn = pii/float(nn)
-      ji = 3
-      jl = 2
-      jr = 2
-c
-      do 120 j1=2,l1,2
-      do 120 j2=j1,l2,l1
-      do 120 j3=j2,l3,l2
-      do 120 j4=j3,l4,l3
-      do 120 j5=j4,l5,l4
-      do 120 j6=j5,l6,l5
-      do 120 j7=j6,l7,l6
-      do 120 j8=j7,l8,l7
-      do 120 j9=j8,l9,l8
-      do 120 j10=j9,l10,l9
-      do 120 j11=j10,l11,l10
-      do 120 j12=j11,l12,l11
-      do 120 j13=j12,l13,l12
-      do 120 j14=j13,l14,l13
-      do 120 jthet=j14,l15,l14
-        th2 = jthet - 2
-        if (th2) 50, 50, 90
-  50    do 60 k=1,int
-          t0 = b0(k) + b1(k)
-          t1 = b0(k) - b1(k)
-          t2 = b2(k)*2.0
-          t3 = b3(k)*2.0
-          b0(k) = t0 + t2
-          b2(k) = t0 - t2
-          b1(k) = t1 + t3
-          b3(k) = t1 - t3
-  60    continue
-c
-        if (nn-4) 120, 120, 70
-  70    k0 = int*4 + 1
-        kl = k0 + int - 1
-        do 80 k=k0,kl
-          t2 = b0(k) - b2(k)
-          t3 = b1(k) + b3(k)
-          b0(k) = (b0(k)+b2(k))*2.0
-          b2(k) = (b3(k)-b1(k))*2.0
-          b1(k) = (t2+t3)*p7two
-          b3(k) = (t3-t2)*p7two
-  80    continue
-        go to 120
-  90    arg = th2*piovn
-        c1 = cos(arg)
-        s1 = -sin(arg)
-        c2 = c1**2 - s1**2
-        s2 = c1*s1 + c1*s1
-        c3 = c1*c2 - s1*s2
-        s3 = c2*s1 + s2*c1
-c
-        int4 = int*4
-        j0 = jr*int4 + 1
-        k0 = ji*int4 + 1
-        jlast = j0 + int - 1
-        do 100 j=j0,jlast
-          k = k0 + j - j0
-          t0 = b0(j) + b6(k)
-          t1 = b7(k) - b1(j)
-          t2 = b0(j) - b6(k)
-          t3 = b7(k) + b1(j)
-          t4 = b2(j) + b4(k)
-          t5 = b5(k) - b3(j)
-          t6 = b5(k) + b3(j)
-          t7 = b4(k) - b2(j)
-          b0(j) = t0 + t4
-          b4(k) = t1 + t5
-          b1(j) = (t2+t6)*c1 - (t3+t7)*s1
-          b5(k) = (t2+t6)*s1 + (t3+t7)*c1
-          b2(j) = (t0-t4)*c2 - (t1-t5)*s2
-          b6(k) = (t0-t4)*s2 + (t1-t5)*c2
-          b3(j) = (t2-t6)*c3 - (t3-t7)*s3
-          b7(k) = (t2-t6)*s3 + (t3-t7)*c3
- 100    continue
-        jr = jr + 2
-        ji = ji - 2
-        if (ji-jl) 110, 110, 120
- 110    ji = 2*jr - 1
-        jl = jr
- 120  continue
-      return
-      end
diff --git a/math/ieee/chap1/fr4tr.f b/math/ieee/chap1/fr4tr.f
deleted file mode 100644
index 16ba3fe6a..000000000
--- a/math/ieee/chap1/fr4tr.f
+++ /dev/null
@@ -1,118 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  fr4tr
-c radix 4 iteration subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine fr4tr(int, nn, b0, b1, b2, b3, b4, b5, b6, b7)
-      dimension l(15), b0(2), b1(2), b2(2), b3(2), b4(2), b5(2), b6(2),
-     *    b7(2)
-      common /cons/ pii, p7, p7two, c22, s22, pi2
-      equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
-     *    (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
-     *    (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
-     *    (l1,l(15))
-c
-c jthet is a reversed binary counter, jr steps two at a time to
-c locate the real parts of intermediate results, and ji locates
-c the imaginary part corresponding to jr.
-c
-      l(1) = nn/4
-      do 40 k=2,15
-        if (l(k-1)-2) 10, 20, 30
-  10    l(k-1) = 2
-  20    l(k) = 2
-        go to 40
-  30    l(k) = l(k-1)/2
-  40  continue
-c
-      piovn = pii/float(nn)
-      ji = 3
-      jl = 2
-      jr = 2
-c
-      do 120 j1=2,l1,2
-      do 120 j2=j1,l2,l1
-      do 120 j3=j2,l3,l2
-      do 120 j4=j3,l4,l3
-      do 120 j5=j4,l5,l4
-      do 120 j6=j5,l6,l5
-      do 120 j7=j6,l7,l6
-      do 120 j8=j7,l8,l7
-      do 120 j9=j8,l9,l8
-      do 120 j10=j9,l10,l9
-      do 120 j11=j10,l11,l10
-      do 120 j12=j11,l12,l11
-      do 120 j13=j12,l13,l12
-      do 120 j14=j13,l14,l13
-      do 120 jthet=j14,l15,l14
-        th2 = jthet - 2
-        if (th2) 50, 50, 90
-  50    do 60 k=1,int
-          t0 = b0(k) + b2(k)
-          t1 = b1(k) + b3(k)
-          b2(k) = b0(k) - b2(k)
-          b3(k) = b1(k) - b3(k)
-          b0(k) = t0 + t1
-          b1(k) = t0 - t1
-  60    continue
-c
-        if (nn-4) 120, 120, 70
-  70    k0 = int*4 + 1
-        kl = k0 + int - 1
-        do 80 k=k0,kl
-          pr = p7*(b1(k)-b3(k))
-          pi = p7*(b1(k)+b3(k))
-          b3(k) = b2(k) + pi
-          b1(k) = pi - b2(k)
-          b2(k) = b0(k) - pr
-          b0(k) = b0(k) + pr
-  80    continue
-        go to 120
-c
-  90    arg = th2*piovn
-        c1 = cos(arg)
-        s1 = sin(arg)
-        c2 = c1**2 - s1**2
-        s2 = c1*s1 + c1*s1
-        c3 = c1*c2 - s1*s2
-        s3 = c2*s1 + s2*c1
-c
-        int4 = int*4
-        j0 = jr*int4 + 1
-        k0 = ji*int4 + 1
-        jlast = j0 + int - 1
-        do 100 j=j0,jlast
-          k = k0 + j - j0
-          r1 = b1(j)*c1 - b5(k)*s1
-          r5 = b1(j)*s1 + b5(k)*c1
-          t2 = b2(j)*c2 - b6(k)*s2
-          t6 = b2(j)*s2 + b6(k)*c2
-          t3 = b3(j)*c3 - b7(k)*s3
-          t7 = b3(j)*s3 + b7(k)*c3
-          t0 = b0(j) + t2
-          t4 = b4(k) + t6
-          t2 = b0(j) - t2
-          t6 = b4(k) - t6
-          t1 = r1 + t3
-          t5 = r5 + t7
-          t3 = r1 - t3
-          t7 = r5 - t7
-          b0(j) = t0 + t1
-          b7(k) = t4 + t5
-          b6(k) = t0 - t1
-          b1(j) = t5 - t4
-          b2(j) = t2 - t7
-          b5(k) = t6 + t3
-          b4(k) = t2 + t7
-          b3(j) = t3 - t6
- 100    continue
-c
-        jr = jr + 2
-        ji = ji - 2
-        if (ji-jl) 110, 110, 120
- 110    ji = 2*jr - 1
-        jl = jr
- 120  continue
-      return
-      end
diff --git a/math/ieee/chap1/fsst.f b/math/ieee/chap1/fsst.f
deleted file mode 100644
index 75ce56cbc..000000000
--- a/math/ieee/chap1/fsst.f
+++ /dev/null
@@ -1,71 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  fsst
-c fourier synthesis subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine fsst(b, n)
-c
-c this subroutine synthesizes the real vector b(k), for
-c k=1,2,...,n, from the fourier coefficients stored in the
-c b array of size n+2.  the dc term is in b(1) with b(2) equal
-c to  0.  the jth harmonic is stored as b(2*j+1) + i b(2*j+2).
-c the n/2 harmonic is in b(n+1) with b(n+2) equal to 0.
-c the subroutine is called as fsst(b,n) where n=2**m and
-c b is the real array discussed above.
-c
-      dimension b(2)
-      common /const/ pii, p7, p7two, c22, s22, pi2
-c
-c iw is a machine dependent write device number
-c
-      iw = i1mach(2)
-c
-      pii = 4.*atan(1.)
-      pi8 = pii/8.
-      p7 = 1./sqrt(2.)
-      p7two = 2.*p7
-      c22 = cos(pi8)
-      s22 = sin(pi8)
-      pi2 = 2.*pii
-      do 10 i=1,15
-        m = i
-        nt = 2**i
-        if (n.eq.nt) go to 20
-  10  continue
-      write (iw,9999)
-9999  format (33h n is not a power of two for fsst)
-      stop
-  20  b(2) = b(n+1)
-      do 30 i=4,n,2
-        b(i) = -b(i)
-  30  continue
-c
-c scale the input by n
-c
-      do 40 i=1,n
-        b(i) = b(i)/float(n)
-  40  continue
-      n4pow = m/2
-c
-c scramble the inputs
-c
-      call ford2(m, b)
-      call ford1(m, b)
-c
-      if (n4pow.eq.0) go to 60
-      nn = 4*n
-      do 50 it=1,n4pow
-        nn = nn/4
-        int = n/nn
-        call fr4syn(int, nn, b(1), b(int+1), b(2*int+1), b(3*int+1),
-     *      b(1), b(int+1), b(2*int+1), b(3*int+1))
-  50  continue
-c
-c do a radix 2 iteration if one is required
-c
-  60  if (m-n4pow*2) 80, 80, 70
-  70  int = n/2
-      call fr2tr(int, b(1), b(int+1))
-  80  return
-      end
diff --git a/math/ieee/chap1/iftaoh.f b/math/ieee/chap1/iftaoh.f
deleted file mode 100644
index 8cb0a9083..000000000
--- a/math/ieee/chap1/iftaoh.f
+++ /dev/null
@@ -1,87 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: iftaoh
-c compute idft for real, antisymmetric, odd harmonic, n-point sequence
-c using n/4-point fft
-c antisymmetric sequence means x(m)=-x(n-m), m=1,...,n/2-1
-c odd harmonic means x(2*k)=0, all k, where x(k) is the dft of x(m)
-c x(m)has the property x(m)=x(n/2-m), m=0,1,...,n/4-1,  x(0)=0
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine iftaoh(x, n, y)
-      dimension x(1), y(1)
-c
-c  x = real array which on input contains the n/4 imaginary points
-c      of the odd harmonics of the transform of the original time
-c      sequence--i.e. the zero valued real parts are not input nor
-c      are the zero-valued even harmonics
-c      on output x contains the first (n/4+1) points of the original
-c      time sequence (antisymmetrical)
-c  n = true size of input
-c  y = scratch array of size n/4+2
-c
-c
-c handle n = 2 and n = 4 cases separately
-c
-      if (n.gt.4) go to 20
-      if (n.eq.4) go to 10
-c
-c for n=2  assume x(1)=0, x(2)=0, compute idft directly
-c
-      x(1) = 0.
-      return
-c
-c for n=4, assume x(1)=x(3)=0, x(2)=-x(4)=x0, compute idft directly
-c
-  10  x(2) = -x(1)/2.
-      x(1) = 0.
-      return
-c
-c code for values of n which are multiples of 8
-c
-  20  twopi = 8.*atan(1.0)
-      no2 = n/2
-      no4 = n/4
-      no8 = n/8
-      tpn = twopi/float(n)
-c
-c scramble original dft (x(k)) to give y(k)
-c use recursion to give sin multipliers
-c
-      cosi = cos(tpn)
-      sini = sin(tpn)
-      cosd = cos(tpn*2.)
-      sind = sin(tpn*2.)
-      do 30 i=1,no8
-        ind = 2*i
-        ind1 = no4 + 1 - i
-        ak = (x(i)-x(ind1))/2.
-        bk = -(x(i)+x(ind1))
-        y(ind) = ak
-        y(ind-1) = bk*sini
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-  30  continue
-c
-c the sequence y(k) is an odd harmonic sequence
-c use subroutine iftohm to give y(m)
-c
-      call iftohm(y, no2)
-c
-c form x sequence from y sequence
-c
-      x(2) = y(1)/2.
-      x(1) = 0.
-      if (n.eq.8) return
-      do 40 i=2,no8
-        ind = 2*i
-        ind1 = no4 + 2 - i
-        x(ind-1) = (y(i)+y(ind1))/2.
-        t1 = (y(i)-y(ind1))/2.
-        x(ind) = t1 + x(ind-2)
-  40  continue
-      x(no4+1) = y(no8+1)
-      return
-      end
diff --git a/math/ieee/chap1/iftasm.f b/math/ieee/chap1/iftasm.f
deleted file mode 100644
index 6ec12c7d7..000000000
--- a/math/ieee/chap1/iftasm.f
+++ /dev/null
@@ -1,77 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: iftasm
-c compute idft for real, antisymmetric, n-point sequence x(m) using
-c n/2-point fft
-c antisymmetric sequence means x(m)=-x(n-m), m=1,...,n/2-1
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine iftasm(x, n, y)
-      dimension x(1), y(1)
-c
-c x = imaginary array which on input contains the n/2+1 real points of
-c     the transform of the input--i.e. the zero valued real parts
-c     are not given as input
-c     on output x contains the n/2 points of the time sequence
-c     (antisymmetrical)
-c n = true size of input
-c y = scratch array of size n/2+2
-c
-c
-c for n = 2, assume x(1)=0, x(2)=0
-c
-      if (n.gt.2) go to 10
-      x(1) = 0
-      x(2) = 0
-      return
-  10  twopi = 8.*atan(1.0)
-c
-c first compute x1=x(1) term directly
-c use recursion on the sine cosine terms
-c
-      no2 = n/2
-      no4 = n/4
-      tpn = twopi/float(n)
-c
-c scramble original dft (x(k)) to give y(k)
-c use recursion relation to give sin(tpn*i) multiplier
-c
-      cosi = cos(tpn)
-      sini = sin(tpn)
-      cosd = cosi
-      sind = sini
-      nind = no4 + 1
-      do 20 i=2,nind
-        ind = 2*i
-        ind1 = no2 + 2 - i
-        ak = (x(i)-x(ind1))/2.
-        bk = -(x(i)+x(ind1))
-        y(ind) = ak
-        y(ind-1) = bk*sini
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-  20  continue
-      y(1) = 0.
-      y(2) = 0.
-c
-c take n/2 point idft of y
-c
-      call fsst(y, no2)
-c
-c form x sequence from y sequence
-c
-      x(2) = y(1)/2.
-      x(1) = 0.
-      if (n.eq.4) go to 40
-      do 30 i=2,no4
-        ind = 2*i
-        ind1 = no2 + 2 - i
-        x(ind-1) = (y(i)-y(ind1))/2.
-        t1 = (y(i)+y(ind1))/2.
-        x(ind) = t1 + x(ind-2)
-  30  continue
-  40  x(no2) = -y(no4+1)/2.
-      return
-      end
diff --git a/math/ieee/chap1/iftohm.f b/math/ieee/chap1/iftohm.f
deleted file mode 100644
index df084e870..000000000
--- a/math/ieee/chap1/iftohm.f
+++ /dev/null
@@ -1,83 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: iftohm
-c compute idft for real, n-point, odd harmonic sequences using an
-c n/2 point fft
-c odd harmonic means x(2*k)=0, all k where x(k) is the dft of x(m)
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine iftohm(x, n)
-      dimension x(1)
-c
-c x = real array which on input contains the n/4 complex values of the
-c     odd harmonics of the input--stored in the sequence re(x(1)),
-c     im(x(1)),re(x(2)),im(x(2)),...
-c     on output x contains the first n/2 points of the input
-c ****note: x must be dimensioned to size n/2+2 for fft routine
-c n = true size of x sequence
-c
-c first compute real(x(1)) and real(x(n/2-1)) separately
-c also simultaneously multiply original sequence by sin(twopi*(m-1)/n)
-c sin and cos are computed recursively
-c
-c
-c for n = 2, assume x(1)=x0, x(2)=-x0, compute idft directly
-c
-      if (n.gt.2) go to 10
-      x(1) = 0.5*x(1)
-      x(2) = -x(1)
-      return
-  10  twopi = 8.*atan(1.0)
-      tpn = twopi/float(n)
-      no2 = n/2
-      no4 = n/4
-      nind = no2
-c
-c solve for x(0)=x0 directly
-c
-      x0 = 0.
-      do 20 i=1,no2,2
-        x0 = x0 + 2.*x(i)
-  20  continue
-      x0 = x0/float(n)
-c
-c form y(k)=j*(x(2k+1)-x(2k-1))
-c overwrite x array with y sequence
-c
-      xpr = x(1)
-      xpi = x(2)
-      x(1) = -2.*x(2)
-      x(2) = 0.
-      if (no4.eq.1) go to 40
-      do 30 i=3,nind,2
-        ti = x(i) - xpr
-        tr = -x(i+1) + xpi
-        xpr = x(i)
-        xpi = x(i+1)
-        x(i) = tr
-        x(i+1) = ti
-  30  continue
-  40  x(no2+1) = 2.*xpi
-      x(no2+2) = 0.
-c
-c take n/2 point (real) ifft of preprocessed sequence x
-c
-      call fsst(x, no2)
-c
-c solve for x(m) by dividing by 4*sin(twopi*m/n) for m=1,2,...,n/2-1
-c for m=0 substitute precomputed value x0
-c
-      cosi = 4.
-      sini = 0.
-      cosd = cos(tpn)
-      sind = sin(tpn)
-      do 50 i=2,no2
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-        x(i) = x(i)/sini
-  50  continue
-      x(1) = x0
-      return
-      end
diff --git a/math/ieee/chap1/iftsoh.f b/math/ieee/chap1/iftsoh.f
deleted file mode 100644
index 6098e6e4d..000000000
--- a/math/ieee/chap1/iftsoh.f
+++ /dev/null
@@ -1,94 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: iftsoh
-c compute idft for real, symmetric, odd harmonic, n-point sequence
-c using n/4-point fft
-c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1
-c odd harmonic means x(2*k)=0, all k, where x(k) is the dft of x(m)
-c x(m) has the property x(m)=-x(n/2-m), m=0,1,...,n/4-1,  x(n/4)=0
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine iftsoh(x, n, y)
-      dimension x(1), y(1)
-c
-c  x = real array which on input contains the n/4 real points of
-c      the odd harmonics of the transform of the original time sequence
-c      i.e. the zero valued imaginary parts are not given nor are the
-c      zero valued even harmonics
-c      on output x contains the first n/4 points of the original input
-c      sequence (symmetrical)
-c  n = true size of input
-c  y = scratch array of size n/4+2
-c
-c
-c handle n = 2 and n = 4 cases separately
-c
-      if (n.gt.4) go to 10
-c
-c for n=2, 4 assume x(1)=x0, x(2)=-x0, compute idft directly
-c
-      x(1) = x(1)/2.
-      return
-c
-c code for values of n which are multiples of 8
-c
-  10  twopi = 8.*atan(1.0)
-      no2 = n/2
-      no4 = n/4
-      no8 = n/8
-      tpn = twopi/float(n)
-c
-c first compute x1=x(1) term directly
-c use recursion on the sine cosine terms
-c
-      cosd = cos(tpn*2.)
-      sind = sin(tpn*2.)
-      cosi = 2.*cos(tpn)
-      sini = 2.*sin(tpn)
-      x1 = 0.
-      do 20 i=1,no4
-        x1 = x1 + x(i)*cosi
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-  20  continue
-      x1 = x1/float(n)
-c
-c scramble original dft (x(k)) to give y(k)
-c use recursion relation to give sin multipliers
-c
-      cosi = cos(tpn)
-      sini = sin(tpn)
-      do 30 i=1,no8
-        ind = 2*i
-        ind1 = no4 + 1 - i
-        ak = (x(i)+x(ind1))/2.
-        bk = (x(i)-x(ind1))
-        y(ind-1) = ak
-        y(ind) = bk*sini
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-  30  continue
-c
-c the sequence y(k) is the odd harmonics dft output
-c use subroutine iftohm to obtain y(m), the inverse transform
-c
-      call iftohm(y, no2)
-c
-c form x(m) sequence from y(m) sequence
-c use x1 initial condition on the recursion
-c
-      x(1) = y(1)
-      x(2) = x1
-      if (no8.eq.1) return
-      do 40 i=2,no8
-        ind = 2*i
-        ind1 = no4 + 2 - i
-        t1 = (y(i)+y(ind1))/2.
-        x(ind-1) = (y(i)-y(ind1))/2.
-        x(ind) = t1 + x(ind-2)
-  40  continue
-      return
-      end
diff --git a/math/ieee/chap1/iftsym.f b/math/ieee/chap1/iftsym.f
deleted file mode 100644
index e45c9a37c..000000000
--- a/math/ieee/chap1/iftsym.f
+++ /dev/null
@@ -1,90 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: iftsym
-c compute idft for real, symmetric, n-point sequence x(m) using
-c n/2-point fft
-c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1
-c note: index m is sequence index--not fortran index
-c-----------------------------------------------------------------------
-c
-      subroutine iftsym(x, n, y)
-      dimension x(1), y(1)
-c
-c x = real array which on input contains the n/2+1 real points of the
-c     transform of the input--i.e. the zero valued imaginary parts
-c     are not given as input
-c     on output x contains the n/2+1 points of the time sequence
-c     (symmetrical)
-c n = true size of input
-c y = scratch array of size n/2+2
-c
-c
-c for n = 2, compute idft directly
-c
-      if (n.gt.2) go to 10
-      t = (x(1)+x(2))/2.
-      x(2) = (x(1)-x(2))/2.
-      x(1) = t
-      return
-  10  twopi = 8.*atan(1.0)
-c
-c first compute x1=x(1) term directly
-c use recursion on the sine cosine terms
-c
-      no2 = n/2
-      no4 = n/4
-      tpn = twopi/float(n)
-      cosd = cos(tpn)
-      sind = sin(tpn)
-      cosi = 2.
-      sini = 0.
-      x1 = x(1) - x(no2+1)
-      do 20 i=2,no2
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-        x1 = x1 + x(i)*cosi
-  20  continue
-      x1 = x1/float(n)
-c
-c scramble original dft (x(k)) to give y(k)
-c use recursion relation to generate sin(tpn*i) multiplier
-c
-      cosi = cos(tpn)
-      sini = sin(tpn)
-      cosd = cosi
-      sind = sini
-      y(1) = (x(1)+x(no2+1))/2.
-      y(2) = 0.
-      nind = no4 + 1
-      do 30 i=2,nind
-        ind = 2*i
-        nind1 = no2 + 2 - i
-        ak = (x(i)+x(nind1))/2.
-        bk = (x(i)-x(nind1))
-        y(ind-1) = ak
-        y(ind) = bk*sini
-        temp = cosi*cosd - sini*sind
-        sini = cosi*sind + sini*cosd
-        cosi = temp
-  30  continue
-c
-c take n/2 point idft of y
-c
-      call fsst(y, no2)
-c
-c form x sequence from y sequence
-c
-      x(1) = y(1)
-      x(2) = x1
-      if (n.eq.4) go to 50
-      do 40 i=2,no4
-        ind = 2*i
-        ind1 = no2 + 2 - i
-        x(ind-1) = (y(i)+y(ind1))/2.
-        t1 = (y(i)-y(ind1))/2.
-        x(ind) = t1 + x(ind-2)
-  40  continue
-  50  x(no2+1) = y(no4+1)
-      return
-      end
diff --git a/math/ieee/chap1/inishl.f b/math/ieee/chap1/inishl.f
deleted file mode 100644
index 925cd657e..000000000
--- a/math/ieee/chap1/inishl.f
+++ /dev/null
@@ -1,179 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: inishl
-c   this subroutine initializes the wfta routine for a given
-c   value of the transform length n.  the factors of n are
-c   determined, the multiplication coefficients are calculated
-c   and stored in the array coef(.), the input and output
-c   permutation vectors are computed and stored in the arrays
-c   indx1(.) and indx2(.)
-c
-c-----------------------------------------------------------------------
-c
-      subroutine inishl(n,coef,xr,xi,indx1,indx2,ierr)
-      dimension coef(1),xr(1),xi(1)
-      integer s1,s2,s3,s4,indx1(1),indx2(1),p1
-      dimension co3(3),co4(4),co8(8),co9(11),co16(18),cda(18),cdb(11),
-     1cdc(9),cdd(6)
-      common na,nb,nc,nd,nd1,nd2,nd3,nd4
-c
-c   data statements assign short dft coefficients.
-c
-      data co4(1),co4(2),co4(3),co4(4)/4*1.0/
-c
-      data cda(1),cda(2),cda(3),cda(4),cda(5),cda(6),cda(7),
-     1 cda(8),cda(9),cda(10),cda(11),cda(12),cda(13),cda(14),
-     2 cda(15),cda(16),cda(17),cda(18)/18*1.0/
-c
-      data cdb(1),cdb(2),cdb(3),cdb(4),cdb(5),cdb(6),cdb(7),cdb(8),
-     1 cdb(9),cdb(10),cdb(11)/11*1.0/
-c
-      data ionce/1/
-c
-c   get multiplier constants
-c
-      if(ionce.ne.1) go to 20
-      call const(co3,co8,co16,co9,cdc,cdd)
-20    ionce=-1
-c
-c   following segment determines factors of n and chooses
-c   the appropriate short dft coefficients.
-c
-      iout=i1mach(2)
-      ierr=0
-      nd1=1
-      na=1
-      nb=1
-      nd2=1
-      nc=1
-      nd3=1
-      nd=1
-      nd4=1
-      if(n.le.0 .or. n.gt.5040) go to 190
-      if(16*(n/16).eq.n) go to 30
-      if(8*(n/8).eq.n) go to 40
-      if(4*(n/4).eq.n) go to 50
-      if(2*(n/2).ne.n) go to 70
-      nd1=2
-      na=2
-      cda(2)=1.0
-      go to 70
-30    nd1=18
-      na=16
-      do 31 j=1,18
-31    cda(j)=co16(j)
-      go to 70
-40    nd1=8
-      na=8
-      do 41 j=1,8
-41    cda(j)=co8(j)
-      go to 70
-50    nd1=4
-      na=4
-      do 51 j=1,4
-51    cda(j)=co4(j)
-70    if(3*(n/3).ne.n) go to 120
-      if(9*(n/9).eq.n) go to 100
-      nd2=3
-      nb=3
-      do 71 j=1,3
-71    cdb(j)=co3(j)
-      go to 120
-100   nd2=11
-      nb=9
-      do 110 j=1,11
-110   cdb(j)=co9(j)
-120   if(7*(n/7).ne.n) go to 160
-      nd3=9
-      nc=7
-160   if(5*(n/5).ne.n) go to 190
-      nd4=6
-      nd=5
-190   m=na*nb*nc*nd
-      if(m.eq.n) go to 250
-      write(iout,210)
-210   format(21h this n does not work)
-      ierr=-1
-      return
-c
-c   next segment generates the dft coefficients by
-c   multiplying together the short dft coefficients
-c
-250   j=1
-      do 300 n4=1,nd4
-      do 300 n3=1,nd3
-      do 300 n2=1,nd2
-      do 300 n1=1,nd1
-      coef(j)=cda(n1)*cdb(n2)*cdc(n3)*cdd(n4)
-      j=j+1
-300   continue
-c
-c   following segment forms the input indexing vector
-c
-      j=1
-      nu=nb*nc*nd
-      nv=na*nc*nd
-      nw=na*nb*nd
-      ny=na*nb*nc
-      k=1
-      do 440 n4=1,nd
-      do 430 n3=1,nc
-      do 420 n2=1,nb
-      do 410 n1=1,na
-405   if(k.le.n) go to 408
-      k=k-n
-      go to 405
-408   indx1(j)=k
-      j=j+1
-410   k=k+nu
-420   k=k+nv
-430   k=k+nw
-440   k=k+ny
-c
-c   following segment forms the output indexing vector
-c
-      m=1
-      s1=0
-      s2=0
-      s3=0
-      s4=0
-      if(na.eq.1) go to 530
-520   p1=m*nu-1
-      if((p1/na)*na.eq.p1) go to 510
-      m=m+1
-      go to 520
-510   s1=p1+1
-530   if(nb.eq.1) go to 540
-      m=1
-550   p1=m*nv-1
-      if((p1/nb)*nb.eq.p1) go to 560
-      m=m+1
-      go to 550
-560   s2=p1+1
-540   if(nc.eq.1) go to 630
-      m=1
-620   p1=m*nw-1
-      if((p1/nc)*nc.eq.p1) go to 610
-      m=m+1
-      go to 620
-610   s3=p1+1
-630   if(nd.eq.1) go to 660
-      m=1
-640   p1=m*ny-1
-      if((p1/nd)*nd.eq.p1) go to 650
-      m=m+1
-      go to 640
-650   s4=p1+1
-660   j=1
-      do 810 n4=1,nd
-      do 810 n3=1,nc
-      do 810 n2=1,nb
-      do 810 n1=1,na
-      indx2(j)=s1*(n1-1)+s2*(n2-1)+s3*(n3-1)+s4*(n4-1)+1
-900   if(indx2(j).le.n) go to 910
-      indx2(j)=indx2(j)-n
-      go to 900
-910   j=j+1
-810   continue
-      return
-      end
diff --git a/math/ieee/chap1/ord1.f b/math/ieee/chap1/ord1.f
deleted file mode 100644
index 21e4494e3..000000000
--- a/math/ieee/chap1/ord1.f
+++ /dev/null
@@ -1,24 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  ord1
-c in-place reordering subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine ord1(m, b)
-      dimension b(2)
-c
-      k = 4
-      kl = 2
-      n = 2**m
-      do 40 j=4,n,2
-        if (k-j) 20, 20, 10
-  10    t = b(j)
-        b(j) = b(k)
-        b(k) = t
-  20    k = k - 2
-        if (k-kl) 30, 30, 40
-  30    k = 2*j
-        kl = j
-  40  continue
-      return
-      end
diff --git a/math/ieee/chap1/ord2.f b/math/ieee/chap1/ord2.f
deleted file mode 100644
index 10f03ec7a..000000000
--- a/math/ieee/chap1/ord2.f
+++ /dev/null
@@ -1,46 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  ord2
-c in-place reordering subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine ord2(m, b)
-      dimension l(15), b(2)
-      equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
-     *    (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
-     *    (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
-     *    (l1,l(15))
-      n = 2**m
-      l(1) = n
-      do 10 k=2,m
-        l(k) = l(k-1)/2
-  10  continue
-      do 20 k=m,14
-        l(k+1) = 2
-  20  continue
-      ij = 2
-      do 40 j1=2,l1,2
-      do 40 j2=j1,l2,l1
-      do 40 j3=j2,l3,l2
-      do 40 j4=j3,l4,l3
-      do 40 j5=j4,l5,l4
-      do 40 j6=j5,l6,l5
-      do 40 j7=j6,l7,l6
-      do 40 j8=j7,l8,l7
-      do 40 j9=j8,l9,l8
-      do 40 j10=j9,l10,l9
-      do 40 j11=j10,l11,l10
-      do 40 j12=j11,l12,l11
-      do 40 j13=j12,l13,l12
-      do 40 j14=j13,l14,l13
-      do 40 ji=j14,l15,l14
-        if (ij-ji) 30, 40, 40
-  30    t = b(ij-1)
-        b(ij-1) = b(ji-1)
-        b(ji-1) = t
-        t = b(ij)
-        b(ij) = b(ji)
-        b(ji) = t
-  40    ij = ij + 2
-      return
-      end
diff --git a/math/ieee/chap1/r2tr.f b/math/ieee/chap1/r2tr.f
deleted file mode 100644
index 510763fe3..000000000
--- a/math/ieee/chap1/r2tr.f
+++ /dev/null
@@ -1,16 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  r2tr
-c radix 2 iteration subroutine
-c-----------------------------------------------------------------------
-c
-c
-      subroutine r2tr(int, b0, b1)
-      dimension b0(2), b1(2)
-      do 10 k=1,int
-        t = b0(k) + b1(k)
-        b1(k) = b0(k) - b1(k)
-        b0(k) = t
-  10  continue
-      return
-      end
diff --git a/math/ieee/chap1/r2tx.f b/math/ieee/chap1/r2tx.f
deleted file mode 100644
index 425c10f8d..000000000
--- a/math/ieee/chap1/r2tx.f
+++ /dev/null
@@ -1,18 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  r2tx
-c radix 2 iteration subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine r2tx(nthpo, cr0, cr1, ci0, ci1)
-      dimension cr0(2), cr1(2), ci0(2), ci1(2)
-      do 10 k=1,nthpo,2
-        r1 = cr0(k) + cr1(k)
-        cr1(k) = cr0(k) - cr1(k)
-        cr0(k) = r1
-        fi1 = ci0(k) + ci1(k)
-        ci1(k) = ci0(k) - ci1(k)
-        ci0(k) = fi1
-  10  continue
-      return
-      end
diff --git a/math/ieee/chap1/r4syn.f b/math/ieee/chap1/r4syn.f
deleted file mode 100644
index 77c97bba1..000000000
--- a/math/ieee/chap1/r4syn.f
+++ /dev/null
@@ -1,20 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  r4syn
-c radix 4 synthesis
-c-----------------------------------------------------------------------
-c
-      subroutine r4syn(int, b0, b1, b2, b3)
-      dimension b0(2), b1(2), b2(2), b3(2)
-      do 10 k=1,int
-        t0 = b0(k) + b1(k)
-        t1 = b0(k) - b1(k)
-        t2 = b2(k) + b2(k)
-        t3 = b3(k) + b3(k)
-        b0(k) = t0 + t2
-        b2(k) = t0 - t2
-        b1(k) = t1 + t3
-        b3(k) = t1 - t3
-  10  continue
-      return
-      end
diff --git a/math/ieee/chap1/r4tr.f b/math/ieee/chap1/r4tr.f
deleted file mode 100644
index 5fc46eab8..000000000
--- a/math/ieee/chap1/r4tr.f
+++ /dev/null
@@ -1,18 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  r4tr
-c radix 4 iteration subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine r4tr(int, b0, b1, b2, b3)
-      dimension b0(2), b1(2), b2(2), b3(2)
-      do 10 k=1,int
-        r0 = b0(k) + b2(k)
-        r1 = b1(k) + b3(k)
-        b2(k) = b0(k) - b2(k)
-        b3(k) = b1(k) - b3(k)
-        b0(k) = r0 + r1
-        b1(k) = r0 - r1
-  10  continue
-      return
-      end
diff --git a/math/ieee/chap1/r4tx.f b/math/ieee/chap1/r4tx.f
deleted file mode 100644
index 4e5649e82..000000000
--- a/math/ieee/chap1/r4tx.f
+++ /dev/null
@@ -1,29 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  r4tx
-c radix 4 iteration subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine r4tx(nthpo, cr0, cr1, cr2, cr3, ci0, ci1, ci2, ci3)
-      dimension cr0(2), cr1(2), cr2(2), cr3(2), ci0(2), ci1(2), ci2(2),
-     *    ci3(2)
-      do 10 k=1,nthpo,4
-        r1 = cr0(k) + cr2(k)
-        r2 = cr0(k) - cr2(k)
-        r3 = cr1(k) + cr3(k)
-        r4 = cr1(k) - cr3(k)
-        fi1 = ci0(k) + ci2(k)
-        fi2 = ci0(k) - ci2(k)
-        fi3 = ci1(k) + ci3(k)
-        fi4 = ci1(k) - ci3(k)
-        cr0(k) = r1 + r3
-        ci0(k) = fi1 + fi3
-        cr1(k) = r1 - r3
-        ci1(k) = fi1 - fi3
-        cr2(k) = r2 - fi4
-        ci2(k) = fi2 + r4
-        cr3(k) = r2 + fi4
-        ci3(k) = fi2 - r4
-  10  continue
-      return
-      end
diff --git a/math/ieee/chap1/r8syn.f b/math/ieee/chap1/r8syn.f
deleted file mode 100644
index 054413d71..000000000
--- a/math/ieee/chap1/r8syn.f
+++ /dev/null
@@ -1,186 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  r8syn
-c radix 8 synthesis subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine r8syn(int, nn, br0, br1, br2, br3, br4, br5, br6, br7,
-     *    bi0, bi1, bi2, bi3, bi4, bi5, bi6, bi7)
-      dimension l(15), br0(2), br1(2), br2(2), br3(2), br4(2), br5(2),
-     *    br6(2), br7(2), bi0(2), bi1(2), bi2(2), bi3(2), bi4(2),
-     *    bi5(2), bi6(2), bi7(2)
-      common /con1/ pii, p7, p7two, c22, s22, pi2
-      equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
-     *    (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
-     *    (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
-     *    (l1,l(15))
-      l(1) = nn/8
-      do 40 k=2,15
-        if (l(k-1)-2) 10, 20, 30
-  10    l(k-1) = 2
-  20    l(k) = 2
-        go to 40
-  30    l(k) = l(k-1)/2
-  40  continue
-      piovn = pii/float(nn)
-      ji = 3
-      jl = 2
-      jr = 2
-c
-      do 120 j1=2,l1,2
-      do 120 j2=j1,l2,l1
-      do 120 j3=j2,l3,l2
-      do 120 j4=j3,l4,l3
-      do 120 j5=j4,l5,l4
-      do 120 j6=j5,l6,l5
-      do 120 j7=j6,l7,l6
-      do 120 j8=j7,l8,l7
-      do 120 j9=j8,l9,l8
-      do 120 j10=j9,l10,l9
-      do 120 j11=j10,l11,l10
-      do 120 j12=j11,l12,l11
-      do 120 j13=j12,l13,l12
-      do 120 j14=j13,l14,l13
-      do 120 jthet=j14,l15,l14
-        th2 = jthet - 2
-        if (th2) 50, 50, 90
-  50    do 60 k=1,int
-          t0 = br0(k) + br1(k)
-          t1 = br0(k) - br1(k)
-          t2 = br2(k) + br2(k)
-          t3 = br3(k) + br3(k)
-          t4 = br4(k) + br6(k)
-          t6 = br7(k) - br5(k)
-          t5 = br4(k) - br6(k)
-          t7 = br7(k) + br5(k)
-          pr = p7*(t7+t5)
-          pi = p7*(t7-t5)
-          tt0 = t0 + t2
-          tt1 = t1 + t3
-          t2 = t0 - t2
-          t3 = t1 - t3
-          t4 = t4 + t4
-          t5 = pr + pr
-          t6 = t6 + t6
-          t7 = pi + pi
-          br0(k) = tt0 + t4
-          br1(k) = tt1 + t5
-          br2(k) = t2 + t6
-          br3(k) = t3 + t7
-          br4(k) = tt0 - t4
-          br5(k) = tt1 - t5
-          br6(k) = t2 - t6
-          br7(k) = t3 - t7
-  60    continue
-        if (nn-8) 120, 120, 70
-  70    k0 = int*8 + 1
-        kl = k0 + int - 1
-        do 80 k=k0,kl
-          t1 = bi0(k) + bi6(k)
-          t2 = bi7(k) - bi1(k)
-          t3 = bi0(k) - bi6(k)
-          t4 = bi7(k) + bi1(k)
-          pr = t3*c22 + t4*s22
-          pi = t4*c22 - t3*s22
-          t5 = bi2(k) + bi4(k)
-          t6 = bi5(k) - bi3(k)
-          t7 = bi2(k) - bi4(k)
-          t8 = bi5(k) + bi3(k)
-          rr = t8*c22 - t7*s22
-          ri = -t8*s22 - t7*c22
-          bi0(k) = (t1+t5) + (t1+t5)
-          bi4(k) = (t2+t6) + (t2+t6)
-          bi1(k) = (pr+rr) + (pr+rr)
-          bi5(k) = (pi+ri) + (pi+ri)
-          t5 = t1 - t5
-          t6 = t2 - t6
-          bi2(k) = p7two*(t6+t5)
-          bi6(k) = p7two*(t6-t5)
-          rr = pr - rr
-          ri = pi - ri
-          bi3(k) = p7two*(ri+rr)
-          bi7(k) = p7two*(ri-rr)
-  80    continue
-        go to 120
-  90    arg = th2*piovn
-        c1 = cos(arg)
-        s1 = -sin(arg)
-        c2 = c1**2 - s1**2
-        s2 = c1*s1 + c1*s1
-        c3 = c1*c2 - s1*s2
-        s3 = c2*s1 + s2*c1
-        c4 = c2**2 - s2**2
-        s4 = c2*s2 + c2*s2
-        c5 = c2*c3 - s2*s3
-        s5 = c3*s2 + s3*c2
-        c6 = c3**2 - s3**2
-        s6 = c3*s3 + c3*s3
-        c7 = c3*c4 - s3*s4
-        s7 = c4*s3 + s4*c3
-        int8 = int*8
-        j0 = jr*int8 + 1
-        k0 = ji*int8 + 1
-        jlast = j0 + int - 1
-        do 100 j=j0,jlast
-          k = k0 + j - j0
-          tr0 = br0(j) + bi6(k)
-          ti0 = bi7(k) - br1(j)
-          tr1 = br0(j) - bi6(k)
-          ti1 = bi7(k) + br1(j)
-          tr2 = br2(j) + bi4(k)
-          ti2 = bi5(k) - br3(j)
-          tr3 = bi5(k) + br3(j)
-          ti3 = bi4(k) - br2(j)
-          tr4 = br4(j) + bi2(k)
-          ti4 = bi3(k) - br5(j)
-          t0 = br4(j) - bi2(k)
-          t1 = bi3(k) + br5(j)
-          tr5 = p7*(t0+t1)
-          ti5 = p7*(t1-t0)
-          tr6 = br6(j) + bi0(k)
-          ti6 = bi1(k) - br7(j)
-          t0 = br6(j) - bi0(k)
-          t1 = bi1(k) + br7(j)
-          tr7 = -p7*(t0-t1)
-          ti7 = -p7*(t1+t0)
-          t0 = tr0 + tr2
-          t1 = ti0 + ti2
-          t2 = tr1 + tr3
-          t3 = ti1 + ti3
-          tr2 = tr0 - tr2
-          ti2 = ti0 - ti2
-          tr3 = tr1 - tr3
-          ti3 = ti1 - ti3
-          t4 = tr4 + tr6
-          t5 = ti4 + ti6
-          t6 = tr5 + tr7
-          t7 = ti5 + ti7
-          ttr6 = ti4 - ti6
-          ti6 = tr6 - tr4
-          ttr7 = ti5 - ti7
-          ti7 = tr7 - tr5
-          br0(j) = t0 + t4
-          bi0(k) = t1 + t5
-          br1(j) = c1*(t2+t6) - s1*(t3+t7)
-          bi1(k) = c1*(t3+t7) + s1*(t2+t6)
-          br2(j) = c2*(tr2+ttr6) - s2*(ti2+ti6)
-          bi2(k) = c2*(ti2+ti6) + s2*(tr2+ttr6)
-          br3(j) = c3*(tr3+ttr7) - s3*(ti3+ti7)
-          bi3(k) = c3*(ti3+ti7) + s3*(tr3+ttr7)
-          br4(j) = c4*(t0-t4) - s4*(t1-t5)
-          bi4(k) = c4*(t1-t5) + s4*(t0-t4)
-          br5(j) = c5*(t2-t6) - s5*(t3-t7)
-          bi5(k) = c5*(t3-t7) + s5*(t2-t6)
-          br6(j) = c6*(tr2-ttr6) - s6*(ti2-ti6)
-          bi6(k) = c6*(ti2-ti6) + s6*(tr2-ttr6)
-          br7(j) = c7*(tr3-ttr7) - s7*(ti3-ti7)
-          bi7(k) = c7*(ti3-ti7) + s7*(tr3-ttr7)
- 100    continue
-        jr = jr + 2
-        ji = ji - 2
-        if (ji-jl) 110, 110, 120
- 110    ji = 2*jr - 1
-        jl = jr
- 120  continue
-      return
-      end
diff --git a/math/ieee/chap1/r8tr.f b/math/ieee/chap1/r8tr.f
deleted file mode 100644
index d49ef58cc..000000000
--- a/math/ieee/chap1/r8tr.f
+++ /dev/null
@@ -1,201 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: r8tr
-c radix 8 iteration subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine r8tr(int, nn, br0, br1, br2, br3, br4, br5, br6, br7,
-     *    bi0, bi1, bi2, bi3, bi4, bi5, bi6, bi7)
-      dimension l(15), br0(2), br1(2), br2(2), br3(2), br4(2), br5(2),
-     *    br6(2), br7(2), bi0(2), bi1(2), bi2(2), bi3(2), bi4(2),
-     *    bi5(2), bi6(2), bi7(2)
-      common /con/ pii, p7, p7two, c22, s22, pi2
-      equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
-     *    (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
-     *    (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
-     *    (l1,l(15))
-c
-c set up counters such that jthet steps through the arguments
-c of w, jr steps through starting locations for the real part of the
-c intermediate results and ji steps through starting locations
-c of the imaginary part of the intermediate results.
-c
-      l(1) = nn/8
-      do 40 k=2,15
-        if (l(k-1)-2) 10, 20, 30
-  10    l(k-1) = 2
-  20    l(k) = 2
-        go to 40
-  30    l(k) = l(k-1)/2
-  40  continue
-      piovn = pii/float(nn)
-      ji = 3
-      jl = 2
-      jr = 2
-      do 120 j1=2,l1,2
-      do 120 j2=j1,l2,l1
-      do 120 j3=j2,l3,l2
-      do 120 j4=j3,l4,l3
-      do 120 j5=j4,l5,l4
-      do 120 j6=j5,l6,l5
-      do 120 j7=j6,l7,l6
-      do 120 j8=j7,l8,l7
-      do 120 j9=j8,l9,l8
-      do 120 j10=j9,l10,l9
-      do 120 j11=j10,l11,l10
-      do 120 j12=j11,l12,l11
-      do 120 j13=j12,l13,l12
-      do 120 j14=j13,l14,l13
-      do 120 jthet=j14,l15,l14
-        th2 = jthet - 2
-        if (th2) 50, 50, 90
-  50    do 60 k=1,int
-          t0 = br0(k) + br4(k)
-          t1 = br1(k) + br5(k)
-          t2 = br2(k) + br6(k)
-          t3 = br3(k) + br7(k)
-          t4 = br0(k) - br4(k)
-          t5 = br1(k) - br5(k)
-          t6 = br2(k) - br6(k)
-          t7 = br3(k) - br7(k)
-          br2(k) = t0 - t2
-          br3(k) = t1 - t3
-          t0 = t0 + t2
-          t1 = t1 + t3
-          br0(k) = t0 + t1
-          br1(k) = t0 - t1
-          pr = p7*(t5-t7)
-          pi = p7*(t5+t7)
-          br4(k) = t4 + pr
-          br7(k) = t6 + pi
-          br6(k) = t4 - pr
-          br5(k) = pi - t6
-  60    continue
-        if (nn-8) 120, 120, 70
-  70    k0 = int*8 + 1
-        kl = k0 + int - 1
-        do 80 k=k0,kl
-          pr = p7*(bi2(k)-bi6(k))
-          pi = p7*(bi2(k)+bi6(k))
-          tr0 = bi0(k) + pr
-          ti0 = bi4(k) + pi
-          tr2 = bi0(k) - pr
-          ti2 = bi4(k) - pi
-          pr = p7*(bi3(k)-bi7(k))
-          pi = p7*(bi3(k)+bi7(k))
-          tr1 = bi1(k) + pr
-          ti1 = bi5(k) + pi
-          tr3 = bi1(k) - pr
-          ti3 = bi5(k) - pi
-          pr = tr1*c22 - ti1*s22
-          pi = ti1*c22 + tr1*s22
-          bi0(k) = tr0 + pr
-          bi6(k) = tr0 - pr
-          bi7(k) = ti0 + pi
-          bi1(k) = pi - ti0
-          pr = -tr3*s22 - ti3*c22
-          pi = tr3*c22 - ti3*s22
-          bi2(k) = tr2 + pr
-          bi4(k) = tr2 - pr
-          bi5(k) = ti2 + pi
-          bi3(k) = pi - ti2
-  80    continue
-        go to 120
-  90    arg = th2*piovn
-        c1 = cos(arg)
-        s1 = sin(arg)
-        c2 = c1**2 - s1**2
-        s2 = c1*s1 + c1*s1
-        c3 = c1*c2 - s1*s2
-        s3 = c2*s1 + s2*c1
-        c4 = c2**2 - s2**2
-        s4 = c2*s2 + c2*s2
-        c5 = c2*c3 - s2*s3
-        s5 = c3*s2 + s3*c2
-        c6 = c3**2 - s3**2
-        s6 = c3*s3 + c3*s3
-        c7 = c3*c4 - s3*s4
-        s7 = c4*s3 + s4*c3
-        int8 = int*8
-        j0 = jr*int8 + 1
-        k0 = ji*int8 + 1
-        jlast = j0 + int - 1
-        do 100 j=j0,jlast
-          k = k0 + j - j0
-          tr1 = br1(j)*c1 - bi1(k)*s1
-          ti1 = br1(j)*s1 + bi1(k)*c1
-          tr2 = br2(j)*c2 - bi2(k)*s2
-          ti2 = br2(j)*s2 + bi2(k)*c2
-          tr3 = br3(j)*c3 - bi3(k)*s3
-          ti3 = br3(j)*s3 + bi3(k)*c3
-          tr4 = br4(j)*c4 - bi4(k)*s4
-          ti4 = br4(j)*s4 + bi4(k)*c4
-          tr5 = br5(j)*c5 - bi5(k)*s5
-          ti5 = br5(j)*s5 + bi5(k)*c5
-          tr6 = br6(j)*c6 - bi6(k)*s6
-          ti6 = br6(j)*s6 + bi6(k)*c6
-          tr7 = br7(j)*c7 - bi7(k)*s7
-          ti7 = br7(j)*s7 + bi7(k)*c7
-c
-          t0 = br0(j) + tr4
-          t1 = bi0(k) + ti4
-          tr4 = br0(j) - tr4
-          ti4 = bi0(k) - ti4
-          t2 = tr1 + tr5
-          t3 = ti1 + ti5
-          tr5 = tr1 - tr5
-          ti5 = ti1 - ti5
-          t4 = tr2 + tr6
-          t5 = ti2 + ti6
-          tr6 = tr2 - tr6
-          ti6 = ti2 - ti6
-          t6 = tr3 + tr7
-          t7 = ti3 + ti7
-          tr7 = tr3 - tr7
-          ti7 = ti3 - ti7
-c
-          tr0 = t0 + t4
-          ti0 = t1 + t5
-          tr2 = t0 - t4
-          ti2 = t1 - t5
-          tr1 = t2 + t6
-          ti1 = t3 + t7
-          tr3 = t2 - t6
-          ti3 = t3 - t7
-          t0 = tr4 - ti6
-          t1 = ti4 + tr6
-          t4 = tr4 + ti6
-          t5 = ti4 - tr6
-          t2 = tr5 - ti7
-          t3 = ti5 + tr7
-          t6 = tr5 + ti7
-          t7 = ti5 - tr7
-          br0(j) = tr0 + tr1
-          bi7(k) = ti0 + ti1
-          bi6(k) = tr0 - tr1
-          br1(j) = ti1 - ti0
-          br2(j) = tr2 - ti3
-          bi5(k) = ti2 + tr3
-          bi4(k) = tr2 + ti3
-          br3(j) = tr3 - ti2
-          pr = p7*(t2-t3)
-          pi = p7*(t2+t3)
-          br4(j) = t0 + pr
-          bi3(k) = t1 + pi
-          bi2(k) = t0 - pr
-          br5(j) = pi - t1
-          pr = -p7*(t6+t7)
-          pi = p7*(t6-t7)
-          br6(j) = t4 + pr
-          bi1(k) = t5 + pi
-          bi0(k) = t4 - pr
-          br7(j) = pi - t5
- 100    continue
-        jr = jr + 2
-        ji = ji - 2
-        if (ji-jl) 110, 110, 120
- 110    ji = 2*jr - 1
-        jl = jr
- 120  continue
-      return
-      end
diff --git a/math/ieee/chap1/r8tx.f b/math/ieee/chap1/r8tx.f
deleted file mode 100644
index 9cc4f591f..000000000
--- a/math/ieee/chap1/r8tx.f
+++ /dev/null
@@ -1,107 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  r8tx
-c radix 8 iteration subroutine
-c-----------------------------------------------------------------------
-c
-      subroutine r8tx(nxtlt, nthpo, lengt, cr0, cr1, cr2, cr3, cr4,
-     *    cr5, cr6, cr7, ci0, ci1, ci2, ci3, ci4, ci5, ci6, ci7)
-      dimension cr0(2), cr1(2), cr2(2), cr3(2), cr4(2), cr5(2), cr6(2),
-     *    cr7(2), ci1(2), ci2(2), ci3(2), ci4(2), ci5(2), ci6(2),
-     *    ci7(2), ci0(2)
-      common /con2/ pi2, p7
-c
-      scale = pi2/float(lengt)
-      do 30 j=1,nxtlt
-        arg = float(j-1)*scale
-        c1 = cos(arg)
-        s1 = sin(arg)
-        c2 = c1**2 - s1**2
-        s2 = c1*s1 + c1*s1
-        c3 = c1*c2 - s1*s2
-        s3 = c2*s1 + s2*c1
-        c4 = c2**2 - s2**2
-        s4 = c2*s2 + c2*s2
-        c5 = c2*c3 - s2*s3
-        s5 = c3*s2 + s3*c2
-        c6 = c3**2 - s3**2
-        s6 = c3*s3 + c3*s3
-        c7 = c3*c4 - s3*s4
-        s7 = c4*s3 + s4*c3
-        do 20 k=j,nthpo,lengt
-          ar0 = cr0(k) + cr4(k)
-          ar1 = cr1(k) + cr5(k)
-          ar2 = cr2(k) + cr6(k)
-          ar3 = cr3(k) + cr7(k)
-          ar4 = cr0(k) - cr4(k)
-          ar5 = cr1(k) - cr5(k)
-          ar6 = cr2(k) - cr6(k)
-          ar7 = cr3(k) - cr7(k)
-          ai0 = ci0(k) + ci4(k)
-          ai1 = ci1(k) + ci5(k)
-          ai2 = ci2(k) + ci6(k)
-          ai3 = ci3(k) + ci7(k)
-          ai4 = ci0(k) - ci4(k)
-          ai5 = ci1(k) - ci5(k)
-          ai6 = ci2(k) - ci6(k)
-          ai7 = ci3(k) - ci7(k)
-          br0 = ar0 + ar2
-          br1 = ar1 + ar3
-          br2 = ar0 - ar2
-          br3 = ar1 - ar3
-          br4 = ar4 - ai6
-          br5 = ar5 - ai7
-          br6 = ar4 + ai6
-          br7 = ar5 + ai7
-          bi0 = ai0 + ai2
-          bi1 = ai1 + ai3
-          bi2 = ai0 - ai2
-          bi3 = ai1 - ai3
-          bi4 = ai4 + ar6
-          bi5 = ai5 + ar7
-          bi6 = ai4 - ar6
-          bi7 = ai5 - ar7
-          cr0(k) = br0 + br1
-          ci0(k) = bi0 + bi1
-          if (j.le.1) go to 10
-          cr1(k) = c4*(br0-br1) - s4*(bi0-bi1)
-          ci1(k) = c4*(bi0-bi1) + s4*(br0-br1)
-          cr2(k) = c2*(br2-bi3) - s2*(bi2+br3)
-          ci2(k) = c2*(bi2+br3) + s2*(br2-bi3)
-          cr3(k) = c6*(br2+bi3) - s6*(bi2-br3)
-          ci3(k) = c6*(bi2-br3) + s6*(br2+bi3)
-          tr = p7*(br5-bi5)
-          ti = p7*(br5+bi5)
-          cr4(k) = c1*(br4+tr) - s1*(bi4+ti)
-          ci4(k) = c1*(bi4+ti) + s1*(br4+tr)
-          cr5(k) = c5*(br4-tr) - s5*(bi4-ti)
-          ci5(k) = c5*(bi4-ti) + s5*(br4-tr)
-          tr = -p7*(br7+bi7)
-          ti = p7*(br7-bi7)
-          cr6(k) = c3*(br6+tr) - s3*(bi6+ti)
-          ci6(k) = c3*(bi6+ti) + s3*(br6+tr)
-          cr7(k) = c7*(br6-tr) - s7*(bi6-ti)
-          ci7(k) = c7*(bi6-ti) + s7*(br6-tr)
-          go to 20
-  10      cr1(k) = br0 - br1
-          ci1(k) = bi0 - bi1
-          cr2(k) = br2 - bi3
-          ci2(k) = bi2 + br3
-          cr3(k) = br2 + bi3
-          ci3(k) = bi2 - br3
-          tr = p7*(br5-bi5)
-          ti = p7*(br5+bi5)
-          cr4(k) = br4 + tr
-          ci4(k) = bi4 + ti
-          cr5(k) = br4 - tr
-          ci5(k) = bi4 - ti
-          tr = -p7*(br7+bi7)
-          ti = p7*(br7-bi7)
-          cr6(k) = br6 + tr
-          ci6(k) = bi6 + ti
-          cr7(k) = br6 - tr
-          ci7(k) = bi6 - ti
-  20    continue
-  30  continue
-      return
-      end
diff --git a/math/ieee/chap1/rad4sb.f b/math/ieee/chap1/rad4sb.f
deleted file mode 100644
index dfebf0cca..000000000
--- a/math/ieee/chap1/rad4sb.f
+++ /dev/null
@@ -1,38 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  rad4sb
-c used by subroutine radix4. never directly accessed by user.
-c-----------------------------------------------------------------------
-c
-       subroutine rad4sb(ntype)
-c
-c      input: ntype = type of butterfly invoked
-c      output: parameters used by subroutine radix4
-c
-       dimension ix(2996)
-       common /xx/ix
-       common ntypl,kkp,index,ixc
-       if(ntype.eq.ntypl) go to 7
-       ix(ixc)=0
-       ix(ixc+1)=ntype
-       ixc=ixc+2
-       if(ntype.ne.4) go to 4
-       indexp=(index-1)*9
-       ix(ixc)=kkp+1
-       ix(ixc+1)=indexp+1
-       ixc=ixc+2
-       go to 6
-4      ix(ixc)=kkp+1
-       ixc=ixc+1
-6      ntypl=ntype
-       return
-7      if(ntype.ne.4) go to 8
-       indexp=(index-1)*9
-       ix(ixc)=kkp+1
-       ix(ixc+1)=indexp+1
-       ixc=ixc+2
-       return
-8      ix(ixc)=kkp+1
-       ixc=ixc+1
-       return
-       end
diff --git a/math/ieee/chap1/radix4.f b/math/ieee/chap1/radix4.f
deleted file mode 100644
index cd7033059..000000000
--- a/math/ieee/chap1/radix4.f
+++ /dev/null
@@ -1,488 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine:  radix4
-c computes forward or inverse complex dft via radix-4 fft.
-c uses autogen technique to yield time efficient program.
-c-----------------------------------------------------------------------
-c
-        subroutine radix4(mm,iflag,jflag)
-c
-c       input:
-c             mm = power of 4 (i.e., n = 4**mm complex point transform)
-c             (mm.ge.2 and mm.le.5)
-c
-c          iflag = 1 on first pass for given n
-c                = 0 on subsequent passes for given n
-c
-c          jflag = -1 for forward transform
-c                = +1 for inverse transform
-c
-c       input/output:
-c              a = array of dimensions 2*n with real and imaginary parts
-c                  of dft input/output in odd, even array components.
-c
-c       for optimal time efficiency, common is used to pass arrays.
-c       this means that dimensions of arrays a, ix, and t can be
-c       modified to reflect maximum value of n = 4**mm to be used. note
-c       that array "ix" is also dimensioned in subroutine "rad4sb".
-c
-c    i.e.,    a(    ) ix(  ) t(   )
-c
-c       m =2      32     38     27
-c       m<=3     128    144    135
-c       m<=4     512    658    567
-c       m<=5    2048   2996   2295
-c
-       dimension a(2048),ix(2996),t(2295)
-       dimension nfac(11),np(209)
-       common ntypl,kkp,index,ixc
-       common /aa/a
-       common /xx/ix
-c
-c      check for mm<2 or mm>5
-c
-       if(mm.lt.2.or.mm.gt.5)stop
-c
-c      initialize on first pass """"""""""""""""""""""""""""""""""""""""
-c
-       if(iflag.eq.1) go to 9999
-c
-c      fast fourier transform start ####################################
-c
-8885   kspan=2*4**mm
-       if(jflag.eq.1) go to 8887
-c
-c      conjugate data for forward transform
-c
-       do 8886 j=2,n2,2
-8886   a(j)=-a(j)
-       go to 8889
-c
-c      multiply data by n**(-1) if inverse transform
-c
-8887   do 8888 j=1,n2,2
-       a(j)=a(j)*xp
-8888   a(j+1)=a(j+1)*xp
-8889   i=3
-       it=ix(i-1)
-       go to (1,2,3,4,5,6,7,8),it
-c***********************************************************************
-c
-c                                    8 multiply butterfly
-c
-1      kk=ix(i)
-c
-11     k1=kk+kspan
-       k2=k1+kspan
-       k3=k2+kspan
-c
-       akp=a(kk)+a(k2)
-       akm=a(kk)-a(k2)
-       ajp=a(k1)+a(k3)
-       ajm=a(k1)-a(k3)
-       a(kk)=akp+ajp
-c
-       bkp=a(kk+1)+a(k2+1)
-       bkm=a(kk+1)-a(k2+1)
-       bjp=a(k1+1)+a(k3+1)
-       bjm=a(k1+1)-a(k3+1)
-       a(kk+1)=bkp+bjp
-c
-       bjp=bkp-bjp
-c
-       a(k2+1)=(akp+bjp-ajp)*c707
-       a(k2)=a(k2+1)+bjp*cm141
-c
-       bkp=bkm+ajm
-       akp=akm-bjm
-c
-       ac0=(akp+bkp)*c924
-       a(k1+1)=ac0+akp*cm541
-       a(k1)  =ac0+bkp*cm131
-c
-       bkm=bkm-ajm
-       akm=akm+bjm
-c
-       ac0=(akm+bkm)*c383
-       a(k3+1)=ac0+akm*c541
-       a(k3)  =ac0+bkm*cm131
-c
-       i=i+1
-       kk=ix(i)
-       if (kk) 111,111,11
-111    i=i+2
-       it=ix(i-1)
-       go to (1,2,3,4,5,6,7,8), it
-c***********************************************************************
-c
-c                                    4 multiply butterfly
-c
-2      kk=ix(i)
-c
-22     k1=kk+kspan
-       k2=k1+kspan
-       k3=k2+kspan
-c
-       akp=a(kk)+a(k2)
-       akm=a(kk)-a(k2)
-       ajp=a(k1)+a(k3)
-       ajm=a(k1)-a(k3)
-       a(kk)=akp+ajp
-c
-       bkp=a(kk+1)+a(k2+1)
-       bkm=a(kk+1)-a(k2+1)
-       bjp=a(k1+1)+a(k3+1)
-       bjm=a(k1+1)-a(k3+1)
-       a(kk+1)=bkp+bjp
-       a(k2)=-bkp+bjp
-       a(k2+1)=akp-ajp
-c
-       bkp=bkm+ajm
-c
-       a(k1+1)=(bkp+akm-bjm)*c707
-       a(k1)=a(k1+1)+bkp*cm141
-c
-       akm=akm+bjm
-c
-       a(k3+1)=(akm+ajm-bkm)*c707
-       a(k3)=a(k3+1)+akm*cm141
-c
-       i=i+1
-       kk=ix(i)
-       if (kk) 222,222,22
-222    i=i+2
-       it=ix(i-1)
-       go to (1,2,3,4,5,6,7,8), it
-c***********************************************************************
-c
-c                                    8 multiply butterfly
-c
-3      kk=ix(i)
-c
-33     k1=kk+kspan
-       k2=k1+kspan
-       k3=k2+kspan
-c
-       akp=a(kk)+a(k2)
-       akm=a(kk)-a(k2)
-       ajp=a(k1)+a(k3)
-       ajm=a(k1)-a(k3)
-       a(kk)=akp+ajp
-c
-       bkp=a(kk+1)+a(k2+1)
-       bkm=a(kk+1)-a(k2+1)
-       bjp=a(k1+1)+a(k3+1)
-       bjm=a(k1+1)-a(k3+1)
-       a(kk+1)=bkp+bjp
-c
-       ajp=akp-ajp
-c
-       a(k2+1)=(ajp+bjp-bkp)*c707
-       a(k2)=a(k2+1)+ajp*cm141
-c
-       bkp=bkm+ajm
-       akp=akm-bjm
-c
-       ac0=(akp+bkp)*c383
-       a(k1+1)=ac0+akp*c541
-       a(k1)  =ac0+bkp*cm131
-c
-       bkm=bkm-ajm
-       akm=akm+bjm
-c
-       ac0=(akm+bkm)*cm924
-       a(k3+1)=ac0+akm*c541
-       a(k3)  =ac0+bkm*c131
-c
-       i=i+1
-       kk=ix(i)
-       if (kk) 333,333,33
-333    i=i+2
-       it=ix(i-1)
-       go to (1,2,3,4,5,6,7,8), it
-c***********************************************************************
-c
-c                                    general 9 multiply butterfly
-c
-4      kk=ix(i)
-c
-44     k1=kk+kspan
-       k2=k1+kspan
-       k3=k2+kspan
-c
-       akp=a(kk)+a(k2)
-       akm=a(kk)-a(k2)
-       ajp=a(k1)+a(k3)
-       ajm=a(k1)-a(k3)
-       a(kk)=akp+ajp
-c
-       bkp=a(kk+1)+a(k2+1)
-       bkm=a(kk+1)-a(k2+1)
-       bjp=a(k1+1)+a(k3+1)
-       bjm=a(k1+1)-a(k3+1)
-       a(kk+1)=bkp+bjp
-c
-       ajp=akp-ajp
-       bjp=bkp-bjp
-c
-       j=ix(i+1)
-c
-       ac0=(ajp+bjp)*t(j+8)
-       a(k2+1)=ac0+ajp*t(j+6)
-       a(k2)  =ac0+bjp*t(j+7)
-c
-       bkp=bkm+ajm
-       akp=akm-bjm
-c
-       ac0=(akp+bkp)*t(j+5)
-       a(k1+1)=ac0+akp*t(j+3)
-       a(k1)  =ac0+bkp*t(j+4)
-c
-       bkm=bkm-ajm
-       akm=akm+bjm
-c
-       ac0=(akm+bkm)*t(j+2)
-       a(k3+1)=ac0+akm*t(j)
-       a(k3)  =ac0+bkm*t(j+1)
-c
-       i=i+2
-       kk=ix(i)
-       if (kk) 444,444,44
-444    i=i+2
-       it=ix(i-1)
-       go to (1,2,3,4,5,6,7,8), it
-c***********************************************************************
-c
-c                                    0 multiply butterfly
-c
-5      kk=ix(i)
-c
-55     k1=kk+kspan
-       k2=k1+kspan
-       k3=k2+kspan
-c
-       akp=a(kk)+a(k2)
-       akm=a(kk)-a(k2)
-       ajp=a(k1)+a(k3)
-       ajm=a(k1)-a(k3)
-       a(kk)=akp+ajp
-       a(k2)=akp-ajp
-c
-       bkp=a(kk+1)+a(k2+1)
-       bkm=a(kk+1)-a(k2+1)
-       bjp=a(k1+1)+a(k3+1)
-       bjm=a(k1+1)-a(k3+1)
-       a(kk+1)=bkp+bjp
-       a(k2+1)=bkp-bjp
-c
-       a(k3+1)=bkm-ajm
-       a(k1+1)=bkm+ajm
-       a(k3)=akm+bjm
-       a(k1)=akm-bjm
-c
-       i=i+1
-       kk=ix(i)
-       if (kk) 555,555,55
-555    i=i+2
-       it=ix(i-1)
-       go to (1,2,3,4,5,6,7,8), it
-c***********************************************************************
-c
-c                                       offset reduced
-c
-6      kspan=kspan/4
-       i=i+2
-       it=ix(i-1)
-       go to (1,2,3,4,5,6,7,8), it
-c***********************************************************************
-c
-c                                       bit reversal (shuffling)
-c
-7      ip1=ix(i)
-77     ip2=ix(i+1)
-       t1=a(ip2)
-       a(ip2)=a(ip1)
-       a(ip1)=t1
-       t1=a(ip2+1)
-       a(ip2+1)=a(ip1+1)
-       a(ip1+1)=t1
-       i=i+2
-       ip1=ix(i)
-       if (ip1) 777,777,77
-777    i=i+2
-       it=ix(i-1)
-       go to (1,2,3,4,5,6,7,8), it
-c***********************************************************************
-8      if(jflag.eq.1) go to 888
-c
-c      conjugate output if forward transform
-c
-       do 88 j=2,n2,2
-88     a(j)=-a(j)
-888    return
-c
-c      fast fourier transform ends #####################################
-c
-c      initialization phase starts. done only once
-c
-9999   ixc=1
-       n=4**mm
-       xp=n
-       xp=1./xp
-       ntot=n
-       n2=n*2
-       nspan=n
-       n1test=n/16
-       n2test=n/8
-       n3test=(3*n)/16
-       nspan4=nspan/4
-       ibase=0
-       isn=1
-       inc=isn
-       rad=8.0*atan(1.0)
-       pi=4.*atan(1.0)
-       c707=sin(pi/4.)
-       cm141=-2.*c707
-       c383=sin(pi/8.)
-       c924=cos(pi/8.)
-       cm924=-c924
-       c541=c924-c383
-       cm541=-c541
-       c131=c924+c383
-       cm131=-c131
-10     nt=inc*ntot
-       ks=inc*nspan
-       kspan=ks
-       jc=ks/n
-       radf=rad*float(jc)*.5
-       i=0
-c
-c      determine the factors of n
-c      all factors must be 4 for this version
-c
-       m=0
-       k=n
-15     m=m+1
-       nfac(m)=4
-       k=k/4
-20     if(k-(k/4)*4.eq.0) go to 15
-       kt=1
-       if(n.ge.256) kt=2
-       kspan0=kspan
-       ntypl=0
-c
-100    ndelta=kspan0/kspan
-       index=0
-       sd=radf/float(kspan)
-       cd=2.0*sin(sd)**2
-       sd=sin(sd+sd)
-       kk=1
-       i=i+1
-c
-c      transform for a factor of 4
-c
-       kspan=kspan/4
-       ix(ixc)=0
-       ix(ixc+1)=6
-       ixc=ixc+2
-c
-410    c1=1.0
-       s1=0.0
-420    k1=kk+kspan
-       k2=k1+kspan
-       k3=k2+kspan
-       if(s1.eq.0.0) go to 460
-430    if(kspan.ne.nspan4) go to 431
-       t(ibase+5)=-(s1+c1)
-       t(ibase+6)=c1
-       t(ibase+4)=s1-c1
-       t(ibase+8)=-(s2+c2)
-       t(ibase+9)=c2
-       t(ibase+7)=s2-c2
-       t(ibase+2)=-(s3+c3)
-       t(ibase+3)=c3
-       t(ibase+1)=s3-c3
-       ibase=ibase+9
-c
-431    kkp=(kk-1)*2
-       if(index.ne.n1test) go to 150
-       call rad4sb(1)
-       go to 5035
-150    if(index.ne.n2test) go to 160
-       call rad4sb(2)
-       go to 5035
-160    if(index.ne.n3test) go to 170
-       call rad4sb(3)
-       go to 5035
-170    call rad4sb(4)
-5035   kk=k3+kspan
-       if(kk.le.nt) go to 420
-440    index=index+ndelta
-       c2=c1-(cd*c1+sd*s1)
-       s1=(sd*c1-cd*s1)+s1
-       c1=c2
-       c2=c1*c1-s1*s1
-       s2=c1*s1+c1*s1
-       c3=c2*c1-s2*s1
-       s3=c2*s1+s2*c1
-       kk=kk-nt+jc
-       if(kk.le.kspan) go to 420
-       kk=kk-kspan+inc
-       if(kk.le.jc) go to 410
-       if(kspan.eq.jc) go to 800
-       go to 100
-460    kkp=(kk-1)*2
-       call rad4sb(5)
-5050   kk=k3+kspan
-       if(kk.le.nt) go to 420
-       go to 440
-c
-800    ix(ixc)=0
-       ix(ixc+1)=7
-       ixc=ixc+2
-c
-c      compute parameters to permute the results to normal order
-c      done in two steps
-c      permutation for square factors of n
-c
-       np(1)=ks
-       k=kt+kt+1
-       if(m.lt.k) k=k-1
-       j=1
-       np(k+1)=jc
-810    np(j+1)=np(j)/nfac(j)
-       np(k)=np(k+1)*nfac(j)
-       j=j+1
-       k=k-1
-       if(j.lt.k) go to 810
-       k3=np(k+1)
-       kspan=np(2)
-       kk=jc+1
-       k2=kspan+1
-       j=1
-c
-c      permutation for single variate transform
-c
-820    kkp=(kk-1)*2
-       k2p=(k2-1)*2
-       ix(ixc)=kkp+1
-       ix(ixc+1)=k2p+1
-       ixc=ixc+2
-       kk=kk+inc
-       k2=kspan+k2
-       if(k2.lt.ks) go to 820
-830    k2=k2-np(j)
-       j=j+1
-       k2=np(j+1)+k2
-       if(k2.gt.np(j)) go to 830
-       j=1
-840    if(kk.lt.k2) go to 820
-       kk=kk+inc
-       k2=kspan+k2
-       if(k2.lt.ks) go to 840
-       if(kk.lt.ks) go to 830
-       jc=k3
-       ix(ixc)=0
-       ix(ixc+1)=8
-       go to 8885
-       end
diff --git a/math/ieee/chap1/test/test12.f b/math/ieee/chap1/test/test12.f
deleted file mode 100644
index a830ab2cb..000000000
--- a/math/ieee/chap1/test/test12.f
+++ /dev/null
@@ -1,90 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c main program: fastmain  -  fast fourier transforms
-c authors:      g. d. bergland and m. t. dolan
-c               bell laboratories, murray hill, new jersey 07974
-c
-c input:        the program calls on a random number
-c               generator for input and checks dft and
-c               idft with a 32-point sequence
-c-----------------------------------------------------------------------
-c
-      dimension x(32), y(32), b(34)
-c
-c generate random numbers and store array in b so
-c the same sequence can be used in all tests.
-c note that b is dimensioned to size n+2.
-c
-c iw is a machine dependent write device number
-c
-      iw = i1mach(2)
-c
-      do 10 i=1,32
-        x(i) = uni(0)
-        b(i) = x(i)
-  10  continue
-      m = 5
-      n = 2**m
-      np1 = n + 1
-      np2 = n + 2
-      knt = 1
-c
-c test fast-fsst then ffa-ffs
-c
-      write (iw,9999)
-  20  write (iw,9998) (b(i),i=1,n)
-      if (knt.eq.1) call fast(b, n)
-      if (knt.eq.2) call ffa(b, n)
-      write (iw,9997) (b(i),i=1,np1,2)
-      write (iw,9996) (b(i),i=2,np2,2)
-      if (knt.eq.1) call fsst(b, n)
-      if (knt.eq.2) call ffs(b, n)
-      write (iw,9995) (b(i),i=1,n)
-      knt = knt + 1
-      if (knt.eq.3) go to 40
-c
-      write (iw,9994)
-      do 30 i=1,n
-        b(i) = x(i)
-  30  continue
-      go to 20
-c
-c test fft842 with real input then complex
-c
-  40  write (iw,9993)
-      do 50 i=1,n
-        b(i) = x(i)
-        y(i) = 0.
-  50  continue
-  60  write (iw,9992) (b(i),i=1,n)
-      write (iw,9991) (y(i),i=1,n)
-      call fft842(0, n, b, y)
-      write (iw,9997) (b(i),i=1,n)
-      write (iw,9996) (y(i),i=1,n)
-      call fft842(1, n, b, y)
-      write (iw,9990) (b(i),i=1,n)
-      write (iw,9989) (y(i),i=1,n)
-      knt = knt + 1
-      if (knt.eq.5) go to 80
-c
-      write (iw,9988)
-      do 70 i=1,n
-        b(i) = x(i)
-        y(i) = uni(0)
-  70  continue
-      go to 60
-c
-9999  format (19h1test fast and fsst)
-9998  format (20h0real input sequence/(4e17.8))
-9997  format (29h0real components of transform/(4e17.8))
-9996  format (29h0imag components of transform/(4e17.8))
-9995  format (23h0real inverse transform/(4e17.8))
-9994  format (17h1test ffa and ffs)
-9993  format (37h1test fft842 with real input sequence/(4e17.8))
-9992  format (34h0real components of input sequence/(4e17.8))
-9991  format (34h0imag components of input sequence/(4e17.8))
-9990  format (37h0real components of inverse transform/(4e17.8))
-9989  format (37h0imag components of inverse transform/(4e17.8))
-9988  format (40h1test fft842 with complex input sequence)
-  80  stop
-      end
diff --git a/math/ieee/chap1/test/test13.f b/math/ieee/chap1/test/test13.f
deleted file mode 100644
index d71570695..000000000
--- a/math/ieee/chap1/test/test13.f
+++ /dev/null
@@ -1,260 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c main program: test program for fft subroutines
-c author:      l r rabiner
-c      bell laboratories, murray hill, new jersey, 07974
-c input:       randomly chosen sequences to test fft subroutines
-c      for sequences with special properties
-c      n is the fft length (n must be a power of 2)
-c          2<= n <= 4096
-c-----------------------------------------------------------------------
-c
-      dimension x(4098), y(4098)
-c
-c  define i/0 device codes
-c  input: input to this program is user-interactive
-c        that is - a question is written on the user
-c        terminal (iout1) ad the user types in the answer.
-c
-c output: all output is written on the standard
-c         output unit (iout2).
-c
-      ind = i1mach(1)
-      iout1 = i1mach(4)
-      iout2 = i1mach(2)
-c
-c read in analysis size for fft
-c
-  10  write (iout1,9999)
-9999  format (30h fft size(2.le.n.le.4096)(i4)=)
-      read (ind,9998) n
-9998  format (i4)
-      if (n.eq.0) stop
-      do 20 i=1,12
-        itest = 2**i
-        if (n.eq.itest) go to 30
-  20  continue
-      write (iout1,9997)
-9997  format (45h n is not a power of 2 in the range 2 to 4096)
-      go to 10
-  30  write (iout2,9996) n
-9996  format (11h testing n=, i5, 17h random sequences)
-      write (iout2,9992)
-      np2 = n + 2
-      no2 = n/2
-      no2p1 = no2 + 1
-      no4 = n/4
-      no4p1 = no4 + 1
-c
-c create symmetrical sequence of size n
-c
-      do 40 i=2,no2
-        x(i) = uni(0) - 0.5
-        ind1 = np2 - i
-        x(ind1) = x(i)
-  40  continue
-      x(1) = uni(0) - 0.5
-      x(no2p1) = uni(0) - 0.5
-      do 50 i=1,no2p1
-        y(i) = x(i)
-  50  continue
-      write (iout2,9995)
-9995  format (28h original symmetric sequence)
-      write (iout2,9993) (x(i),i=1,n)
-      write (iout2,9992)
-c
-c compute true fft of n point sequence
-c
-      call fast(x, n)
-      write (iout2,9994) n
-9994  format (1h , i4, 32h point fft of symmetric sequence)
-      write (iout2,9993) (x(i),i=1,np2)
-9993  format (1h , 5e13.5)
-      write (iout2,9992)
-9992  format (1h /1h )
-c
-c use subroutine fftsym to obtain dft from no2 point fft
-c
-      do 60 i=1,no2p1
-        x(i) = y(i)
-  60  continue
-      call fftsym(x, n, y)
-      write (iout2,9991)
-9991  format (17h output of fftsym)
-      write (iout2,9993) (x(i),i=1,no2p1)
-      write (iout2,9992)
-c
-c use subroutine iftsym to obtain original sequence from no2 point dft
-c
-      call iftsym(x, n, y)
-      write (iout2,9990)
-9990  format (17h output of iftsym)
-      write (iout2,9993) (x(i),i=1,no2p1)
-      write (iout2,9992)
-c
-c create antisymmetric n point sequence
-c
-      do 70 i=2,no2
-        x(i) = uni(0) - 0.5
-        ind1 = np2 - i
-        x(ind1) = -x(i)
-  70  continue
-      x(1) = 0.
-      x(no2p1) = 0.
-      do 80 i=1,no2p1
-        y(i) = x(i)
-  80  continue
-      write (iout2,9989)
-9989  format (32h original antisymmetric sequence)
-      write (iout2,9993) (x(i),i=1,n)
-      write (iout2,9992)
-c
-c obtain n point dft of antisymmetric sequence
-c
-      call fast(x, n)
-      write (iout2,9988) n
-9988  format (1h , i4, 36h point fft of antisymmetric sequence)
-      write (iout2,9993) (x(i),i=1,np2)
-      write (iout2,9992)
-c
-c use subroutine fftasm to obtain dft from no2 point fft
-c
-      do 90 i=1,no2
-        x(i) = y(i)
-  90  continue
-      call fftasm(x, n, y)
-      write (iout2,9987)
-9987  format (17h output of fftasm)
-      write (iout2,9993) (x(i),i=1,no2p1)
-      write (iout2,9992)
-c
-c use subroutine iftasm to obtain original sequence from no2 point dft
-c
-      call iftasm(x, n, y)
-      write (iout2,9986)
-9986  format (17h output of iftasm)
-      write (iout2,9993) (x(i),i=1,no2)
-      write (iout2,9992)
-c
-c create sequence with only odd harmonics--begin in frequency domain
-c
-      do 100 i=1,np2,2
-        x(i) = 0.
-        x(i+1) = 0.
-        if (mod(i,4).eq.1) go to 100
-        x(i) = uni(0) - 0.5
-        x(i+1) = uni(0) - 0.5
-        if (n.eq.2) x(i+1) = 0.
- 100  continue
-      write (iout2,9985) n
-9985  format (1h , i4, 35h point fft of odd harmonic sequence)
-      write (iout2,9993) (x(i),i=1,np2)
-      write (iout2,9992)
-c
-c transform back to time sequence
-c
-      call fsst(x, n)
-      write (iout2,9984)
-9984  format (31h original odd harmonic sequence)
-      write (iout2,9993) (x(i),i=1,n)
-      write (iout2,9992)
-c
-c use subroutine fftohm to obtain dft from no2 point fft
-c
-      call fftohm(x, n)
-      write (iout2,9983)
-9983  format (17h output of fftohm)
-      write (iout2,9993) (x(i),i=1,no2)
-      write (iout2,9992)
-c
-c use subroutine iftohm to obtain original sequence from no2 point dft
-c
-      call iftohm(x, n)
-      write (iout2,9982)
-9982  format (17h output of iftohm)
-      write (iout2,9993) (x(i),i=1,no2)
-      write (iout2,9992)
-c
-c create sequence with only real valued odd harmonics
-c
-      do 110 i=1,np2,2
-        x(i) = 0.
-        x(i+1) = 0.
-        if (mod(i,4).eq.1) go to 110
-        x(i) = uni(0) - 0.5
- 110  continue
-      write (iout2,9981) n
-9981  format (1h , i4, 45h point fft of odd harmonic, symmetric sequenc,
-     *    1he)
-      write (iout2,9993) (x(i),i=1,np2)
-      write (iout2,9992)
-c
-c transform back to time sequence
-c
-      call fsst(x, n)
-      write (iout2,9980)
-9980  format (42h original odd harmonic, symmetric sequence)
-      write (iout2,9993) (x(i),i=1,n)
-      write (iout2,9992)
-c
-c use subroutine fftsoh to obtain dft from no4 point fft
-c
-      call fftsoh(x, n, y)
-      write (iout2,9979)
-9979  format (17h output of fftsoh)
-      write (iout2,9993) (x(i),i=1,no4)
-      write (iout2,9992)
-c
-c use subroutine iftsoh to obtain original sequence from no4 point dft
-c
-      call iftsoh(x, n, y)
-      write (iout2,9978)
-9978  format (17h output of iftsoh)
-      write (iout2,9993) (x(i),i=1,no4)
-      write (iout2,9992)
-c
-c create sequence with only imaginary valued odd harmonics--begin
-c in frequency domain
-c
-      do 120 i=1,np2,2
-        x(i) = 0.
-        x(i+1) = 0.
-        if (mod(i,4).eq.1) go to 120
-        x(i+1) = uni(0) - 0.5
- 120  continue
-      write (iout2,9977) n
-9977  format (1h , i4, 41h point fft of odd harmonic, antisymmetric,
-     *    9h sequence)
-      write (iout2,9993) (x(i),i=1,np2)
-      write (iout2,9992)
-c
-c transform back to time sequence
-c
-      call fsst(x, n)
-      write (iout2,9976)
-9976  format (46h original odd harmonic, antisymmetric sequence)
-      write (iout2,9993) (x(i),i=1,n)
-      write (iout2,9992)
-c
-c use subroutine fftaoh to obtain dft from no4 point fft
-c
-      call fftaoh(x, n, y)
-      write (iout2,9975)
-9975  format (17h output of fftaoh)
-      write (iout2,9993) (x(i),i=1,no4)
-      write (iout2,9992)
-c
-c use subroutine iftaoh to obtain original sequence from n/4 point dft
-c
-      call iftaoh(x, n, y)
-      write (iout2,9974)
-9974  format (17h output of iftaoh)
-      write (iout2,9993) (x(i),i=1,no4p1)
-      write (iout2,9992)
-c
-c begin a new page
-c
-      write (iout2,9973)
-9973  format (1h1)
-      go to 10
-      end
diff --git a/math/ieee/chap1/test/test17.f b/math/ieee/chap1/test/test17.f
deleted file mode 100644
index 15a34788f..000000000
--- a/math/ieee/chap1/test/test17.f
+++ /dev/null
@@ -1,147 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c main program: test program to exercise the wfta subroutine
-c   the test waveform is a complex exponential a**i whose
-c   transform is known analytically to be (1 - a**n)/(1 - a*w**k).
-c
-c authors:
-c   james h. mcclellan     and     hamid nawab
-c   department of electrical engineering and computer science
-c   massachusetts of technology
-c   cambridge, mass.  02139
-c
-c inputs:
-c   n-- transform length. it must be formed as the product of
-c       relatively prime integers from the set:
-c           2,3,4,5,7,8,9,16
-c   invrs is the flag for forward or inverse transform.
-c           invrs = 1 yields inverse transform
-c           invrs .ne. 1 gives forward transform
-c   rad and phi are the magnitude and angle (as a fraction of
-c       2*pi/n) of the complex exponential test signal.
-c          suggestion: rad = 0.98, phi = 0.5.
-c-----------------------------------------------------------------------
-c
-      double precision pi2,pin,xn,xj,xt
-      dimension xr(1260),xi(1260)
-      complex cone,ca,can,cnum,cden
-c
-c   output will be punched
-c
-      iout=i1mach(3)
-      input=i1mach(1)
-      cone=cmplx(1.0,0.0)
-      pi2=8.0d0*datan(1.0d0)
-50    continue
-      read(input,130)n
-130   format(i5)
-      write(iout,150) n
-150   format(10h length = ,i5)
-      if(n.le.0 .or. n.gt.1260) stop
-c
-c   enter a 1 to perform the inverse
-c
-c     read(input,130) invrs
-      invrs = 0
-c
-c   enter magnitude and angle (in fraction of 2*pi/n)
-c   avoid multiples of n for the angle if the radius is
-c   close to one.  suggestion: rad = 0.98, phi = 0.5.
-c
-c     read(input,160) rad,phi
-      rad = 0.98
-      phi = 0.5
-160   format(2f15.10)
-      xn=float(n)
-      pin=phi
-      pin=pin*pi2/xn
-c
-c   generate z**j
-c
-      init=0
-      do 200 j=1,n
-      an=rad**(j-1)
-      xj=j-1
-      xj=xj*pin
-      xt=dcos(xj)
-      xr(j)=xt
-      xr(j)=xr(j)*an
-      xt=dsin(xj)
-      xi(j)=xt
-      xi(j)=xi(j)*an
-200   continue
-      can=cmplx(xr(n),xi(n))
-      ca=cmplx(xr(2),xi(2))
-      can=can*ca
-c
-c   print first 50 values of input sequence
-c
-      max=50
-      if(n.lt.50)max=n
-      write(iout,300)(j,xr(j),xi(j),j=1,max)
-c
-c   call the winograd fourier transform algorithm
-c
-      call wfta(xr,xi,n,invrs,init,ierr)
-c
-c   check for error return
-c
-      if(ierr.lt.0) write(iout,250) ierr
-250   format(1x,5herror,i5)
-      if(ierr.lt.0) go to 50
-c
-c   print first 50 values of the transformed sequence
-c
-      write(iout,300)(j,xr(j),xi(j),j=1,max)
-300   format(1x,3hj =,i3,6hreal =,e20.12,6himag =,e20.12)
-c
-c   calculate absolute and relative deviations
-c
-      devabs=0.0
-      devrel=0.0
-      cnum=cone-can
-      pin=pi2/xn
-      do 350 j=1,n
-      xj=j-1
-      xj=-xj*pin
-      if(invrs.eq.1) xj=-xj
-      tr=dcos(xj)
-      ti=dsin(xj)
-      can=cmplx(tr,ti)
-      cden=cone-ca*can
-      cden=cnum/cden
-c
-c   true value of the transform (1. - a**n)/(1. - a*w**k),
-c   where a = rad*exp(j*phi*(2*pi/n)), w = exp(-j*2*pi/n).
-c   for the inverse transform the complex exponential w
-c   is conjugated.
-c
-      tr=real(cden)
-      ti=aimag(cden)
-      if(invrs.ne.1) go to 330
-c
-c   scale inverse transform by 1/n
-c
-      tr=tr/float(n)
-      ti=ti/float(n)
-330   tr=xr(j)-tr
-      ti=xi(j)-ti
-      devabs=sqrt(tr*tr+ti*ti)
-      xmag=sqrt(xr(j)*xr(j)+xi(j)*xi(j))
-      devrel=100.0*devabs/xmag
-      if(devabs.le.devmx1)go to 340
-      devmx1=devabs
-      labs=j-1
-340   if(devrel.le.devmx2)go to 350
-      devmx2=devrel
-      lrel=j-1
-350   continue
-c
-c   print the absolute and relative deviations together
-c   with their locations.
-c
-      write(iout,380) devabs,labs,devrel,lrel
-380   format(1x,21habsolute deviation = ,e20.12,9h at index,i5/
-     1 1x,21hrelative deviation = ,f11.7,8h percent,1x,9h at index,i5)
-      go to 50
-      end
diff --git a/math/ieee/chap1/test/test18.f b/math/ieee/chap1/test/test18.f
deleted file mode 100644
index cf550e015..000000000
--- a/math/ieee/chap1/test/test18.f
+++ /dev/null
@@ -1,71 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c main program: time-efficient radix-4 fast fourier transform
-c author:       l. robert morris
-c               department of systems engineering and computing science
-c               carleton university, ottawa, canada k1s 5b6
-c
-c input:        the array "a" contains the data to be transformed
-c-----------------------------------------------------------------------
-c
-c       test program for autogen radix-4 fft
-c
-        dimension a(2048),b(2048)
-        common /aa/a
-c
-        ioutd=i1mach(2)
-c
-c       compute dft and idft for n = 64, 256, and 1024 complex points
-c
-        do 1 mm=3,5
-        n=4**mm
-        do 2 j=1,n
-        a(2*j-1)=uni(0)
-        a(2*j  )=uni(0)
-        b(2*j-1)=a(2*j-1)
-2       b(2*j  )=a(2*j)
-c
-c       forward dft
-c
-        call radix4(mm,1,-1)
-c
-        if(mm.ne.3) go to 5
-c
-c       list dft input, output for n = 64 only
-c
-        write(ioutd,98)
-        write(ioutd,100)
-        do 3 j=1,n
-        write(ioutd,96) b(2*j-1),b(2*j),a(2*j-1),a(2*j)
-3       continue
-c
-c       inverse dft
-c
-5       call radix4(mm,0, 1)
-c
-c       list dft input, idft output for n = 64 only
-c
-        if(mm.ne.3) go to 7
-c
-        write(ioutd,99)
-        write(ioutd,100)
-        do 6 j=1,n
-        write(ioutd,96) b(2*j-1),b(2*j),a(2*j-1),a(2*j)
-6       continue
-c
-c       calculate rms error
-c
-7       err=0.0
-        do 8 j=1,n
-8       err=err+(a(2*j-1)-b(2*j-1))**2+(a(2*j)-b(2*j))**2
-        err=sqrt(err/float(n))
-        write(ioutd,97) mm,err
-1       continue
-c
-96      format(1x,4(f10.6,2x))
-97      format(1x,20h   rms error for m =,i2,4h is ,e14.6/)
-98      format(1x,43h       dft input                dft  output/)
-99      format(1x,43h       dft input                idft output/)
-100     format(1x,44h   real        imag       real          imag/)
-        stop
-        end
diff --git a/math/ieee/chap1/time/time12.f b/math/ieee/chap1/time/time12.f
deleted file mode 100644
index 774c2cc9e..000000000
--- a/math/ieee/chap1/time/time12.f
+++ /dev/null
@@ -1,53 +0,0 @@
-c-----------------------------------------------------------------------
-c
-c-----------------------------------------------------------------------
-c
-      parameter  (SIZE = 1024, ILOOP = 100)
-      complex  a, w
-      real     breal(SIZE), bimag(SIZE), qbreal(SIZE), qbimag(SIZE)
-c
-      ioutd = i1mach(2)
-      nn = SIZE
-      tpi = 8.*atan(1.)
-      tpion = tpi/float(nn)
-      w = cmplx(cos(tpion),-sin(tpion))
-c
-c generate a**k as test function
-c result to b(i) for modification by dft and idft subroutines and
-c a copy to qb(i) to compare final result with for error difference.
-c
-      a = (.9,.3)
-      breal(1) = 1.0
-      bimag(1) = 0.0
-      qbreal(1) = 1.0
-      qbimag(1) = 0.0
-      do 10 k=2,nn
-        w = a**(k-1)
-        breal(k) = real(w)
-        bimag(k) = aimag(w)
-        qbreal(k) = breal(k)
-        qbimag(k) = bimag(k)
-  10  continue
-c
-c now compute dft, idft, dft, idft, ...
-c first dft is computed specially, in case subroutine needs to be started.
-c
-      call fft842(0, SIZE, breal, bimag)
-      do 25 icount = 1, ILOOP
-         call fft842(1, SIZE, breal, bimag)
-         call fft842(0, SIZE, breal, bimag)
-   25 continue
-      call fft842(1, SIZE, breal, bimag)
-c
-c calculate rms error between b(i) and qb(i).
-c
-      err = 0.
-      do 30 i=1,nn
-        err = err + (breal(i)-qbreal(i))**2
-     *            + (bimag(i)-qbimag(i))**2
-  30  continue
-      err = sqrt(err / float(nn))
-      write (ioutd,9994) ILOOP, err
- 9994 format(' rms error, after ', i5, ' loops = ', e15.8)
-      stop
-      end
diff --git a/math/ieee/chap1/time/time17.f b/math/ieee/chap1/time/time17.f
deleted file mode 100644
index adcce4c2d..000000000
--- a/math/ieee/chap1/time/time17.f
+++ /dev/null
@@ -1,53 +0,0 @@
-c-----------------------------------------------------------------------
-c
-c-----------------------------------------------------------------------
-c
-      parameter  (SIZE = 1008, ILOOP = 100)
-      complex  a, w
-      real     breal(SIZE), bimag(SIZE), qbreal(SIZE), qbimag(SIZE)
-c
-      ioutd = i1mach(2)
-      nn = SIZE
-      tpi = 8.*atan(1.)
-      tpion = tpi/float(nn)
-      w = cmplx(cos(tpion),-sin(tpion))
-c
-c generate a**k as test function
-c result to b(i) for modification by dft and idft subroutines and
-c a copy to qb(i) to compare final result with for error difference.
-c
-      a = (.9,.3)
-      breal(1) = 1.0
-      bimag(1) = 0.0
-      qbreal(1) = 1.0
-      qbimag(1) = 0.0
-      do 10 k=2,nn
-        w = a**(k-1)
-        breal(k) = real(w)
-        bimag(k) = aimag(w)
-        qbreal(k) = breal(k)
-        qbimag(k) = bimag(k)
-  10  continue
-c
-c now compute dft, idft, dft, idft, ...
-c first dft is computed specially, in case subroutine needs to be started.
-c
-      call wfta(breal, bimag, SIZE, 0, 0, ierr)
-      do 25 icount = 1, ILOOP
-         call wfta(breal, bimag, SIZE, 1, 1, ierr)
-         call wfta(breal, bimag, SIZE, 0, 1, ierr)
-   25 continue
-      call wfta(breal, bimag, SIZE, 1, 1, ierr)
-c
-c calculate rms error between b(i) and qb(i).
-c
-      err = 0.
-      do 30 i=1,nn
-        err = err + (breal(i)-qbreal(i))**2
-     *            + (bimag(i)-qbimag(i))**2
-  30  continue
-      err = sqrt(err / float(nn))
-      write (ioutd,9994) ILOOP, err
- 9994 format(' rms error, after ', i5, ' loops = ', e15.8)
-      stop
-      end
diff --git a/math/ieee/chap1/time/time18.f b/math/ieee/chap1/time/time18.f
deleted file mode 100644
index 876805540..000000000
--- a/math/ieee/chap1/time/time18.f
+++ /dev/null
@@ -1,48 +0,0 @@
-c-----------------------------------------------------------------------
-c
-c-----------------------------------------------------------------------
-c
-      parameter  (SIZE = 1024, ILOOP = 100)
-      complex  w, b(SIZE), qb(SIZE), a
-      common   /aa/ b
-c
-      ioutd = i1mach(2)
-      nn = SIZE
-      tpi = 8.*atan(1.)
-      tpion = tpi/float(nn)
-      w = cmplx(cos(tpion),-sin(tpion))
-c
-c generate a**k as test function
-c result to b(i) for modification by dft and idft subroutines and
-c a copy to qb(i) to compare final result with for error difference.
-c
-      a = (.9,.3)
-      b(1) = (1.,0.)
-      qb(1) = b(1)
-      do 10 k=2,nn
-        b(k) = a**(k-1)
-        qb(k) = b(k)
-  10  continue
-c
-c now compute dft, idft, dft, idft, ...
-c first dft is computed specially, in case subroutine needs to be started.
-c
-      call radix4(5, 1, -1)
-      do 25 icount = 1, ILOOP
-         call radix4(5, 0,  1)
-         call radix4(5, 0, -1)
-   25 continue
-      call radix4(5, 0,  1)
-c
-c calculate rms error between b(i) and qb(i).
-c
-      err = 0.
-      do 30 i=1,nn
-        de = cabs(qb(i)-b(i))
-        err = err + de**2
-  30  continue
-      err = sqrt(err / float(nn))
-      write (ioutd,9994) ILOOP, err
- 9994 format(' rms error, after ', i5, ' loops = ', e15.8)
-      stop
-      end
diff --git a/math/ieee/chap1/weave1.f b/math/ieee/chap1/weave1.f
deleted file mode 100644
index 9c83d0eae..000000000
--- a/math/ieee/chap1/weave1.f
+++ /dev/null
@@ -1,371 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: weave1
-c   this subroutine implements the different pre-weave
-c   modules of the wfta.  the working arrays are sr and si.
-c   the routine checks to see which factors are present
-c   in the transform length n = na*nb*nc*nd and executes
-c   the pre-weave code for these factors.
-c
-c-----------------------------------------------------------------------
-c
-      subroutine weave1(sr,si)
-      common na,nb,nc,nd,nd1,nd2,nd3,nd4
-      dimension q(8),t(16)
-      dimension sr(1),si(1)
-      if(na.eq.1) go to 300
-      if(na.ne.2) go to 800
-c
-c **********************************************************************
-c
-c     the following code implements the 2 point pre-weave module
-c
-c **********************************************************************
-c
-      nlup2=2*(nd2-nb)
-      nlup23=2*nd2*(nd3-nc)
-      nbase=1
-      do 240 n4=1,nd
-      do 230 n3=1,nc
-      do 220 n2=1,nb
-      nr1=nbase+1
-      t0=sr(nbase)+sr(nr1)
-      sr(nr1)=sr(nbase)-sr(nr1)
-      sr(nbase)=t0
-      t0=si(nbase)+si(nr1)
-      si(nr1)=si(nbase)-si(nr1)
-      si(nbase)=t0
-220   nbase=nbase+2
-230   nbase=nbase+nlup2
-240   nbase=nbase+nlup23
-800   if(na.ne.8) go to 1600
-c
-c **********************************************************************
-c
-c     the following code implements the 8 point pre-weave module
-c
-c **********************************************************************
-c
-      nlup2=8*(nd2-nb)
-      nlup23=8*nd2*(nd3-nc)
-      nbase=1
-      do 840 n4=1,nd
-      do 830 n3=1,nc
-      do 820 n2=1,nb
-      nr1=nbase+1
-      nr2=nr1+1
-      nr3=nr2+1
-      nr4=nr3+1
-      nr5=nr4+1
-      nr6=nr5+1
-      nr7=nr6+1
-      t3=sr(nr3)+sr(nr7)
-      t7=sr(nr3)-sr(nr7)
-      t0=sr(nbase)+sr(nr4)
-      sr(nr4)=sr(nbase)-sr(nr4)
-      t1=sr(nr1)+sr(nr5)
-      t5=sr(nr1)-sr(nr5)
-      t2=sr(nr2)+sr(nr6)
-      sr(nr6)=sr(nr2)-sr(nr6)
-      sr(nbase)=t0+t2
-      sr(nr2)=t0-t2
-      sr(nr1)=t1+t3
-      sr(nr3)=t1-t3
-      sr(nr5)=t5+t7
-      sr(nr7)=t5-t7
-      t3=si(nr3)+si(nr7)
-      t7=si(nr3)-si(nr7)
-      t0=si(nbase)+si(nr4)
-      si(nr4)=si(nbase)-si(nr4)
-      t1=si(nr1)+si(nr5)
-      t5=si(nr1)-si(nr5)
-      t2=si(nr2)+si(nr6)
-      si(nr6)=si(nr2)-si(nr6)
-      si(nbase)=t0+t2
-      si(nr2)=t0-t2
-      si(nr1)=t1+t3
-      si(nr3)=t1-t3
-      si(nr5)=t5+t7
-      si(nr7)=t5-t7
-820   nbase=nbase+8
-830   nbase=nbase+nlup2
-840   nbase=nbase+nlup23
-1600  if(na.ne.16) go to 300
-c
-c **********************************************************************
-c
-c     the following code implements the 16 point pre-weave module
-c
-c **********************************************************************
-c
-      nlup2=18*(nd2-nb)
-      nlup23=18*nd2*(nd3-nc)
-      nbase=1
-      do 1640 n4=1,nd
-      do 1630 n3=1,nc
-      do 1620 n2=1,nb
-      nr1=nbase+1
-      nr2=nr1+1
-      nr3=nr2+1
-      nr4=nr3+1
-      nr5=nr4+1
-      nr6=nr5+1
-      nr7=nr6+1
-      nr8=nr7+1
-      nr9=nr8+1
-      nr10=nr9+1
-      nr11=nr10+1
-      nr12=nr11+1
-      nr13=nr12+1
-      nr14=nr13+1
-      nr15=nr14+1
-      nr16=nr15+1
-      nr17=nr16+1
-      jbase=nbase
-      do 1645 j=1,8
-      t(j)=sr(jbase)+sr(jbase+8)
-      t(j+8)=sr(jbase)-sr(jbase+8)
-      jbase=jbase+1
-1645  continue
-      do 1650 j=1,4
-      q(j)=t(j)+t(j+4)
-      q(j+4)=t(j)-t(j+4)
-1650  continue
-      sr(nbase)=q(1)+q(3)
-      sr(nr2)=q(1)-q(3)
-      sr(nr1)=q(2)+q(4)
-      sr(nr3)=q(2)-q(4)
-      sr(nr5)=q(6)+q(8)
-      sr(nr7)=q(6)-q(8)
-      sr(nr4)=q(5)
-      sr(nr6)=q(7)
-      sr(nr8)=t(9)
-      sr(nr9)=t(10)+t(16)
-      sr(nr15)=t(10)-t(16)
-      sr(nr13)=t(14)+t(12)
-      sr(nr11)=t(14)-t(12)
-      sr(nr17)=sr(nr11)+sr(nr15)
-      sr(nr16)=sr(nr9)+sr(nr13)
-      sr(nr10)=t(11)+t(15)
-      sr(nr14)=t(11)-t(15)
-      sr(nr12)=t(13)
-      jbase=nbase
-      do 1745 j=1,8
-      t(j)=si(jbase)+si(jbase+8)
-      t(j+8)=si(jbase)-si(jbase+8)
-      jbase=jbase+1
-1745  continue
-      do 1750 j=1,4
-      q(j)=t(j)+t(j+4)
-      q(j+4)=t(j)-t(j+4)
-1750  continue
-      si(nbase)=q(1)+q(3)
-      si(nr2)=q(1)-q(3)
-      si(nr1)=q(2)+q(4)
-      si(nr3)=q(2)-q(4)
-      si(nr5)=q(6)+q(8)
-      si(nr7)=q(6)-q(8)
-      si(nr4)=q(5)
-      si(nr6)=q(7)
-      si(nr8)=t(9)
-      si(nr9)=t(10)+t(16)
-      si(nr15)=t(10)-t(16)
-      si(nr13)=t(14)+t(12)
-      si(nr11)=t(14)-t(12)
-      si(nr17)=si(nr11)+si(nr15)
-      si(nr16)=si(nr9)+si(nr13)
-      si(nr10)=t(11)+t(15)
-      si(nr14)=t(11)-t(15)
-      si(nr12)=t(13)
-1620  nbase=nbase+18
-1630  nbase=nbase+nlup2
-1640  nbase=nbase+nlup23
-300   if(nb.eq.1) go to 700
-      if(nb.ne.3) go to 900
-c
-c **********************************************************************
-c
-c     the following code implements the 3 point pre-weave module
-c
-c **********************************************************************
-c
-      nlup2=2*nd1
-      nlup23=3*nd1*(nd3-nc)
-      nbase=1
-      noff=nd1
-      do 340 n4=1,nd
-      do 330 n3=1,nc
-      do 310 n2=1,nd1
-      nr1=nbase+noff
-      nr2=nr1+noff
-      t1=sr(nr1)+sr(nr2)
-      sr(nbase)=sr(nbase)+t1
-      sr(nr2)=sr(nr1)-sr(nr2)
-      sr(nr1)=t1
-      t1=si(nr1)+si(nr2)
-      si(nbase)=si(nbase)+t1
-      si(nr2)=si(nr1)-si(nr2)
-      si(nr1)=t1
-310   nbase=nbase+1
-330   nbase=nbase+nlup2
-340   nbase=nbase+nlup23
-900   if(nb.ne.9) go to 700
-c
-c **********************************************************************
-c
-c     the following code implements the 9 point pre-weave module
-c
-c **********************************************************************
-c
-      nlup2=10*nd1
-      nlup23=11*nd1*(nd3-nc)
-      nbase=1
-      noff=nd1
-      do 940 n4=1,nd
-      do 930 n3=1,nc
-      do 910 n2=1,nd1
-      nr1=nbase+noff
-      nr2=nr1+noff
-      nr3=nr2+noff
-      nr4=nr3+noff
-      nr5=nr4+noff
-      nr6=nr5+noff
-      nr7=nr6+noff
-      nr8=nr7+noff
-      nr9=nr8+noff
-      nr10=nr9+noff
-      t3=sr(nr3)+sr(nr6)
-      t6=sr(nr3)-sr(nr6)
-      sr(nbase)=sr(nbase)+t3
-      t7=sr(nr7)+sr(nr2)
-      t2=sr(nr7)-sr(nr2)
-      sr(nr2)=t6
-      t1=sr(nr1)+sr(nr8)
-      t8=sr(nr1)-sr(nr8)
-      sr(nr1)=t3
-      t4=sr(nr4)+sr(nr5)
-      t5=sr(nr4)-sr(nr5)
-      sr(nr3)=t1+t4+t7
-      sr(nr4)=t1-t7
-      sr(nr5)=t4-t1
-      sr(nr6)=t7-t4
-      sr(nr10)=t2+t5+t8
-      sr(nr7)=t8-t2
-      sr(nr8)=t5-t8
-      sr(nr9)=t2-t5
-      t3=si(nr3)+si(nr6)
-      t6=si(nr3)-si(nr6)
-      si(nbase)=si(nbase)+t3
-      t7=si(nr7)+si(nr2)
-      t2=si(nr7)-si(nr2)
-      si(nr2)=t6
-      t1=si(nr1)+si(nr8)
-      t8=si(nr1)-si(nr8)
-      si(nr1)=t3
-      t4=si(nr4)+si(nr5)
-      t5=si(nr4)-si(nr5)
-      si(nr3)=t1+t4+t7
-      si(nr4)=t1-t7
-      si(nr5)=t4-t1
-      si(nr6)=t7-t4
-      si(nr10)=t2+t5+t8
-      si(nr7)=t8-t2
-      si(nr8)=t5-t8
-      si(nr9)=t2-t5
-910   nbase=nbase+1
-930   nbase=nbase+nlup2
-940   nbase=nbase+nlup23
-700   if(nc.ne.7) go to 500
-c
-c **********************************************************************
-c
-c     the following code implements the 7 point pre-weave module
-c
-c **********************************************************************
-c
-      noff=nd1*nd2
-      nbase=1
-      nlup2=8*noff
-      do 740 n4=1,nd
-      do 710 n1=1,noff
-      nr1=nbase+noff
-      nr2=nr1+noff
-      nr3=nr2+noff
-      nr4=nr3+noff
-      nr5=nr4+noff
-      nr6=nr5+noff
-      nr7=nr6+noff
-      nr8=nr7+noff
-      t1=sr(nr1)+sr(nr6)
-      t6=sr(nr1)-sr(nr6)
-      t4=sr(nr4)+sr(nr3)
-      t3=sr(nr4)-sr(nr3)
-      t2=sr(nr2)+sr(nr5)
-      t5=sr(nr2)-sr(nr5)
-      sr(nr5)=t6-t3
-      sr(nr2)=t5+t3+t6
-      sr(nr6)=t5-t6
-      sr(nr8)=t3-t5
-      sr(nr3)=t2-t1
-      sr(nr4)=t1-t4
-      sr(nr7)=t4-t2
-      t1=t1+t4+t2
-      sr(nbase)=sr(nbase)+t1
-      sr(nr1)=t1
-      t1=si(nr1)+si(nr6)
-      t6=si(nr1)-si(nr6)
-      t4=si(nr4)+si(nr3)
-      t3=si(nr4)-si(nr3)
-      t2=si(nr2)+si(nr5)
-      t5=si(nr2)-si(nr5)
-      si(nr5)=t6-t3
-      si(nr2)=t5+t3+t6
-      si(nr6)=t5-t6
-      si(nr8)=t3-t5
-      si(nr3)=t2-t1
-      si(nr4)=t1-t4
-      si(nr7)=t4-t2
-      t1=t1+t4+t2
-      si(nbase)=si(nbase)+t1
-      si(nr1)=t1
-710   nbase=nbase+1
-740   nbase=nbase+nlup2
-500   if(nd.ne.5) return
-c
-c **********************************************************************
-c
-c     the following code implements the 5 point pre-weave module
-c
-c **********************************************************************
-c
-      noff=nd1*nd2*nd3
-      nbase=1
-      do 510 n1=1,noff
-      nr1=nbase+noff
-      nr2=nr1+noff
-      nr3=nr2+noff
-      nr4=nr3+noff
-      nr5=nr4+noff
-      t4=sr(nr1)-sr(nr4)
-      t1=sr(nr1)+sr(nr4)
-      t3=sr(nr3)+sr(nr2)
-      t2=sr(nr3)-sr(nr2)
-      sr(nr3)=t1-t3
-      sr(nr1)=t1+t3
-      sr(nbase)=sr(nbase)+sr(nr1)
-      sr(nr5)=t2+t4
-      sr(nr2)=t4
-      sr(nr4)=t2
-      t4=si(nr1)-si(nr4)
-      t1=si(nr1)+si(nr4)
-      t3=si(nr3)+si(nr2)
-      t2=si(nr3)-si(nr2)
-      si(nr3)=t1-t3
-      si(nr1)=t1+t3
-      si(nbase)=si(nbase)+si(nr1)
-      si(nr5)=t2+t4
-      si(nr2)=t4
-      si(nr4)=t2
-510   nbase=nbase+1
-      return
-      end
diff --git a/math/ieee/chap1/weave2.f b/math/ieee/chap1/weave2.f
deleted file mode 100644
index 5c04ff914..000000000
--- a/math/ieee/chap1/weave2.f
+++ /dev/null
@@ -1,412 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c
-c subroutine: weave2
-c   this subroutine implements the post-weave modules
-c   of the wfta.  the working arrays are sr and si.
-c   the routine checks to see which factors are present
-c   in the transform length n = na*nb*nc*nd and executes
-c   the post-weave code for these factors.
-c
-c-----------------------------------------------------------------------
-c
-      subroutine weave2(sr,si)
-      common na,nb,nc,nd,nd1,nd2,nd3,nd4
-      dimension sr(1),si(1)
-      dimension q(8),t(16)
-      if(nd.ne.5) go to 700
-c
-c **********************************************************************
-c
-c     the following code implements the 5 point post-weave module
-c
-c **********************************************************************
-c
-      noff=nd1*nd2*nd3
-      nbase=1
-      do 510 n1=1,noff
-      nr1=nbase+noff
-      nr2=nr1+noff
-      nr3=nr2+noff
-      nr4=nr3+noff
-      nr5=nr4+noff
-      t1=sr(nbase)+sr(nr1)
-      t3=t1-sr(nr3)
-      t1=t1+sr(nr3)
-      t4=si(nr2)+si(nr5)
-      t2=si(nr4)+si(nr5)
-      sr(nr1)=t1-t4
-      sr4=t1+t4
-      sr2=t3+t2
-      sr(nr3)=t3-t2
-      t1=si(nbase)+si(nr1)
-      t3=t1-si(nr3)
-      t1=t1+si(nr3)
-      t4=sr(nr2)+sr(nr5)
-      t2=sr(nr4)+sr(nr5)
-      si(nr1)=t1+t4
-      si(nr4)=t1-t4
-      si(nr2)=t3-t2
-      si(nr3)=t3+t2
-      sr(nr2)=sr2
-      sr(nr4)=sr4
-510   nbase=nbase+1
-700   if(nc.ne.7) go to 300
-c
-c **********************************************************************
-c
-c     the following code implements the 7 point post-weave module
-c
-c **********************************************************************
-c
-      noff=nd1*nd2
-      nbase=1
-      nlup2=8*noff
-      do 740 n4=1,nd
-      do 710 n1=1,noff
-      nr1=nbase+noff
-      nr2=nr1+noff
-      nr3=nr2+noff
-      nr4=nr3+noff
-      nr5=nr4+noff
-      nr6=nr5+noff
-      nr7=nr6+noff
-      nr8=nr7+noff
-      t1=sr(nr1)+sr(nbase)
-      t2=t1-sr(nr3)-sr(nr4)
-      t4=t1+sr(nr3)-sr(nr7)
-      t1=t1+sr(nr4)+sr(nr7)
-      t6=si(nr2)+si(nr5)+si(nr8)
-      t5=si(nr2)-si(nr5)-si(nr6)
-      t3=si(nr2)+si(nr6)-si(nr8)
-      sr(nr1)=t1-t6
-      sr6=t1+t6
-      sr2=t2-t5
-      sr5=t2+t5
-      sr(nr4)=t4-t3
-      sr(nr3)=t4+t3
-      t1=si(nr1)+si(nbase)
-      t2=t1-si(nr3)-si(nr4)
-      t4=t1+si(nr3)-si(nr7)
-      t1=t1+si(nr4)+si(nr7)
-      t6=sr(nr2)+sr(nr5)+sr(nr8)
-      t5=sr(nr2)-sr(nr5)-sr(nr6)
-      t3=sr(nr2)+sr(nr6)-sr(nr8)
-      si(nr1)=t1+t6
-      si(nr6)=t1-t6
-      si(nr2)=t2+t5
-      si(nr5)=t2-t5
-      si(nr4)=t4+t3
-      si(nr3)=t4-t3
-      sr(nr2)=sr2
-      sr(nr5)=sr5
-      sr(nr6)=sr6
-710   nbase=nbase+1
-740   nbase=nbase+nlup2
-300   if(nb.eq.1) go to 400
-      if(nb.ne.3) go to 900
-c
-c **********************************************************************
-c
-c     the following code implements the 3 point post-weave module
-c
-c **********************************************************************
-c
-      nlup2=2*nd1
-      nlup23=3*nd1*(nd3-nc)
-      nbase=1
-      noff=nd1
-      do 340 n5=1,nd
-      do 330 n4=1,nc
-      do 310 n2=1,nd1
-      nr1=nbase+noff
-      nr2=nr1+noff
-      t1=sr(nbase)+sr(nr1)
-      sr(nr1)=t1-si(nr2)
-      sr2=t1+si(nr2)
-      t1=si(nbase)+si(nr1)
-      si(nr1)=t1+sr(nr2)
-      si(nr2)=t1-sr(nr2)
-      sr(nr2)=sr2
-310   nbase=nbase+1
-330   nbase=nbase+nlup2
-340   nbase=nbase+nlup23
-900   if(nb.ne.9) go to 400
-c
-c **********************************************************************
-c
-c     the following code implements the 9 point post-weave module
-c
-c **********************************************************************
-c
-      nlup2=10*nd1
-      nlup23=11*nd1*(nd3-nc)
-      nbase=1
-      noff=nd1
-      do 940 n4=1,nd
-      do 930 n3=1,nc
-      do 910 n2=1,nd1
-      nr1=nbase+noff
-      nr2=nr1+noff
-      nr3=nr2+noff
-      nr4=nr3+noff
-      nr5=nr4+noff
-      nr6=nr5+noff
-      nr7=nr6+noff
-      nr8=nr7+noff
-      nr9=nr8+noff
-      nr10=nr9+noff
-      t3=sr(nbase)-sr(nr3)
-      t7=sr(nbase)+sr(nr1)
-      sr(nbase)=sr(nbase)+sr(nr3)+sr(nr3)
-      t6=t3+si(nr10)
-      sr(nr3)=t3-si(nr10)
-      t4=t7+sr(nr5)-sr(nr6)
-      t1=t7-sr(nr4)-sr(nr5)
-      t7=t7+sr(nr4)+sr(nr6)
-      sr(nr6)=t6
-      t8=si(nr2)-si(nr7)-si(nr8)
-      t5=si(nr2)+si(nr8)-si(nr9)
-      t2=si(nr2)+si(nr7)+si(nr9)
-      sr(nr1)=t7-t2
-      sr8=t7+t2
-      sr(nr4)=t1-t8
-      sr(nr5)=t1+t8
-      sr7=t4-t5
-      sr2=t4+t5
-      t3=si(nbase)-si(nr3)
-      t7=si(nbase)+si(nr1)
-      si(nbase)=si(nbase)+si(nr3)+si(nr3)
-      t6=t3-sr(nr10)
-      si(nr3)=t3+sr(nr10)
-      t4=t7+si(nr5)-si(nr6)
-      t1=t7-si(nr4)-si(nr5)
-      t7=t7+si(nr4)+si(nr6)
-      si(nr6)=t6
-      t8=sr(nr2)-sr(nr7)-sr(nr8)
-      t5=sr(nr2)+sr(nr8)-sr(nr9)
-      t2=sr(nr2)+sr(nr7)+sr(nr9)
-      si(nr1)=t7+t2
-      si(nr8)=t7-t2
-      si(nr4)=t1+t8
-      si(nr5)=t1-t8
-      si(nr7)=t4+t5
-      si(nr2)=t4-t5
-      sr(nr2)=sr2
-      sr(nr7)=sr7
-      sr(nr8)=sr8
-910   nbase=nbase+1
-930   nbase=nbase+nlup2
-940   nbase=nbase+nlup23
-400   if(na.eq.1) return
-      if(na.ne.4) go to 800
-c
-c **********************************************************************
-c
-c     the following code implements the 4 point post-weave module
-c
-c **********************************************************************
-c
-      nlup2=4*(nd2-nb)
-      nlup23=4*nd2*(nd3-nc)
-      nbase=1
-      do 440 n4=1,nd
-      do 430 n3=1,nc
-      do 420 n2=1,nb
-      nr1=nbase+1
-      nr2=nr1+1
-      nr3=nr2+1
-      tr0=sr(nbase)+sr(nr2)
-      tr2=sr(nbase)-sr(nr2)
-      tr1=sr(nr1)+sr(nr3)
-      tr3=sr(nr1)-sr(nr3)
-      ti1=si(nr1)+si(nr3)
-      ti3=si(nr1)-si(nr3)
-      sr(nbase)=tr0+tr1
-      sr(nr2)=tr0-tr1
-      sr(nr1)=tr2+ti3
-      sr(nr3)=tr2-ti3
-      ti0=si(nbase)+si(nr2)
-      ti2=si(nbase)-si(nr2)
-      si(nbase)=ti0+ti1
-      si(nr2)=ti0-ti1
-      si(nr1)=ti2-tr3
-      si(nr3)=ti2+tr3
-420   nbase=nbase+4
-430   nbase=nbase+nlup2
-440   nbase=nbase+nlup23
-800   if(na.ne.8) go to 1600
-c
-c **********************************************************************
-c
-c     the following code implements the 8 point post-weave module
-c
-c **********************************************************************
-c
-      nlup2=8*(nd2-nb)
-      nlup23=8*nd2*(nd3-nc)
-      nbase=1
-      do 840 n4=1,nd
-      do 830 n3=1,nc
-      do 820 n2=1,nb
-      nr1=nbase+1
-      nr2=nr1+1
-      nr3=nr2+1
-      nr4=nr3+1
-      nr5=nr4+1
-      nr6=nr5+1
-      nr7=nr6+1
-      t1=sr(nbase)-sr(nr1)
-      sr(nbase)=sr(nbase)+sr(nr1)
-      sr6=sr(nr2)+si(nr3)
-      sr(nr2)=sr(nr2)-si(nr3)
-      t4=sr(nr4)-si(nr5)
-      t5=sr(nr4)+si(nr5)
-      t6=sr(nr7)-si(nr6)
-      t7=sr(nr7)+si(nr6)
-      sr(nr4)=t1
-      sr(nr1)=t4+t6
-      sr3=t4-t6
-      sr5=t5-t7
-      sr(nr7)=t5+t7
-      t1=si(nbase)-si(nr1)
-      si(nbase)=si(nbase)+si(nr1)
-      t3=si(nr2)-sr(nr3)
-      si(nr2)=si(nr2)+sr(nr3)
-      t4=si(nr4)+sr(nr5)
-      t5=si(nr4)-sr(nr5)
-      si(nr6)=t3
-      t6=sr(nr6)+si(nr7)
-      t7=sr(nr6)-si(nr7)
-      si(nr4)=t1
-      si(nr1)=t4+t6
-      si(nr3)=t4-t6
-      si(nr5)=t5+t7
-      si(nr7)=t5-t7
-      sr(nr3)=sr3
-      sr(nr5)=sr5
-      sr(nr6)=sr6
-820   nbase=nbase+8
-830   nbase=nbase+nlup2
-840   nbase=nbase+nlup23
-1600  if(na.ne.16) return
-c
-c **********************************************************************
-c
-c     the following code implements the 16 point post-weave module
-c
-c **********************************************************************
-c
-      nlup2=18*(nd2-nb)
-      nlup23=18*nd2*(nd3-nc)
-      nbase=1
-      do 1640 n4=1,nd
-      do 1630 n3=1,nc
-      do 1620 n2=1,nb
-      nr1=nbase+1
-      nr2=nr1+1
-      nr3=nr2+1
-      nr4=nr3+1
-      nr5=nr4+1
-      nr6=nr5+1
-      nr7=nr6+1
-      nr8=nr7+1
-      nr9=nr8+1
-      nr10=nr9+1
-      nr11=nr10+1
-      nr12=nr11+1
-      nr13=nr12+1
-      nr14=nr13+1
-      nr15=nr14+1
-      nr16=nr15+1
-      nr17=nr16+1
-      t(2)=sr(nbase)-sr(nr1)
-      sr(nbase)=sr(nr1)+sr(nbase)
-      t(4)=sr(nr2)+si(nr3)
-      t(3)=sr(nr2)-si(nr3)
-      t(6)=sr(nr4)+si(nr5)
-      t(5)=sr(nr4)-si(nr5)
-      t(8)=-si(nr6)-sr(nr7)
-      t(7)=-si(nr6)+sr(nr7)
-      t(9)=sr(nr8)+sr(nr14)
-      t(15)=sr(nr8)-sr(nr14)
-      t(13)=-si(nr10)-si(nr12)
-      t(11)=si(nr10)-si(nr12)
-      t(16)=sr(nr15)-sr(nr17)
-      t(12)=sr(nr11)-sr(nr17)
-      t(10)=-si(nr9)-si(nr16)
-      t(14)=-si(nr16)+si(nr13)
-      sr(nr2)=t(5)+t(7)
-      sr6=t(5)-t(7)
-      sr10=t(6)+t(8)
-      sr(nr14)=t(6)-t(8)
-      q(7)=t(9)+t(10)
-      q(8)=t(9)-t(10)
-      q(1)=t(11)+t(12)
-      q(2)=t(11)-t(12)
-      q(4)=t(14)+t(15)
-      q(5)=t(15)-t(14)
-      q(3)=t(13)+t(16)
-      q(6)=t(13)-t(16)
-      sr(nr1)=q(3)+q(7)
-      sr(nr7)=q(7)-q(3)
-      sr9=q(8)+q(6)
-      sr(nr15)=q(8)-q(6)
-      sr5=q(1)+q(4)
-      sr3=q(4)-q(1)
-      sr13=q(2)+q(5)
-      sr11=q(5)-q(2)
-      sr(nr8)=t(2)
-      sr(nr4)=t(3)
-      sr12=t(4)
-      t(2)=si(nbase)-si(nr1)
-      si(nbase)=si(nr1)+si(nbase)
-      t(4)=si(nr2)-sr(nr3)
-      t(3)=si(nr2)+sr(nr3)
-      t(6)=si(nr4)-sr(nr5)
-      t(5)=si(nr4)+sr(nr5)
-      t(8)=sr(nr6)-si(nr7)
-      t(7)=sr(nr6)+si(nr7)
-      t(9)=si(nr8)+si(nr14)
-      t(15)=si(nr8)-si(nr14)
-      t(13)=sr(nr10)+sr(nr12)
-      t(11)=sr(nr12)-sr(nr10)
-      t(16)=si(nr15)-si(nr17)
-      t(12)=si(nr11)-si(nr17)
-      t(10)=sr(nr9)+sr(nr16)
-      si(nr2)=t(5)+t(7)
-      si(nr6)=t(5)-t(7)
-      si(nr10)=t(6)+t(8)
-      si(nr14)=t(6)-t(8)
-      q(7)=t(9)+t(10)
-      q(8)=t(9)-t(10)
-      q(1)=t(11)+t(12)
-      q(2)=t(11)-t(12)
-      q(4)=t(14)+t(15)
-      q(5)=t(15)-t(14)
-      q(3)=t(13)+t(16)
-      q(6)=t(13)-t(16)
-      si(nr1)=q(3)+q(7)
-      si(nr7)=q(7)-q(3)
-      si(nr9)=q(8)+q(6)
-      si(nr15)=q(8)-q(6)
-      si(nr5)=q(1)+q(4)
-      si(nr3)=q(4)-q(1)
-      si(nr13)=q(2)+q(5)
-      si(nr11)=q(5)-q(2)
-      si(nr8)=t(2)
-      si(nr4)=t(3)
-      si(nr12)=t(4)
-      sr(nr3)=sr3
-      sr(nr5)=sr5
-      sr(nr6)=sr6
-      sr(nr9)=sr9
-      sr(nr10)=sr10
-      sr(nr11)=sr11
-      sr(nr12)=sr12
-      sr(nr13)=sr13
-1620  nbase=nbase+18
-1630  nbase=nbase+nlup2
-1640  nbase=nbase+nlup23
-      return
-      end
diff --git a/math/ieee/chap1/wfta.f b/math/ieee/chap1/wfta.f
deleted file mode 100644
index 9e819941a..000000000
--- a/math/ieee/chap1/wfta.f
+++ /dev/null
@@ -1,150 +0,0 @@
-c
-c-----------------------------------------------------------------------
-c subroutine: wfta
-c winograd fourier transform algorithm
-c-----------------------------------------------------------------------
-c
-      subroutine wfta(xr,xi,n,invrs,init,ierr)
-      dimension xr(1),xi(1)
-c
-c inputs:
-c   n-- transform length.  must be formed as the product of
-c       relatively prime integers from the set:
-c           2,3,4,5,7,8,9,16
-c       thus the largest possible value of n is 5040.
-c   xr(.)-- array that holds the real part of the data
-c           to be transformed.
-c   xi(.)-- array that holds the imaginary part of the
-c           data to be transformed.
-c   invrs-- parameter that flags whether or not the inverse
-c           transform is to be calculated.  a division by n
-c           is included in the inverse.
-c           invrs = 1 yields inverse transform
-c           invrs .ne. 1 gives forward transform
-c   init-- parameter that flags whether or not the program
-c          is to be initialized for this value of n.  the
-c          initialization is performed only once in order to
-c          to speed up the computation on succeeding calls
-c          to the wfta routine, when n is held fixed.
-c          init = 0 results in initialization.
-c   ierr-- error code that is negative when the wfta
-c          terminates incorrectly.
-c           0 = successful completion
-c          -1 = this value of n does not factor properly
-c          -2 = an initialization has not been done for
-c               this value of n.
-c
-c
-c   the following two cards may be changed if the maximum
-c   desired transform length is less than 5040
-c
-c  *********************************************************************
-      dimension sr(1782),si(1782),coef(1782)
-      integer indx1(1008),indx2(1008)
-c  *********************************************************************
-c
-      common na,nb,nc,nd,nd1,nd2,nd3,nd4
-c
-c   test for initial run
-c
-      if(init.eq.0) call inishl(n,coef,xr,xi,indx1,indx2,ierr)
-c
-      if(ierr.lt.0) return
-      m=na*nb*nc*nd
-      if(m.eq.n) go to 100
-      ierr=-2
-      return
-c
-c  error(-2)-- program not initialized for this value of n
-c
-100   nmult=nd1*nd2*nd3*nd4
-c
-c   the following code maps the data arrays xr and xi to
-c   the working arrays sr and si via the mapping vector
-c   indx1(.).  the permutation of the data follows the
-c   sino correspondence of the chinese remainder theorem.
-c
-      j=1
-      k=1
-      inc1=nd1-na
-      inc2=nd1*(nd2-nb)
-      inc3=nd1*nd2*(nd3-nc)
-      do 140 n4=1,nd
-      do 130 n3=1,nc
-      do 120 n2=1,nb
-      do 110 n1=1,na
-      ind=indx1(k)
-      sr(j)=xr(ind)
-      si(j)=xi(ind)
-      k=k+1
-110   j=j+1
-120   j=j+inc1
-130   j=j+inc2
-140   j=j+inc3
-c
-c   do the pre-weave modules
-c
-      call weave1(sr,si)
-c
-c   the following loop performs all the multiplications of the
-c   winograd fourier transform algorithm.  the multiplication
-c   coefficients are stored on the initialization pass in the
-c   array coef(.).
-c
-      do 200 j=1,nmult
-      sr(j)=sr(j)*coef(j)
-      si(j)=si(j)*coef(j)
-  200  continue
-c
-c   do the post-weave modules
-c
-      call weave2(sr,si)
-c
-c
-c   the following code maps the working arrays sr and si
-c   to the data arrays xr and xi via the mapping vector
-c   indx2(.).  the permutation of the data follows the
-c   chinese remainder theorem.
-c
-      j=1
-      k=1
-      inc1=nd1-na
-      inc2=nd1*(nd2-nb)
-      inc3=nd1*nd2*(nd3-nc)
-c
-c   check for inverse
-c
-      if(invrs.eq.1) go to 400
-      do 340 n4=1,nd
-      do 330 n3=1,nc
-      do 320 n2=1,nb
-      do 310 n1=1,na
-      kndx=indx2(k)
-      xr(kndx)=sr(j)
-      xi(kndx)=si(j)
-      k=k+1
-310   j=j+1
-320   j=j+inc1
-330   j=j+inc2
-340   j=j+inc3
-      return
-c
-c   different permutation for the inverse
-c
-400   fn=float(n)
-      np2=n+2
-      indx2(1)=n+1
-      do 440 n4=1,nd
-      do 430 n3=1,nc
-      do 420 n2=1,nb
-      do 410 n1=1,na
-      kndx=np2-indx2(k)
-      xr(kndx)=sr(j)/fn
-      xi(kndx)=si(j)/fn
-      k=k+1
-410   j=j+1
-420   j=j+inc1
-430   j=j+inc2
-440   j=j+inc3
-      return
-      end
diff --git a/math/ieee/d1mach.f b/math/ieee/d1mach.f
deleted file mode 100644
index 699468ab9..000000000
--- a/math/ieee/d1mach.f
+++ /dev/null
@@ -1,256 +0,0 @@
-c
-c----------------------------------------------------------------------
-c  function:  d1mach
-c  this routine is from the port mathematical subroutine library
-c  it is described in the bell laboratories computing science
-c  technical report #47 by p.a. fox, a.d. hall and n.l. schryer
-c  a modification to the "i out of bounds" error message
-c  has been made by c. a. mcgonegal - april, 1978
-c----------------------------------------------------------------------
-c
-      double precision function d1mach(i)
-c
-c  double-precision machine constants
-c
-c  d1mach( 1) = b**(emin-1), the smallest positive magnitude.
-c
-c  d1mach( 2) = b**emax*(1 - b**(-t)), the largest magnitude.
-c
-c  d1mach( 3) = b**(-t), the smallest relative spacing.
-c
-c  d1mach( 4) = b**(1-t), the largest relative spacing.
-c
-c  d1mach( 5) = log10(b)
-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.
-c
-c  where possible, octal or hexadecimal constants have been used
-c  to specify the constants exactly which has in some cases
-c  required the use of equivalent integer arrays.
-c
-      integer small(4)
-      integer large(4)
-      integer right(4)
-      integer diver(4)
-      integer log10(4)
-c
-      double precision dmach(5)
-c
-      equivalence (dmach(1),small(1))
-      equivalence (dmach(2),large(1))
-      equivalence (dmach(3),right(1))
-      equivalence (dmach(4),diver(1))
-      equivalence (dmach(5),log10(1))
-c
-c     machine constants for the burroughs 1700 system.
-c
-c     data small(1) / zc00800000 /
-c     data small(2) / z000000000 /
-c
-c     data large(1) / zdffffffff /
-c     data large(2) / zfffffffff /
-c
-c     data right(1) / zcc5800000 /
-c     data right(2) / z000000000 /
-c
-c     data diver(1) / zcc6800000 /
-c     data diver(2) / z000000000 /
-c
-c     data log10(1) / zd00e730e7 /
-c     data log10(2) / zc77800dc0 /
-c
-c     machine constants for the burroughs 5700 system.
-c
-c     data small(1) / o1771000000000000 /
-c     data small(2) / o0000000000000000 /
-c
-c     data large(1) / o0777777777777777 /
-c     data large(2) / o0007777777777777 /
-c
-c     data right(1) / o1461000000000000 /
-c     data right(2) / o0000000000000000 /
-c
-c     data diver(1) / o1451000000000000 /
-c     data diver(2) / o0000000000000000 /
-c
-c     data log10(1) / o1157163034761674 /
-c     data log10(2) / o0006677466732724 /
-c
-c     machine constants for the burroughs 6700/7700 systems.
-c
-c     data small(1) / o1771000000000000 /
-c     data small(2) / o7770000000000000 /
-c
-c     data large(1) / o0777777777777777 /
-c     data large(2) / o7777777777777777 /
-c
-c     data right(1) / o1461000000000000 /
-c     data right(2) / o0000000000000000 /
-c
-c     data diver(1) / o1451000000000000 /
-c     data diver(2) / o0000000000000000 /
-c
-c     data log10(1) / o1157163034761674 /
-c     data log10(2) / o0006677466732724 /
-c
-c     machine constants for the cdc 6000/7000 series.
-c
-c     data small(1) / 00604000000000000000b /
-c     data small(2) / 00000000000000000000b /
-c
-c     data large(1) / 37767777777777777777b /
-c     data large(2) / 37167777777777777777b /
-c
-c     data right(1) / 15604000000000000000b /
-c     data right(2) / 15000000000000000000b /
-c
-c     data diver(1) / 15614000000000000000b /
-c     data diver(2) / 15010000000000000000b /
-c
-c     data log10(1) / 17164642023241175717b /
-c     data log10(2) / 16367571421742254654b /
-c
-c     machine constants for the cray 1
-c
-c     data small(1) / 200004000000000000000b /
-c     data small(2) / 000000000000000000000b /
-c
-c     data large(1) / 577767777777777777777b /
-c     data large(2) / 000007777777777777776b /
-c
-c     data right(1) / 376424000000000000000b /
-c     data right(2) / 000000000000000000000b /
-c
-c     data diver(1) / 376434000000000000000b /
-c     data diver(2) / 000000000000000000000b /
-c
-c     data log10(1) / 377774642023241175717b /
-c     data log10(2) / 000007571421742254654b /
-c
-c     machine constants for the data general eclipse s/200
-c
-c     note - it may be appropriate to include the following card -
-c     static dmach(5)
-c
-c     data small/20k,3*0/,large/77777k,3*177777k/
-c     data right/31420k,3*0/,diver/32020k,3*0/
-c     data log10/40423k,42023k,50237k,74776k/
-c
-c     machine constants for the harris slash 6 and slash 7
-c
-c     data small(1),small(2) / '20000000, '00000201 /
-c     data large(1),large(2) / '37777777, '37777577 /
-c     data right(1),right(2) / '20000000, '00000333 /
-c     data diver(1),diver(2) / '20000000, '00000334 /
-c     data log10(1),log10(2) / '23210115, '10237777 /
-c
-c     machine constants for the honeywell 600/6000 series.
-c
-c     data small(1),small(2) / o402400000000, o000000000000 /
-c     data large(1),large(2) / o376777777777, o777777777777 /
-c     data right(1),right(2) / o604400000000, o000000000000 /
-c     data diver(1),diver(2) / o606400000000, o000000000000 /
-c     data log10(1),log10(2) / o776464202324, o117571775714 /
-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 small(1),small(2) / z00100000, z00000000 /
-c     data large(1),large(2) / z7fffffff, zffffffff /
-c     data right(1),right(2) / z33100000, z00000000 /
-c     data diver(1),diver(2) / z34100000, z00000000 /
-c     data log10(1),log10(2) / z41134413, z509f79ff /
-c
-c     machine constants for the pdp-10 (ka processor).
-c
-c     data small(1),small(2) / "033400000000, "000000000000 /
-c     data large(1),large(2) / "377777777777, "344777777777 /
-c     data right(1),right(2) / "113400000000, "000000000000 /
-c     data diver(1),diver(2) / "114400000000, "000000000000 /
-c     data log10(1),log10(2) / "177464202324, "144117571776 /
-c
-c     machine constants for the pdp-10 (ki processor).
-c
-c     data small(1),small(2) / "000400000000, "000000000000 /
-c     data large(1),large(2) / "377777777777, "377777777777 /
-c     data right(1),right(2) / "103400000000, "000000000000 /
-c     data diver(1),diver(2) / "104400000000, "000000000000 /
-c     data log10(1),log10(2) / "177464202324, "476747767461 /
-c
-c     machine constants for pdp-11 fortran's supporting
-c     32-bit integers (expressed in integer and octal).
-c
-      data small(1),small(2) /    8388608,           0 /
-      data large(1),large(2) / 2147483647,          -1 /
-      data right(1),right(2) /  612368384,           0 /
-      data diver(1),diver(2) /  620756992,           0 /
-      data log10(1),log10(2) / 1067065498, -2063872008 /
-c
-c     data small(1),small(2) / o00040000000, o00000000000 /
-c     data large(1),large(2) / o17777777777, o37777777777 /
-c     data right(1),right(2) / o04440000000, o00000000000 /
-c     data diver(1),diver(2) / o04500000000, o00000000000 /
-c     data log10(1),log10(2) / o07746420232, o20476747770 /
-c
-c     machine constants for pdp-11 fortran's supporting
-c     16-bit integers (expressed in integer and octal).
-c
-c     data small(1),small(2) /    128,      0 /
-c     data small(3),small(4) /      0,      0 /
-c
-c     data large(1),large(2) /  32767,     -1 /
-c     data large(3),large(4) /     -1,     -1 /
-c
-c     data right(1),right(2) /   9344,      0 /
-c     data right(3),right(4) /      0,      0 /
-c
-c     data diver(1),diver(2) /   9472,      0 /
-c     data diver(3),diver(4) /      0,      0 /
-c
-c     data log10(1),log10(2) /  16282,   8346 /
-c     data log10(3),log10(4) / -31493, -12296 /
-c
-c     data small(1),small(2) / o000200, o000000 /
-c     data small(3),small(4) / o000000, o000000 /
-c
-c     data large(1),large(2) / o077777, o177777 /
-c     data large(3),large(4) / o177777, o177777 /
-c
-c     data right(1),right(2) / o022200, o000000 /
-c     data right(3),right(4) / o000000, o000000 /
-c
-c     data diver(1),diver(2) / o022400, o000000 /
-c     data diver(3),diver(4) / o000000, o000000 /
-c
-c     data log10(1),log10(2) / o037632, o020232 /
-c     data log10(3),log10(4) / o102373, o147770 /
-c
-c     machine constants for the univac 1100 series.
-c
-c     data small(1),small(2) / o000040000000, o000000000000 /
-c     data large(1),large(2) / o377777777777, o777777777777 /
-c     data right(1),right(2) / o170540000000, o000000000000 /
-c     data diver(1),diver(2) / o170640000000, o000000000000 /
-c     data log10(1),log10(2) / o177746420232, o411757177572 /
-c
-c     machine constants for the vax-11 with
-c     fortran iv-plus compiler
-c
-c     data small(1),small(2) / z00000080, z00000000 /
-c     data large(1),large(2) / zffff7fff, zffffffff /
-c     data right(1),right(2) / z00002480, z00000000 /
-c     data diver(1),diver(2) / z00002500, z00000000 /
-c     data log10(1),log10(2) / z209a3f9a, zcffa84fb /
-c
-      if (i .lt. 1  .or.  i .gt. 5) goto 100
-c
-      d1mach = dmach(i)
-      return
-c
- 100  iwunit = i1mach(4)
-      write(iwunit, 99)
-  99  format(24hd1mach - i out of bounds)
-      stop
-      end
diff --git a/math/ieee/i1mach.f b/math/ieee/i1mach.f
deleted file mode 100644
index e7b9af2b1..000000000
--- a/math/ieee/i1mach.f
+++ /dev/null
@@ -1,382 +0,0 @@
-c
-c----------------------------------------------------------------------
-c  function:  i1mach
-c  this routine is from the port mathematical subroutine library
-c  it is described in the bell laboratories computing science
-c  technical report #47 by p.a. fox, a.d. hall and n.l. schryer
-c---------------------------------------------------------------------
-c
-      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 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.
-c
-      integer imach(16),output
-c
-      equivalence (imach(4),output)
-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 /
-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 /
-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 /
-c
-c     machine constants for 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) /   48 /
-c     data imach(12) / -974 /
-c     data imach(13) / 1070 /
-c     data imach(14) /   96 /
-c     data imach(15) / -927 /
-c     data imach(16) / 1070 /
-c
-c     machine constants for the cray 1
-c
-c     data imach( 1) /   100 /
-c     data imach( 2) /   101 /
-c     data imach( 3) /   102 /
-c     data imach( 4) /   101 /
-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) / -8192 /
-c     data imach(13) /  8190 /
-c     data imach(14) /    95 /
-c     data imach(15) / -8192 /
-c     data imach(16) /  8190 /
-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 /
-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 /
-c
-c     machine constants for the honeywell 600/6000 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) /    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) / -127 /
-c     data imach(13) /  127 /
-c     data imach(14) /   63 /
-c     data imach(15) / -127 /
-c     data imach(16) /  127 /
-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 /
-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) /    5 /
-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 /
-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) /    5 /
-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 /
-c
-c     machine constants for pdp-11 fortran supporting
-c     32-bit integer arithmetic.
-c
-      data imach( 1) /    5 /
-      data imach( 2) /    6 /
-      data imach( 3) /    6 /
-      data imach( 4) /    0 /
-      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 /
-c
-c     machine constants for pdp-11 fortran supporting
-c     16-bit integer arithmetic.
-c
-c     data imach( 1) /    5 /
-c     data imach( 2) /    6 /
-c     data imach( 3) /    5 /
-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 /
-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 /
-c
-c     machine constants for the vax-11 with
-c     fortran iv-plus compiler
-c
-c     data imach( 1) /    5 /
-c     data imach( 2) /    6 /
-c     data imach( 3) /    5 /
-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 /
-c
-      if (i .lt. 1  .or.  i .gt. 16) go to 10
-c
-      i1mach=imach(i)
-      return
-c
- 10   write(output,9000)
- 9000 format(39h1error    1 in i1mach - i out of bounds)
-c
-      stop
-c
-      end
diff --git a/math/ieee/r1mach.f b/math/ieee/r1mach.f
deleted file mode 100644
index db71aa682..000000000
--- a/math/ieee/r1mach.f
+++ /dev/null
@@ -1,177 +0,0 @@
-c
-c----------------------------------------------------------------------
-c  function:  r1mach
-c  this routine is from the port mathematical subroutine library
-c  it is described in the bell laboratories computing science
-c  technical report #47 by p.a. fox, a.d. hall and n.l. schryer
-c  a modification to the "i out of bounds" error message
-c  has been made by c. a. mcgonegal - april, 1978
-c----------------------------------------------------------------------
-c
-      real function r1mach(i)
-c
-c  single-precision machine constants
-c
-c  r1mach(1) = b**(emin-1), the smallest positive magnitude.
-c
-c  r1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude.
-c
-c  r1mach(3) = b**(-t), the smallest relative spacing.
-c
-c  r1mach(4) = b**(1-t), the largest relative spacing.
-c
-c  r1mach(5) = log10(b)
-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.
-c
-c  where possible, octal or hexadecimal constants have been used
-c  to specify the constants exactly which has in some cases
-c  required the use of equivalent integer arrays.
-c
-      integer small(2)
-      integer large(2)
-      integer right(2)
-      integer diver(2)
-      integer log10(2)
-c
-      real rmach(5)
-c
-      equivalence (rmach(1),small(1))
-      equivalence (rmach(2),large(1))
-      equivalence (rmach(3),right(1))
-      equivalence (rmach(4),diver(1))
-      equivalence (rmach(5),log10(1))
-c
-c     machine constants for the burroughs 1700 system.
-c
-c     data rmach(1) / z400800000 /
-c     data rmach(2) / z5ffffffff /
-c     data rmach(3) / z4e9800000 /
-c     data rmach(4) / z4ea800000 /
-c     data rmach(5) / z500e730e8 /
-c
-c     machine constants for the burroughs 5700/6700/7700 systems.
-c
-c     data rmach(1) / o1771000000000000 /
-c     data rmach(2) / o0777777777777777 /
-c     data rmach(3) / o1311000000000000 /
-c     data rmach(4) / o1301000000000000 /
-c     data rmach(5) / o1157163034761675 /
-c
-c     machine constants for the cdc 6000/7000 series.
-c
-c     data rmach(1) / 00014000000000000000b /
-c     data rmach(2) / 37767777777777777777b /
-c     data rmach(3) / 16404000000000000000b /
-c     data rmach(4) / 16414000000000000000b /
-c     data rmach(5) / 17164642023241175720b /
-c
-c     machine constants for the cray 1
-c
-c     data rmach(1) / 200004000000000000000b /
-c     data rmach(2) / 577767777777777777776b /
-c     data rmach(3) / 377224000000000000000b /
-c     data rmach(4) / 377234000000000000000b /
-c     data rmach(5) / 377774642023241175720b /
-c
-c     machine constants for the data general eclipse s/200
-c
-c     note - it may be appropriate to include the following card -
-c     static rmach(5)
-c
-c     data small/20k,0/,large/77777k,177777k/
-c     data right/35420k,0/,diver/36020k,0/
-c     data log10/40423k,42023k/
-c
-c     machine constants for the harris slash 6 and slash 7
-c
-c     data small(1),small(2) / '20000000, '00000201 /
-c     data large(1),large(2) / '37777777, '00000177 /
-c     data right(1),right(2) / '20000000, '00000352 /
-c     data diver(1),diver(2) / '20000000, '00000353 /
-c     data log10(1),log10(2) / '23210115, '00000377 /
-c
-c     machine constants for the honeywell 600/6000 series.
-c
-c     data rmach(1) / o402400000000 /
-c     data rmach(2) / o376777777777 /
-c     data rmach(3) / o714400000000 /
-c     data rmach(4) / o716400000000 /
-c     data rmach(5) / o776464202324 /
-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 rmach(1) / z00100000 /
-c     data rmach(2) / z7fffffff /
-c     data rmach(3) / z3b100000 /
-c     data rmach(4) / z3c100000 /
-c     data rmach(5) / z41134413 /
-c
-c     machine constants for the pdp-10 (ka or ki processor).
-c
-c     data rmach(1) / "000400000000 /
-c     data rmach(2) / "377777777777 /
-c     data rmach(3) / "146400000000 /
-c     data rmach(4) / "147400000000 /
-c     data rmach(5) / "177464202324 /
-c
-c     machine constants for pdp-11 fortran's supporting
-c     32-bit integers (expressed in integer and octal).
-c
-      data small(1) /    8388608 /
-      data large(1) / 2147483647 /
-      data right(1) /  880803840 /
-      data diver(1) /  889192448 /
-      data log10(1) / 1067065499 /
-c
-c     data rmach(1) / o00040000000 /
-c     data rmach(2) / o17777777777 /
-c     data rmach(3) / o06440000000 /
-c     data rmach(4) / o06500000000 /
-c     data rmach(5) / o07746420233 /
-c
-c     machine constants for pdp-11 fortran's supporting
-c     16-bit integers  (expressed in integer and octal).
-c
-c     data small(1),small(2) /   128,     0 /
-c     data large(1),large(2) / 32767,    -1 /
-c     data right(1),right(2) / 13440,     0 /
-c     data diver(1),diver(2) / 13568,     0 /
-c     data log10(1),log10(2) / 16282,  8347 /
-c
-c     data small(1),small(2) / o000200, o000000 /
-c     data large(1),large(2) / o077777, o177777 /
-c     data right(1),right(2) / o032200, o000000 /
-c     data diver(1),diver(2) / o032400, o000000 /
-c     data log10(1),log10(2) / o037632, o020233 /
-c
-c     machine constants for the univac 1100 series.
-c
-c     data rmach(1) / o000400000000 /
-c     data rmach(2) / o377777777777 /
-c     data rmach(3) / o146400000000 /
-c     data rmach(4) / o147400000000 /
-c     data rmach(5) / o177464202324 /
-c
-c     machine constants for the vax-11 with
-c     fortran iv-plus compiler
-c
-c     data rmach(1) / z00000080 /
-c     data rmach(2) / zffff7fff /
-c     data rmach(3) / z00003480 /
-c     data rmach(4) / z00003500 /
-c     data rmach(5) / z209b3f9a /
-c
-      if (i .lt. 1  .or.  i .gt. 5) goto 100
-c
-      r1mach = rmach(i)
-      return
-c
- 100  iwunit = i1mach(4)
-      write(iwunit, 99)
-  99  format(24hr1mach - i out of bounds)
-      stop
-      end
diff --git a/math/ieee/uni.c b/math/ieee/uni.c
deleted file mode 100644
index 8e76c099c..000000000
--- a/math/ieee/uni.c
+++ /dev/null
@@ -1,8 +0,0 @@
-/* Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
- */
-
-float uni_(dummy)
-float *dummy;
-{
-	return(rand()/ 2147483648.0);	/* return a real [0.0, 1.0) */
-}
diff --git a/math/math.hd b/math/math.hd
index b55579978..f6bc421a6 100644
--- a/math/math.hd
+++ b/math/math.hd
@@ -5,7 +5,6 @@ $bev		= "math$bevington/doc/"
 $curfit		= "math$curfit/doc/"
 $deboor		= "math$deboor/doc/"
 $gsurfit	= "math$gsurfit/doc/"
-$ieee		= "math$ieee/doc/"
 $iminterp	= "math$iminterp/doc/"
 $interp		= "math$interp/doc/"
 $llsq		= "math$llsq/doc/"
@@ -25,10 +24,6 @@ deboor		hlp = deboor$deboor.men,
 		sys = deboor$deboor.hlp,
 		pkg = deboor$deboor.hd
 
-ieee		hlp = ieee$ieee.men,
-		sys = ieee$ieee.hlp,
-		pkg = ieee$ieee.hd
-
 iminterp	hlp = iminterp$iminterp.men,
 		sys = iminterp$iminterp.hlp,
 		pkg = iminterp$iminterp.hd
diff --git a/math/math.men b/math/math.men
index 020a7a6ef..eebcc9d6d 100644
--- a/math/math.men
+++ b/math/math.men
@@ -2,7 +2,6 @@
 	 curfit - IRAF curve fitting package
 	 deboor - C.DeBoor's spline package
 	gsurfit - IRAF general surface fitting package
-	   ieee - IEEE signal processing routines
        iminterp	- IRAF image interpolation package
 	 interp - Old 1D image interpolation routines
 	 lapack - Numerical linear algebra package
diff --git a/math/mkpkg b/math/mkpkg
index 89ffb93d7..403aa1fac 100644
--- a/math/mkpkg
+++ b/math/mkpkg
@@ -16,7 +16,6 @@ mathgen:
 	$call curfit
 	$call deboor
 	$call gsurfit
-	# $call ieee
 	$call iminterp
 	$call interp
 	$call llsq
@@ -52,12 +51,6 @@ gsurfit:
 	$update   libgsurfit.a
 	$checkin  libgsurfit.a lib$
 	;
-ieee:
-	$echo "-------------- LIBIEEE ------------------"
-	$checkout libieee.a lib$
-	$update   libieee.a
-	$checkin  libieee.a lib$
-	;
 iminterp:
 	$echo "-------------- LIBIMINTERP --------------"
 	$checkout libiminterp.a lib$
@@ -125,10 +118,6 @@ libgsurfit.a:				# Generalized 2d surface fitting pkg
 	@gsurfit
 	;
 
-libieee.a:				# IEEE signal processing package
-	@ieee
-	;
-
 libiminterp.a:				# Image interpolation package
 	@iminterp
 	;