\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}