Initial population.

1 files changed, 1174 insertions, 0 deletions

diff --git a/src/algebra/ndftip.as.pamphlet b/src/algebra/ndftip.as.pamphlet
new file mode 100644
index 00000000..3c9c0ea5
--- /dev/null
+++ b/src/algebra/ndftip.as.pamphlet
@@ -0,0 +1,1174 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra ndftip.as}
+\author{Michael Richardson}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{NagDiscreteFourierTransformInterfacePackage}
+<<NagDiscreteFourierTransformInterfacePackage>>=
++++ Author: M.G. Richardson
++++ Date Created: 1995 Dec. 08
++++ Date Last Updated:
++++ Basic Functions:
++++ Related Constructors:
++++ Also See:
++++ AMS Classifications:
++++ Keywords:
++++ References:
++++ Description:
++++ This package provides Axiom-like interfaces to the NAG
++++ Finite Fourier Transform routines in the NAGlink.
+
+NagDiscreteFourierTransformInterfacePackage: with {
+
+  nagDFT : VDF -> VCDF ;                                --  test  1
+
+++ nagDFT(seq) calculates the discrete Fourier transform of a sequence
+++ of real data values 
+#if saturn
+++ $x_{1} \ldots x_{n}$
+#else
+++ \spad{x[1] .. x[n]}
+#endif
+++ supplied in the vector seq.
+++ Note that the definition used for the discrete Fourier transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{-i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(x[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06EAF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06eaf.
+
+  nagDFT : VCDF -> VCDF ;                               --  test  3
+
+++ nagDFT(seq) calculates the discrete Fourier transform of a sequence
+++ of complex data values 
+#if saturn
+++ $z_{1} \ldots z_{n}$
+#else
+++ \spad{z[1] .. z[n]}
+#endif
+++ supplied in the vector seq.
+++ Note that the definition used for the discrete Fourier transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{-i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(z[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06ECF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06ecf.
+
+  nagDFT : PHSDF -> VDF ;                               --  test  7
+
+++ nagDFT(seq) calculates the discrete Fourier transform of a Hermitian
+++ sequence of complex data values,
+#if saturn
+++ $z_{1} \ldots z_{n}$
+#else
+++ \spad{z[1] .. z[n]}
+#endif
+++ supplied in the PackedHermitianSequence seq.
+++ Note that the definition used for the discrete Fourier transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{-i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(z[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06EBF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06ebf.
+
+  nagDFT : LVDF -> LVCDF ;                               --  test 10, 19
+
+++ nagDFT(seqs) calculates the discrete Fourier transform of each of a
+++ list of sequences of real data values 
+#if saturn
+++ $x_{1} \ldots x_{n}$
+#else
+++ \spad{x[1] .. x[n]}
+#endif
+++ supplied in the list of vectors, seqs.
+++ Note that the definition used for the discrete Fourier transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{-i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(x[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06FPF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06fpf.
+
+  nagDFT : LVCDF -> LVCDF ;                             --  test 16
+
+++ nagDFT(seqs) calculates the discrete Fourier transform of each of a
+++ list of sequences of complex data values 
+#if saturn
+++ $z_{1} \ldots z_{n}$
+#else
+++ \spad{z[1] .. z[n]}
+#endif
+++ supplied in the list of vectors, seqs.
+++ Note that the definition used for the discrete Fourier transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{-i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(z[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06FRF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06frf.
+
+  nagDFT : LPHSDF -> LVDF ;                              --  test 12, 21
+
+++ nagDFT(seq) calculates the discrete Fourier transform of a each of a
+++ list of Hermitian sequences of complex data values,
+#if saturn
+++ $z_{1} \ldots z_{n}$
+#else
+++ \spad{z[1] .. z[n]}
+#endif
+++ supplied in the List PackedHermitianSequence, seq.
+++ Note that the definition used for the discrete Fourier transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{-i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(z[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06FQF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06fqf.
+
+  nagInverseDFT : VDF -> VCDF ;                         --  test  8
+
+++ nagInverseDFT(seq) calculates the inverse discrete Fourier
+++ transform of a sequence of real data values 
+#if saturn
+++ $x_{1} \ldots x_{n}$
+#else
+++ \spad{x[1] .. x[n]}
+#endif
+++ supplied in the vector seq.
+++ Note that the definition used for the inverse discrete Fourier
+++ transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(x[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06EAF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06eaf.
+
+  nagInverseDFT : VCDF -> VCDF ;                         --  test  2,  4
+
+++ nagInverseDFT(seq) calculates the inverse discrete Fourier
+++ transform of a sequence of complex data values 
+#if saturn
+++ $z_{1} \ldots z_{n}$
+#else
+++ \spad{z[1] .. z[n]}
+#endif
+++ supplied in the vector seq.
+++ Note that the definition used for the inverse discrete Fourier
+++ transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(z[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06ECF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06ecf.
+
+  nagInverseDFT : PHSDF -> VDF ;                        --  test  6
+
+++ nagInverseDFT(seq) calculates the inverse discrete Fourier transform
+++ of a Hermitian sequence of complex data values 
+#if saturn
+++ $z_{1} \ldots z_{n}$
+#else
+++ \spad{z[1] .. z[n]}
+#endif
+++ supplied in the PackedHermitianSequence seq.
+++ Note that the definition used for the inverse discrete Fourier
+++ transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(z[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06EBF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06ebf.
+
+  nagInverseDFT : LVDF -> LVCDF ;                       --  test 13
+
+++ nagInverseDFT(seqs) calculates the inverse discrete Fourier
+++ transform of each of a list of sequences of real data values 
+#if saturn
+++ $x_{1} \ldots x_{n}$
+#else
+++ \spad{x[1] .. x[n]}
+#endif
+++ supplied in the list of vectors, seqs.
+++ Note that the definition used for the inverse discrete Fourier
+++ transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(x[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06FPF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06fpf.
+
+  nagInverseDFT : LVCDF -> LVCDF ;                       --  test 11, 17
+
+++ nagInverseDFT(seqs) calculates the inverse discrete Fourier
+++ transform of each of a list of sequences of complex data values 
+#if saturn
+++ $z_{1} \ldots z_{n}$
+#else
+++ \spad{z[1] .. z[n]}
+#endif
+++ supplied in the list of vectors, seqs.
+++ Note that the definition used for the inverse discrete Fourier
+++ transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(z[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06FRF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06frf.
+
+  nagInverseDFT : LPHSDF -> LVDF ;                      --  test 15
+
+++ nagInverseDFT(seqs) calculates the inverse discrete Fourier transform
+++ of each of a list of Hermitian sequences of complex data values 
+#if saturn
+++ $z_{1} \ldots z_{n}$
+#else
+++ \spad{z[1] .. z[n]}
+#endif
+++ supplied in the List PackedHermitianSequence, seqs.
+++ Note that the definition used for the inverse discrete Fourier
+++ transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(z[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06FQF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06fqf.
+
+  nagHermitianDFT : VDF -> PHSDF ;                      --  test  5
+
+++ nagHermitianDFT(seq) calculates the discrete Fourier transform, in
+++ packed Hermitian form, of a sequence of real data values 
+#if saturn
+++ $x_{1} \ldots x_{n}$
+#else
+++ \spad{x[1] .. x[n]}
+#endif
+++ supplied in the vector seq.
+++ Note that the definition used for the discrete Fourier transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{-i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(x[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06EAF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06eaf.
+
+  nagHermitianDFT : LVDF -> LPHSDF ;                     --  test 14, 20
+
+++ nagHermitianDFT(seqs) calculates the discrete Fourier transform, in
+++ packed Hermitian form, of each of a list of sequences of real data
+++ values 
+#if saturn
+++ $x_{1} \ldots x_{n}$
+#else
+++ \spad{x[1] .. x[n]}
+#endif
+++ supplied in the list of vectors, seqs.
+++ Note that the definition used for the discrete Fourier transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{-i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(x[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06FPF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06fpf.
+
+  nagHermitianInverseDFT : VDF -> PHSDF ;               --  test  9
+
+++ nagHermitianInverseDFT(seq) calculates the inverse discrete Fourier
+++ transform, in packed Hermitian form, of a sequence of real data
+++ values 
+#if saturn
+++ $x_{1} \ldots x_{n}$
+#else
+++ \spad{x[1] .. x[n]}
+#endif
+++ supplied in the vector seq.
+++ Note that the definition used for the inverse discrete Fourier
+++ transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(x[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06EAF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06eaf.
+
+  nagHermitianInverseDFT : LVDF -> LPHSDF ;             --  test 18
+
+++ nagHermitianInverseDFT(seqs) calculates the inverse discrete Fourier
+++ transform, in packed Hermitian form, of each of a list of sequences
+++ of real data values 
+#if saturn
+++ $x_{1} \ldots x_{n}$
+#else
+++ \spad{x[1] .. x[n]}
+#endif
+++ supplied in the list of vectors, seqs.
+++ Note that the definition used for the inverse discrete Fourier
+++ transform is
+#if saturn
+++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{i \frac{2 \pi j k}{n}
+++ \qquad k = 0 \ldots  n - 1 \]
+#else
+++ \spad{1/sqrt(n)*sum(x[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
+++ \spad{k=0..(n-1)}.
+#endif
+++ The numerical calculation is performed by the NAG routine C06FPF.
+++
+++ For more detailed information, please consult the NAG
+++ manual via the Browser page for the operation c06fpf.
+
+} == add {
+
+  import from AnyFunctions1 MDF ;
+  import from CDF;
+  import from ErrorFunctions ;
+  import from LLDF ;
+  import from MCDF ;
+  import from MDF ;
+  import from NagResultChecks ;
+  import from NagSeriesSummationPackage ;
+  import from PHSDF;
+  import from STRG ;
+  import from List STRG ;
+  import from Symbol ;
+  import from VDF ;
+
+  local (..)(a:INT,b:INT):Generator INT == {
+    generate {
+      t := a ;
+      while (t <= b) repeat {
+        yield t ;
+        t := t + 1 ;
+        }
+      }
+    }
+
+  local ipIfail : INT := -1 ;
+
+-- First, the functions corresponding to single NAGlink calls of C06E
+-- routines (single vector transforms):
+
+-- c06eaf:
+
+  nagHermitianDFT(seq : VDF) : PHSDF ; == {
+    local lseq : INT ;
+
+    lseq := ((# seq)@NNI) pretend INT ; -- @ to eliminate SI possibility
+    row(checkMxDF(c06eaf(lseq,matrix [members seq],ipIfail),
+		   "x",
+		    "C06EAF"),
+			      1)
+				 pretend PHSDF
+  }
+
+-- c06ebf:
+
+  nagDFT(seq : PHSDF) : VDF == {
+    local lseq : INT ;
+
+    lseq := ((# seq)@NNI) pretend INT ; -- @ to eliminate SI possibility
+    row(checkMxDF(c06ebf(lseq,matrix [members seq],ipIfail),
+		   "x",
+		    "C06EBF"),
+			      1)
+  }
+
+-- c06ecf:
+
+  nagDFT(seq : VCDF) : VCDF == {
+    local nseq : NNI ;
+    local lseq : INT ;
+    local rvec, ivec : VDF ;
+    local cvec : VCDF ;
+    local c06ecfResult : RSLT ;
+
+    nseq := # seq ;
+    lseq := nseq pretend INT ;
+    rvec := new(nseq,0) ;
+    ivec := new(nseq,0) ;
+    for i in 1..lseq repeat {
+      rvec(i) := real seq(i) ;
+      ivec(i) := imag seq(i) ;
+    }
+    c06ecfResult := c06ecf(lseq,
+			    matrix [members rvec],
+			     matrix [members ivec],
+			      ipIfail) ;
+    rvec := row(checkMxDF(c06ecfResult,"x","C06ECF"),1) ;
+    ivec := row((retract(c06ecfResult."y") @ MDF),1) ;
+    cvec := new(nseq,0) ;
+    for i in 1..lseq repeat cvec(i) := complex(rvec(i),ivec(i)) ;
+    cvec
+  }
+
+-- inverse transforms, in terms of these and functions from PHS:
+
+  nagInverseDFT(seq : PHSDF) : VDF == nagDFT conjHerm seq ;
+
+  nagHermitianInverseDFT(seq : VDF) : PHSDF
+			== conjHerm nagHermitianDFT seq ;
+
+  nagInverseDFT(seq : VCDF) : VCDF == {
+    local nseq : NNI ;
+    local lseq : INT ;
+    local rvec, ivec : VDF ;
+    local cvec : VCDF ;
+    local c06ecfResult : RSLT ;
+
+    nseq := # seq ;
+    lseq := nseq pretend INT ;
+    rvec := new(nseq,0) ;
+    ivec := new(nseq,0) ;
+    for i in 1..lseq repeat {
+      rvec(i) := real seq(i) ;
+      ivec(i) := - imag seq(i) ;
+    }
+    c06ecfResult := c06ecf(lseq,
+			    matrix [members rvec],
+			     matrix [members ivec],
+			      ipIfail) ;
+    rvec := row(checkMxDF(c06ecfResult,"x","C06ECF"),1) ;
+    ivec := row((retract(c06ecfResult."y") @ MDF),1) ;
+    cvec := new(nseq,0) ;
+    for i in 1..lseq repeat cvec(i) := complex(rvec(i), - ivec(i)) ;
+    cvec
+  }
+
+-- "Full form" equivalents of c06eaf and inverse:
+
+  nagDFT(seq : VDF) : VCDF == expand nagHermitianDFT seq ;
+
+  nagInverseDFT(seq : VDF) : VCDF == expand nagHermitianInverseDFT seq ;
+
+
+-- Next, the functions corresponding to single NAGlink calls of C06F
+-- routines (multiple vector transforms):
+
+-- basic routines:
+
+-- c06fpf
+
+  nagHermitianDFT(seqs : LVDF) : LPHSDF ; == {
+
+    local nr, nc : NNI ;
+    local inr, inc : INT ;
+    local seqMat, trig, result : MDF ;
+    local nextSeq : PHSDF ;
+    local hermDFTs : LPHSDF ;
+
+    nr := # seqs ;
+    inr := nr pretend INT ;
+    nc := # (seqs.1) ;
+    inc := nc pretend INT ;
+    seqMat := new(nr,nc,0) ;
+    for j in 1 .. inc repeat seqMat(1,j) := (seqs.1).j ;
+    for i in 2 .. inr repeat
+      if (# seqs.i) ~= nc 
+      then error ["The data sequences in nagHermitianDFT must all",
+		  " have the same length.  ",
+		  "The length of sequence 1 is ",
+		  string(inc),
+		  "that of sequence ",
+		  string(i pretend INT),
+		  " is ",
+		  string((# seqs.i)@NNI pretend INT), -- @ avoids SI
+		  "."]
+      else for j in 1 .. inc repeat seqMat(i,j) := (seqs.i).j ;
+    trig := new(1@NNI,2*nc,0) ;
+    result :=
+      checkMxDF(c06fpf(inr,inc,"i",seqMat,trig,ipIfail),"x","C06FPF") ;
+    hermDFTs := [] ;
+    for i in inr .. 1  by -1 repeat {
+      nextSeq := new(nc,0) ;
+      for j in 1 .. inc repeat nextSeq(j) := result(1,(j-1)*inr + i) ;
+      hermDFTs := cons(nextSeq,hermDFTs) ;
+    }
+    hermDFTs
+  }
+
+-- c06fqf
+
+  nagDFT(seqs : LPHSDF) : LVDF == {
+    
+    local nr, nc : NNI ;
+    local inr, inc : INT ;
+    local seqMat, trig, result : MDF ;
+    local nextSeq : VDF ;
+    local dfts : LVDF ;
+  
+    nr := # seqs ;
+    inr := nr pretend INT ;
+    nc := # (seqs.1) ;
+    inc := nc pretend INT ;
+    seqMat := new(nr,nc,0) ;
+    for j in 1 .. inc repeat seqMat(1,j) := (seqs.1).j ;
+    for i in 2 .. inr repeat
+      if (# seqs.i) ~= nc 
+      then error ["The data sequences in nagDFT must all",
+                  " have the same length.  ",
+                  "The length of sequence 1 is ",
+                  string(inc),
+                  "that of sequence ",
+                  string(i pretend INT),
+                  " is ",
+                  string((# seqs.i)@NNI pretend INT), -- @ avoids SI
+                  "."]
+      else for j in 1 .. inc repeat seqMat(i,j) := (seqs.i).j ;
+    trig := new(1@NNI,2*nc,0) ;
+    result :=
+      checkMxDF(c06fqf(inr,inc,"i",seqMat,trig,ipIfail),"x","C06FQF") ;
+    dfts := [] ;
+    for i in inr .. 1 by -1 repeat {
+      nextSeq := new(nc,0) ;
+      for j in 1 .. inc repeat nextSeq(j) := result(1,(j-1)*inr + i) ;
+      dfts := cons(nextSeq,dfts) ;
+    }
+    dfts
+  }
+
+-- c06frf
+
+  nagDFT(seqs : LVCDF) : LVCDF == {
+
+    local nr, nc : NNI ;
+    local inr, inc : INT ;
+    local trig, rMat, iMat : MDF ;
+    local result : RSLT ;
+    local nextSeq : VCDF ;
+    local dfts : LVCDF ;
+
+    nr := # seqs ;
+    inr := nr pretend INT ;
+    nc := # (seqs.1) ;
+    inc := nc pretend INT ;
+    rMat := new(nr,nc,0) ;
+    iMat := new(nr,nc,0) ;
+    for j in 1 .. inc repeat {
+      rMat(1,j) := real((seqs.1).j) ;
+      iMat(1,j) := imag((seqs.1).j) ;
+    }
+    for i in 2 .. inr repeat {
+      if (# seqs.i) ~= nc
+      then error ["The data sequences in nagDFT must all",
+		  " have the same length.  ",
+		  "The length of sequence 1 is ",
+		  string(inc),
+		  "that of sequence ",
+		  string(i pretend INT),
+		  " is ",
+		  string((# seqs.i)@NNI pretend INT), -- @ avoids SI
+		  "."]
+      else for j in 1 .. inc repeat {
+	rMat(i,j) := real((seqs.i).j) ;
+	iMat(i,j) := imag((seqs.i).j) ;
+      }
+    }
+    trig := new(1@NNI,2*nc,0) ;
+    result := c06frf(inr,inc,"i",rMat,iMat,trig,ipIfail) ;
+    rMat := checkMxDF(result, "x", "C06FRF") ;
+    iMat := retract(result."y") @ MDF ;
+    dfts := [] ;
+    for i in inr .. 1 by -1 repeat {
+      nextSeq := new(nc,0) ;
+      for j in 1 .. inc repeat
+	nextSeq(j) := complex(rMat(1,(j-1)*inr+i),iMat(1,(j-1)*inr+i)) ;
+      dfts := cons(nextSeq,dfts) ;
+    }
+    dfts
+  }
+
+-- inverse transforms, in terms of these and functions from PHS:
+
+  nagInverseDFT(seqs : LVCDF) : LVCDF == {
+
+    local nr, nc : NNI ;
+    local inr, inc : INT ;
+    local conjSeq : VCDF ;
+    local temp, invdfts : LVCDF ;
+
+    nr := # seqs ;
+    inr := nr pretend INT ;
+    temp := [] ;
+    for i in inr .. 1 by -1 repeat {
+      nc := #(seqs.i) ;
+      inc := nc pretend INT ;
+      conjSeq := new(nc,0) ;
+      for j in 1 .. inc repeat
+        conjSeq(j) := conjugate((seqs.i).j) ;
+      temp := cons(conjSeq,temp) ;
+    }
+    temp := nagDFT temp ;
+    invdfts := [] ;
+    for i in inr .. 1 by -1 repeat {
+      conjSeq := new(nc,0) ;
+      for j in 1 .. inc repeat -- know inc is constant after nagDFT call
+	conjSeq(j) := conjugate((temp.i).j) ;
+      invdfts := cons(conjSeq,invdfts) ;
+    }
+    invdfts
+  }
+
+  nagInverseDFT(seqs : LPHSDF) : LVDF == {
+    local nr : NNI ;
+    local inr : INT ;
+    local conjSeqs : LPHSDF ;
+  
+    nr := # seqs ;
+    inr := nr pretend INT ;
+    conjSeqs := [] ;
+    for i in inr .. 1 by -1 repeat
+      conjSeqs := cons(conjHerm(seqs.i),conjSeqs) ;
+    nagDFT conjSeqs ;
+  }
+
+  nagHermitianInverseDFT(seqs : LVDF) : LPHSDF == {
+    local nr : NNI ;
+    local inr : INT ;
+    local conjSeqs, invSeqs : LPHSDF ;
+  
+    nr := # seqs ;
+    inr := nr pretend INT ;
+    conjSeqs := nagHermitianDFT seqs ;
+    invSeqs := [] ;
+    for i in inr .. 1 by -1 repeat
+      invSeqs := cons(conjHerm(conjSeqs.i),invSeqs) ;
+    invSeqs
+  }
+
+-- "Full form" equivalents of c06fpf and inverse:
+
+  nagDFT(seqs : LVDF) : LVCDF == {
+  
+    local nr : NNI ;
+    local inr : INT ;
+    local hermdfts : LPHSDF ;
+    local dfts : LVCDF ;
+
+    nr := # seqs ;
+    inr := nr pretend INT ;
+    hermdfts := nagHermitianDFT seqs ;
+    dfts := [] ;
+    for i in inr .. 1 by -1 repeat
+      dfts := cons(expand(hermdfts.i),dfts) ;
+    dfts
+  }
+
+  nagInverseDFT(seqs : LVDF) : LVCDF == {
+    local nr : NNI ;
+    local inr : INT ;
+    local hermdfts : LPHSDF ;
+    local invdfts : LVCDF ;
+
+    nr := # seqs ;
+    inr := nr pretend INT ;
+    hermdfts := nagHermitianDFT seqs ;
+    invdfts := [] ;
+    for i in inr .. 1 by -1 repeat
+      invdfts := cons(expand conjHerm(hermdfts.i),invdfts) ;
+    invdfts
+  }
+    
+}
+
+#if NeverAssertThis
+
+-- Note that the conversions of results from DoubleFloat to Float
+-- will become unnecessary if outputGeneral is extended to apply to
+-- DoubleFloat quantities.  Those results not converted will, of
+-- course, then be displayed to 6 s.f.
+
+)lib nrc
+)lib herm
+)lib ndftip
+
+outputGeneral 6
+
+seqA := [0.34907,0.54890,0.74776,0.94459,1.1385,1.3285,1.5137];
+
+seqB := [0.34907 - 0.37168*%i,  _
+         0.54890 - 0.35669*%i,  _
+         0.74776 - 0.31175*%i,  _
+         0.94459 - 0.23702*%i,  _
+         1.13850 - 0.13274*%i,  _
+         1.32850 + 0.00074*%i,  _
+         1.51370 + 0.16298*%i];
+
+hseqC : PackedHermitianSequence DoubleFloat
+hseqC := packHS [0.34907,        _
+                 0.54890 + %i*1.51370,  _
+                 0.74776 + %i*1.32850,  _
+                 0.94459 + %i*1.13850,  _
+                 0.94459 - %i*1.13850,  _
+                 0.74776 - %i*1.32850,  _
+                 0.54890 - %i*1.51370];
+
+seqsD : List Vector DoubleFloat;
+seqsD := [vector [0.3854, 0.6772, 0.1138, 0.6751, 0.6362, 0.1424], _
+          vector [0.5417, 0.2983, 0.1181, 0.7255, 0.8638, 0.8723], _
+          vector [0.9172, 0.0644, 0.6037, 0.6430, 0.0428, 0.4815]];
+
+seqsE : List PackedHermitianSequence DoubleFloat;
+seqsE := [pHS [0.3854, 0.6772, 0.1138, 0.6751, 0.6362, 0.1424], _
+          pHS [0.5417, 0.2983, 0.1181, 0.7255, 0.8638, 0.8723], _
+          pHS [0.9172, 0.0644, 0.6037, 0.6430, 0.0428, 0.4815]];
+
+seqsF : List Vector Complex DoubleFloat
+seqsF := [vector [0.3854 + 0.5417*%i, 0.6772 + 0.2983*%i,   _
+                  0.1138 + 0.1181*%i, 0.6751 + 0.7255*%i,   _
+                  0.6362 + 0.8638*%i, 0.1424 + 0.8723*%i],  _
+          vector [0.9172 + 0.9089*%i, 0.0644 + 0.3118*%i,   _
+                  0.6037 + 0.3465*%i, 0.6430 + 0.6198*%i,   _
+                  0.0428 + 0.2668*%i, 0.4815 + 0.1614*%i],  _
+          vector [0.1156 + 0.6214*%i, 0.0685 + 0.8681*%i,   _
+                  0.2060 + 0.7060*%i, 0.8630 + 0.8652*%i,   _
+                  0.6967 + 0.9190*%i, 0.2792 + 0.3355*%i]];
+
+-- test  1
+
+dftA := nagDFT seqA;
+dftA :: Vector Complex Float :: Matrix Complex Float
+                             -- Matrix to force display as a column,
+                             -- Float to allow outputGeneral to work.
+
+--       +         2.48361         +
+--       |                         |
+--       |- 0.265985 + 0.530898 %i |
+--       |                         |
+--       |- 0.257682 + 0.202979 %i |
+--       |                         |
+--       |- 0.256363 + 0.0580623 %i|
+--       |                         |
+--       |- 0.256363 - 0.0580623 %i|
+--       |                         |
+--       |- 0.257682 - 0.202979 %i |
+--       |                         |
+--       +- 0.265985 - 0.530898 %i +
+
+-- test  2
+
+nagInverseDFT dftA :: Vector Float
+ 
+--       [0.34907,0.5489,0.74776,0.94459,1.1385,1.3285,1.5137]
+
+-- test  3
+
+dftB := nagDFT seqB;
+dftB :: Vector Complex Float :: Matrix Complex Float
+
+--       +  2.48361 - 0.471004 %i  +
+--       |                         |
+--       | - 0.5518 + 0.496841 %i  |
+--       |                         |
+--       |- 0.367113 + 0.0975621 %i|
+--       |                         |
+--       |- 0.287669 - 0.0586476 %i|
+--       |                         |
+--       |- 0.225057 - 0.174772 %i |
+--       |                         |
+--       |- 0.148251 - 0.308396 %i |
+--       |                         |
+--       + 0.0198297 - 0.564956 %i +
+ 
+-- test  4
+ 
+(nagInverseDFT dftB) :: Vector Complex Float :: Matrix Complex Float
+ 
+--       +0.34907 - 0.37168 %i+
+--       |                    |
+--       |0.5489 - 0.35669 %i |
+--       |                    |
+--       |0.74776 - 0.31175 %i|
+--       |                    |
+--       |0.94459 - 0.23702 %i|
+--       |                    |
+--       |1.1385 - 0.13274 %i |
+--       |                    |
+--       |1.3285 + 0.00074 %i |
+--       |                    |
+--       +1.5137 + 0.16298 %i +
+
+-- test  5
+
+hdftA := nagHermitianDFT seqA;
+(expand hdftA) :: Vector Complex Float :: Matrix Complex Float
+
+--       +         2.48361         +
+--       |                         |
+--       |- 0.265985 + 0.530898 %i |
+--       |                         |
+--       |- 0.257682 + 0.202979 %i |
+--       |                         |
+--       |- 0.256363 + 0.0580623 %i|
+--       |                         |
+--       |- 0.256363 - 0.0580623 %i|
+--       |                         |
+--       |- 0.257682 - 0.202979 %i |
+--       |                         |
+--       +- 0.265985 - 0.530898 %i +
+ 
+-- test  6
+ 
+(nagInverseDFT hdftA) :: Vector Float
+
+--       [0.34907,0.5489,0.74776,0.94459,1.1385,1.3285,1.5137]
+
+-- test  7
+
+dftC := nagDFT hseqC;
+dftC :: Vector Float
+ 
+-- [1.82616,1.86862,- 0.017503,0.502001,- 0.598725,- 0.0314404,- 2.62557]
+
+-- test  8
+
+(nagInverseDFT dftC) :: Vector Complex Float
+ 
+-- [0.34907, 0.5489 + 1.5137 %i, 0.74776 + 1.3285 %i, 0.94459 + 1.1385 %i,
+--  0.94459 - 1.1385 %i, 0.74776 - 1.3285 %i, 0.5489 - 1.5137 %i]
+
+-- test  9
+
+nagHermitianInverseDFT dftC
+ 
+-- [0.34907000000000005, 0.54889999999999983, 0.74775999999999987,
+--  0.94459000000000004, 1.1385000000000003, 1.3284999999999998,
+--  1.5136999999999998]
+
+-- test 10:
+
+dftsD := nagDFT seqsD;
+
+dftsD :: List Vector Complex Float
+ 
+-- [
+--   [1.07373, - 0.104062 - 0.00438406 %i, 0.112554 - 0.373777 %i, - 0.146684,
+--    0.112554 + 0.373777 %i, - 0.104062 + 0.00438406 %i]
+--   ,
+
+--   [1.39609, - 0.0365178 + 0.466584 %i, 0.077955 - 0.0607051 %i, - 0.152072,
+--    0.077955 + 0.0607051 %i, - 0.0365178 - 0.466584 %i]
+--   ,
+
+--   [1.12374, 0.0914068 - 0.050841 %i, 0.393551 + 0.345775 %i, 0.153011,
+--    0.393551 - 0.345775 %i, 0.0914068 + 0.050841 %i]
+--   ]
+
+-- test 11:
+
+invdftsD := nagInverseDFT dftsD ;
+invdftsD :: List Vector Complex Float
+ 
+-- [[0.3854,0.6772,0.1138,0.6751,0.6362,0.1424],
+--  [0.5417,0.2983,0.1181,0.7255,0.8638,0.8723],
+--  [0.9172,0.0644,0.6037,0.643,0.0428,0.4815]]
+
+-- test 12:
+
+dftsE := nagDFT seqsE;
+dftsE :: List Vector Float
+ 
+-- [[1.0788,0.662291,- 0.239146,- 0.578284,0.459192,- 0.438816],
+--  [0.857321,1.22614,0.353348,- 0.222169,0.341327,- 1.22908],
+--  [1.18245,0.262509,0.674406,0.552278,0.0539906,- 0.478963]]
+
+-- test 13:
+
+invdftsE := nagInverseDFT dftsE;
+invdftsE :: List Vector Complex Float
+ 
+-- [
+--   [0.3854, 0.6772 + 0.1424 %i, 0.1138 + 0.6362 %i, 0.6751,
+--    0.1138 - 0.6362 %i, 0.6772 - 0.1424 %i]
+--   ,
+
+--   [0.5417, 0.2983 + 0.8723 %i, 0.1181 + 0.8638 %i, 0.7255,
+--    0.1181 - 0.8638 %i, 0.2983 - 0.8723 %i]
+--   ,
+
+--   [0.9172, 0.0644 + 0.4815 %i, 0.6037 + 0.0428 %i, 0.643,
+--    0.6037 - 0.0428 %i, 0.0644 - 0.4815 %i]
+--   ]
+
+-- test 14:
+
+hdftsD := nagHermitianDFT seqsD;
+map(expand,hdftsD) :: List Vector Complex Float
+ 
+-- [
+--   [1.07373, - 0.104062 - 0.00438406 %i, 0.112554 - 0.373777 %i, - 0.146684,
+--    0.112554 + 0.373777 %i, - 0.104062 + 0.00438406 %i]
+--   ,
+
+--   [1.39609, - 0.0365178 + 0.466584 %i, 0.077955 - 0.0607051 %i, - 0.152072,
+--    0.077955 + 0.0607051 %i, - 0.0365178 - 0.466584 %i]
+--   ,
+
+--   [1.12374, 0.0914068 - 0.050841 %i, 0.393551 + 0.345775 %i, 0.153011,
+--    0.393551 - 0.345775 %i, 0.0914068 + 0.050841 %i]
+--   ]
+
+-- test 15:
+
+(nagInverseDFT hdftsD) :: List Vector Float
+ 
+-- [[0.3854,0.6772,0.1138,0.6751,0.6362,0.1424],
+--  [0.5417,0.2983,0.1181,0.7255,0.8638,0.8723],
+--  [0.9172,0.0644,0.6037,0.643,0.0428,0.4815]]
+
+-- test 16:
+
+dftsF := nagDFT seqsF;
+dftsF :: List Vector Complex Float
+ 
+-- [
+--   [1.07373 + 1.39609 %i, - 0.570647 - 0.0409019 %i, 0.173259 - 0.295822 %i,
+--    - 0.146684 - 0.152072 %i, 0.0518489 + 0.451732 %i,
+--    0.362522 - 0.0321337 %i]
+--   ,
+
+--   [1.12374 + 1.06765 %i, 0.172759 + 0.0385858 %i, 0.418548 + 0.748083 %i,
+--    0.153011 + 0.17522 %i, 0.368555 + 0.0565331 %i, 0.0100542 + 0.140268 %i]
+--   ,
+
+--   [0.909985 + 1.76167 %i, - 0.305418 + 0.0624335 %i,
+--    0.407884 - 0.0694786 %i, - 0.078547 + 0.0725049 %i,
+--    - 0.119334 + 0.128511 %i, - 0.531409 - 0.433531 %i]
+--   ]
+
+-- test 17:
+
+invdftsF := nagInverseDFT dftsF ;
+invdftsF :: List Vector Complex Float
+ 
+-- [
+--   [0.3854 + 0.5417 %i, 0.6772 + 0.2983 %i, 0.1138 + 0.1181 %i,
+--    0.6751 + 0.7255 %i, 0.6362 + 0.8638 %i, 0.1424 + 0.8723 %i]
+--   ,
+
+--   [0.9172 + 0.9089 %i, 0.0644 + 0.3118 %i, 0.6037 + 0.3465 %i,
+--    0.643 + 0.6198 %i, 0.0428 + 0.2668 %i, 0.4815 + 0.1614 %i]
+--   ,
+
+--   [0.1156 + 0.6214 %i, 0.0685 + 0.8681 %i, 0.206 + 0.706 %i,
+--    0.863 + 0.8652 %i, 0.6967 + 0.919 %i, 0.2792 + 0.3355 %i]
+--   ]
+
+-- test 18:
+
+nagHermitianInverseDFT dftsE
+ 
+-- [
+--   [0.38540000000000013, 0.67720000000000025, 0.11380000000000001,
+--    0.67510000000000014, 0.63620000000000021, 0.14240000000000003]
+--   ,
+
+--   [0.54170000000000018, 0.29830000000000012, 0.1181, 0.72550000000000014,
+--    0.86380000000000023, 0.87230000000000019]
+--   ,
+
+--   [0.91720000000000035, 0.064399999999999999, 0.60370000000000024,
+--    0.64300000000000013, 0.042799999999999991, 0.48150000000000015]
+--   ]
+
+-- error tests:
+
+-- test 19:
+
+nagDFT [vector [0.3854 + 0.5417*%i, 0.6772 + 0.2983*%i,   _
+                0.1138 + 0.1181*%i, 0.6751 + 0.7255*%i,   _
+                0.6362 + 0.8638*%i, 0.1424 + 0.8723*%i],  _
+        vector [0.1156 + 0.6214*%i, 0.0685 + 0.8681*%i,   _
+                0.6967 + 0.9190*%i, 0.2792 + 0.3355*%i]]
+ 
+-- Error signalled from user code:
+--    The data sequences in nagDFT must all have the same length. The 
+--    length of sequence 1 is 6 that of sequence 2 is 4.
+
+-- test 20:
+
+nagHermitianDFT [vector [0.3854, 0.6751, 0.6362, 0.1424], _
+                 vector [0.5417, 0.7255, 0.8638, 0.8723], _
+                 vector [0.9172, 0.0428, 0.4815]]
+ 
+-- Error signalled from user code:
+--    The data sequences in nagHermitianDFT must all have the same 
+--    length. The length of sequence 1 is 4 that of sequence 3 is 3.
+
+-- test 21:
+
+badSeqs : List PackedHermitianSequence DoubleFloat
+badSeqs := [pHS [0.3854, 0.1138, 0.6751, 0.6362, 0.1424],         _
+            pHS [0.5417, 0.2983, 0.1181, 0.7255, 0.8638, 0.8723], _
+            pHS [0.9172, 0.0644, 0.6037, 0.6430, 0.0428, 0.4815]];
+ 
+nagDFT badSeqs
+ 
+-- Error signalled from user code:
+--    The data sequences in nagDFT must all have the same length. The 
+--    length of sequence 1 is 5 that of sequence 2 is 6.
+
+outputGeneral()
+
+output "End of tests"
+
+#endif
+
+@
+\section{License}
+<<license>>=
+--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+--All rights reserved.
+--
+--Redistribution and use in source and binary forms, with or without
+--modification, are permitted provided that the following conditions are
+--met:
+--
+--    - Redistributions of source code must retain the above copyright
+--      notice, this list of conditions and the following disclaimer.
+--
+--    - Redistributions in binary form must reproduce the above copyright
+--      notice, this list of conditions and the following disclaimer in
+--      the documentation and/or other materials provided with the
+--      distribution.
+--
+--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
+--      names of its contributors may be used to endorse or promote products
+--      derived from this software without specific prior written permission.
+--
+--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+@
+<<*>>=
+<<license>>
+
+-- To test:
+-- sed -ne '1,/^#if NeverAssertThis/d;/#endif/d;p' < ndftip.as > ndftip.input
+-- axiom
+-- )set nag host <some machine running nagd>
+-- )r ndftip.input
+
+#unassert saturn
+
+#include "axiom.as"
+
+DF     ==> DoubleFloat ;
+CDF    ==> Complex DoubleFloat ;
+LDF    ==> List DoubleFloat ;
+LLDF   ==> List LDF ;
+VDF    ==> Vector DoubleFloat ;
+LVDF   ==> List VDF ;
+VCDF   ==> Vector Complex DoubleFloat ;
+LVCDF  ==> List VCDF ;
+MDF    ==> Matrix DoubleFloat ;
+MCDF   ==> Matrix Complex DoubleFloat ;
+INT    ==> Integer ;
+NNI    ==> NonNegativeInteger ;
+RSLT   ==> Result ;
+STRG   ==> String ;
+PHSDF  ==> PackedHermitianSequence DF;
+LPHSDF ==> List PackedHermitianSequence DF;
+
+<<NagDiscreteFourierTransformInterfacePackage>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}