\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} <>= +++ 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} <>= --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. @ <<*>>= <> -- To test: -- sed -ne '1,/^#if NeverAssertThis/d;/#endif/d;p' < ndftip.as > ndftip.input -- axiom -- )set nag host -- )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; <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}