Repatch (from linux) of OSX IRAF

51 files changed, 5664 insertions, 0 deletions

diff --git a/math/ieee/README b/math/ieee/README
new file mode 100644
index 00000000..80b30bfd
--- /dev/null
+++ b/math/ieee/README
@@ -0,0 +1,8 @@
+     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
new file mode 100644
index 00000000..bf86b45f
--- /dev/null
+++ b/math/ieee/chap1/README
@@ -0,0 +1,14 @@
+     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
new file mode 100644
index 00000000..54eb4e41
--- /dev/null
+++ b/math/ieee/chap1/const.f
@@ -0,0 +1,114 @@
+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
new file mode 100644
index 00000000..9a0ad713
--- /dev/null
+++ b/math/ieee/chap1/fast.f
@@ -0,0 +1,73 @@
+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
new file mode 100644
index 00000000..dd62ea01
--- /dev/null
+++ b/math/ieee/chap1/ffa.f
@@ -0,0 +1,84 @@
+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
new file mode 100644
index 00000000..2858fdb0
--- /dev/null
+++ b/math/ieee/chap1/ffs.f
@@ -0,0 +1,80 @@
+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
new file mode 100644
index 00000000..37a3a651
--- /dev/null
+++ b/math/ieee/chap1/fft842.f
@@ -0,0 +1,116 @@
+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
new file mode 100644
index 00000000..f703c98b
--- /dev/null
+++ b/math/ieee/chap1/fftaoh.f
@@ -0,0 +1,82 @@
+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
new file mode 100644
index 00000000..30c84b79
--- /dev/null
+++ b/math/ieee/chap1/fftasm.f
@@ -0,0 +1,67 @@
+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
new file mode 100644
index 00000000..6e3ecc9d
--- /dev/null
+++ b/math/ieee/chap1/fftohm.f
@@ -0,0 +1,101 @@
+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
new file mode 100644
index 00000000..afc18c17
--- /dev/null
+++ b/math/ieee/chap1/fftsoh.f
@@ -0,0 +1,81 @@
+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
new file mode 100644
index 00000000..b1e795bb
--- /dev/null
+++ b/math/ieee/chap1/fftsym.f
@@ -0,0 +1,84 @@
+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
new file mode 100644
index 00000000..4982246f
--- /dev/null
+++ b/math/ieee/chap1/ford1.f
@@ -0,0 +1,24 @@
+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
new file mode 100644
index 00000000..8ee26d69
--- /dev/null
+++ b/math/ieee/chap1/ford2.f
@@ -0,0 +1,46 @@
+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
new file mode 100644
index 00000000..d412c0e2
--- /dev/null
+++ b/math/ieee/chap1/fourea.f
@@ -0,0 +1,98 @@
+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
new file mode 100644
index 00000000..24a103ff
--- /dev/null
+++ b/math/ieee/chap1/fr2tr.f
@@ -0,0 +1,15 @@
+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
new file mode 100644
index 00000000..5ce978f0
--- /dev/null
+++ b/math/ieee/chap1/fr4syn.f
@@ -0,0 +1,109 @@
+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
new file mode 100644
index 00000000..16ba3fe6
--- /dev/null
+++ b/math/ieee/chap1/fr4tr.f
@@ -0,0 +1,118 @@
+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
new file mode 100644
index 00000000..75ce56cb
--- /dev/null
+++ b/math/ieee/chap1/fsst.f
@@ -0,0 +1,71 @@
+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
new file mode 100644
index 00000000..8cb0a908
--- /dev/null
+++ b/math/ieee/chap1/iftaoh.f
@@ -0,0 +1,87 @@
+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
new file mode 100644
index 00000000..6ec12c7d
--- /dev/null
+++ b/math/ieee/chap1/iftasm.f
@@ -0,0 +1,77 @@
+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
new file mode 100644
index 00000000..df084e87
--- /dev/null
+++ b/math/ieee/chap1/iftohm.f
@@ -0,0 +1,83 @@
+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
new file mode 100644
index 00000000..6098e6e4
--- /dev/null
+++ b/math/ieee/chap1/iftsoh.f
@@ -0,0 +1,94 @@
+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
new file mode 100644
index 00000000..e45c9a37
--- /dev/null
+++ b/math/ieee/chap1/iftsym.f
@@ -0,0 +1,90 @@
+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
new file mode 100644
index 00000000..925cd657
--- /dev/null
+++ b/math/ieee/chap1/inishl.f
@@ -0,0 +1,179 @@
+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
new file mode 100644
index 00000000..21e4494e
--- /dev/null
+++ b/math/ieee/chap1/ord1.f
@@ -0,0 +1,24 @@
+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
new file mode 100644
index 00000000..10f03ec7
--- /dev/null
+++ b/math/ieee/chap1/ord2.f
@@ -0,0 +1,46 @@
+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
new file mode 100644
index 00000000..510763fe
--- /dev/null
+++ b/math/ieee/chap1/r2tr.f
@@ -0,0 +1,16 @@
+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
new file mode 100644
index 00000000..425c10f8
--- /dev/null
+++ b/math/ieee/chap1/r2tx.f
@@ -0,0 +1,18 @@
+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
new file mode 100644
index 00000000..77c97bba
--- /dev/null
+++ b/math/ieee/chap1/r4syn.f
@@ -0,0 +1,20 @@
+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
new file mode 100644
index 00000000..5fc46eab
--- /dev/null
+++ b/math/ieee/chap1/r4tr.f
@@ -0,0 +1,18 @@
+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
new file mode 100644
index 00000000..4e5649e8
--- /dev/null
+++ b/math/ieee/chap1/r4tx.f
@@ -0,0 +1,29 @@
+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
new file mode 100644
index 00000000..054413d7
--- /dev/null
+++ b/math/ieee/chap1/r8syn.f
@@ -0,0 +1,186 @@
+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
new file mode 100644
index 00000000..d49ef58c
--- /dev/null
+++ b/math/ieee/chap1/r8tr.f
@@ -0,0 +1,201 @@
+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
new file mode 100644
index 00000000..9cc4f591
--- /dev/null
+++ b/math/ieee/chap1/r8tx.f
@@ -0,0 +1,107 @@
+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
new file mode 100644
index 00000000..dfebf0cc
--- /dev/null
+++ b/math/ieee/chap1/rad4sb.f
@@ -0,0 +1,38 @@
+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
new file mode 100644
index 00000000..cd703305
--- /dev/null
+++ b/math/ieee/chap1/radix4.f
@@ -0,0 +1,488 @@
+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
new file mode 100644
index 00000000..a830ab2c
--- /dev/null
+++ b/math/ieee/chap1/test/test12.f
@@ -0,0 +1,90 @@
+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
new file mode 100644
index 00000000..d7157069
--- /dev/null
+++ b/math/ieee/chap1/test/test13.f
@@ -0,0 +1,260 @@
+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
new file mode 100644
index 00000000..15a34788
--- /dev/null
+++ b/math/ieee/chap1/test/test17.f
@@ -0,0 +1,147 @@
+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
new file mode 100644
index 00000000..cf550e01
--- /dev/null
+++ b/math/ieee/chap1/test/test18.f
@@ -0,0 +1,71 @@
+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
new file mode 100644
index 00000000..774c2cc9
--- /dev/null
+++ b/math/ieee/chap1/time/time12.f
@@ -0,0 +1,53 @@
+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
new file mode 100644
index 00000000..adcce4c2
--- /dev/null
+++ b/math/ieee/chap1/time/time17.f
@@ -0,0 +1,53 @@
+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
new file mode 100644
index 00000000..87680554
--- /dev/null
+++ b/math/ieee/chap1/time/time18.f
@@ -0,0 +1,48 @@
+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
new file mode 100644
index 00000000..9c83d0ea
--- /dev/null
+++ b/math/ieee/chap1/weave1.f
@@ -0,0 +1,371 @@
+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
new file mode 100644
index 00000000..5c04ff91
--- /dev/null
+++ b/math/ieee/chap1/weave2.f
@@ -0,0 +1,412 @@
+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
new file mode 100644
index 00000000..9e819941
--- /dev/null
+++ b/math/ieee/chap1/wfta.f
@@ -0,0 +1,150 @@
+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
new file mode 100644
index 00000000..699468ab
--- /dev/null
+++ b/math/ieee/d1mach.f
@@ -0,0 +1,256 @@
+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
new file mode 100644
index 00000000..e7b9af2b
--- /dev/null
+++ b/math/ieee/i1mach.f
@@ -0,0 +1,382 @@
+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
new file mode 100644
index 00000000..db71aa68
--- /dev/null
+++ b/math/ieee/r1mach.f
@@ -0,0 +1,177 @@
+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
new file mode 100644
index 00000000..8e76c099
--- /dev/null
+++ b/math/ieee/uni.c
@@ -0,0 +1,8 @@
+/* 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) */
+}