From ab8cc85adde879fb963c94d15675783f2cf4b183 Mon Sep 17 00:00:00 2001 From: dos-reis Date: Tue, 14 Aug 2007 05:14:52 +0000 Subject: Initial population. --- src/algebra/ndftip.as.pamphlet | 1174 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1174 insertions(+) create mode 100644 src/algebra/ndftip.as.pamphlet (limited to 'src/algebra/ndftip.as.pamphlet') 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} +<>= ++++ 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} -- cgit v1.2.3