|
|
last edited 5 years ago by pagani |
1 | ||
Editor: pagani
Time: 2018/07/16 20:34:29 GMT+0 |
||
Note: |
changed: - \begin{spad} )abbrev package TYPEPKG TypePackage TypePackage (T : Type) : Exports == Implementation where Exports == with typeof : T -> Type Implementation == add typeof(x:T) == T )abbrev package CMPTYPE CompareTypes CompareTypes(T:Type,S:Type):Exports == Implementation where SEX ==> SExpression Exports == with sameType? : (T, S) -> Boolean Implementation == add import TypePackage(T) import TypePackage(S) sameType?(x:T,y:S):Boolean == test ( typeof(x)::SEX = typeof(y):: SEX ) )abbrev domain ROU RootOfUnity ++ Author: Kurt Pagani ++ Date Created: Fri Jun 01 17:24:19 CEST 2018 ++ License: BSD ++ References: ++ https://en.wikipedia.org/wiki/Root_of_unity ++ https://en.wikipedia.org/wiki/Principal_root_of_unity ++ Description: ++ The nth roots of unity are, by definition, the roots of the polynomial ++ $P(z)=z^n−1$, and are therefore algebraic numbers. As the polynomial $P$ ++ is not irreducible - unless $n=1$, the primitive nth roots of unity are ++ roots of an irreducible polynomial of lower degree, called the cyclotomic ++ polynomial. ++ ++ Group of nth roots of unity ++ The product and the multiplicative inverse of two nth roots of unity ++ are also nth roots of unity. Therefore, the nth roots of unity form ++ a group under multiplication. ++ ++ Notes ++ Any algebraically closed field has exactly $n$ nth roots of unity if ++ $n$ is not divisible by the characteristic of the field. ++ ++ The significance of a root of unity being principal is that it is a ++ necessary condition for the theory of the discrete Fourier transform ++ to work out correctly. ++ ++ Usage and Examples ++ X ==> Expression Complex Integer ++ R ==> RootOfUnity(5,X) ++ z:X ++ r:=rootsOf(z^5-1) or zerosOf(z^5-1) or solve(z^5=1,'z) ++ q:=[convert(t)$R for t in r] ++ [primitive?(t) for t in q] ++ [principal?(t) for t in q] ++ RootOfUnity(n,R) : Exports == Implementation where n:PositiveInteger R:Ring CTOF ==> CoercibleTo OutputForm Exports == Join(Group,CTOF) with convert : R -> % ++ Convert r:R to a n-th root of unity if r^n=1$R. retract : % -> R ++ Retract a n-th root of unity to a member of R. 1 : () -> % ++ The ring unit. primitive? : % -> Boolean ++ An nth root of unity is primitive if it is not a kth root of unity ++ for some smaller k. principal? : % -> Boolean ++ A principal n-th root of unity of a ring is an element a:R ++ satisfying the equations a^n=1$R, sum(a^(j*k),j=0..n-1)=0 ++ for all 1<=k<n. coerce : % -> OutputForm ++ Coerce to output form. if R has ExpressionSpace then ExpressionSpace Implementation == R add Rep := R convert(x) == x^n = 1$R => x error "Probably not a root of unity." retract(x:%):R == x@Rep primitive?(x:%):Boolean == b:List Boolean:=[test(x^m=1$R) for m in 1..n-1] not reduce(_and,b) summ(a:R,m:PositiveInteger):R == s:List R:=[a^j for j in 0..m] reduce(_+,s) principal?(x:%):Boolean == n=1 => false a:R:=x@Rep nn:PositiveInteger:=(n-1)::PositiveInteger b:List Boolean:=[test(summ(a^k,nn)=0$R) for k in 1..nn] reduce(_and,b) )abbrev package DFT DiscreteFourierTransform ++ Author: Kurt Pagani ++ Date Created: Wed Mar 02 23:13:52 CET 2016 ++ License: BSD ++ References: ++ https://en.wikipedia.org/wiki/Discrete_Fourier_transform ++ https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general) ++ ++ Description: ++ Discrete Fourier transform (DFT) over any ring, commonly called ++ a number-theoretic transform (NTT) in the case of finite fields. ++ Given a ring R and a principal n-th root of unity a, the package ++ designator DFT(R,n,a) is used to compute various objects: ++ dftMatrix()$DFT(R,n,a) ++ for instance, will compute and display a matrix M (the DFT matrix) ++ such the discrete Fourier transform of the vector v is given by M*v. ++ ++ Usage and Examples ++ X ==> EXPR COMPLEX INT ++ a:=convert(zerosOf(z^5-1).2)$ROU(5,X) ++ F ==> DFT(X,5,a) ++ m:=dftMatrix()$F ++ v:=vector [1::X,2,3,4,5] ++ w:=dft(v)$F ++ test(w=m*v) ++ test(idft(w)$F=v) ++ more? see test_dft.input ++ DiscreteFourierTransform(R,n,a) : Exports == Implementation where R:Ring n:PositiveInteger a:RootOfUnity(n,R) MATR ==> Matrix R Exports == with dftMatrix : () -> MATR ++ dftMatrix() computes and returns the DFT(R,n,a) matrix. dftInvMatrix : () -> MATR ++ dftInvMatrix() computes and returns the inverse DFT(R,n,a) matrix. dft : Vector R -> Vector R ++ dft(v) performs the DFT of the vector v. idft : Vector R -> Vector R ++ idft(v) performs the inverse DFT of the vector v. Implementation == add err1 := "The n-th root of unity provided is not principal." err2 := "Not invertible." err3 := "Wrong dimension of vector." not principal?(a) => error err1 dftMatrix() == w:R:=retract(a) m:MATR:=matrix [[w^(i*j) for i in 0..n-1] for j in 0..n-1] dftInvMatrix() == w:R:=retract(a) b:=recip(w) c:=recip(n*1::R) if b case R and c case R then m:MATR:=c * matrix [[b^(i*j) for i in 0..n-1] for j in 0..n-1] else error err2 dft(v) == #v ~= n => error err3 dftMatrix() * v idft(v) == #v ~= n => error err3 dftInvMatrix() * v )abbrev package UDFT UnitaryDiscreteFourierTransform ++ Description: ++ With unitary normalization constants 1/sqrt(n), the DFT becomes a ++ unitary transformation, defined by a unitary matrix: ++ UDFT = DFT/sqrt(n) ++ Since IDFT = DFT*/N, we have IUDFT = sqrt(n) * IDFT. ++ Moreover, |det(UDFT)|=1, IUDFT=UDFT*. ++ ++ Usage and Examples ++ M:=udftMatrix()$UDFT(100) ++ IM:=udftInvMatrix()$UDFT(100) ++ w:=M*vector([i for i in 1..100]) ++ IM*w ++ E:=M*IM -- needs some time! ++ trace(E) -- 100 ++ ++ Note that udft/iudft also compute udftMatrix internally each time ++ the functions are called, so it is certainly better to compute the ++ matrices once, if several transformations are to be performed. ++ UnitaryDiscreteFourierTransform(n) : Exports == Implementation where n:PositiveInteger R ==> Expression Complex Integer DFT ==> DiscreteFourierTransform MATR ==> Matrix R Exports == with principalNthRootOfUnity : () -> RootOfUnity(n,R) ++ principalNthRootOfUnity() returns the principal n-th root ++ of unity which generates the cyclic group of the primitive ++ roots of unity as well as the Vandermonde matrix UDFT for ++ $\mathbb{C}^n$ ~ Expression Complex Integer. udftMatrix : () -> MATR ++ udftMatrix() computes and returns the UDFT(n) matrix. udftInvMatrix : () -> MATR ++ udftInvMatrix() computes and returns the inverse UDFT(n) matrix. udft : Vector R -> Vector R ++ udft(v) performs the UDFT of the vector v. iudft : Vector R -> Vector R ++ iudft(v) performs the inverse UDFT of the vector v. Implementation == add err1 := "Wrong dimension of vector." principalNthRootOfUnity():RootOfUnity(n,R) == z:R:='z::R p:R:=z^n-1$R zop:List R:=zerosOf(p) a:R:=zop.2 convert(inv a)$RootOfUnity(n,R) a:RootOfUnity(n,R):=principalNthRootOfUnity() udftMatrix() == dftMatrix()$DFT(R,n,a) / sqrt(n::R) udftInvMatrix() == dftInvMatrix()$DFT(R,n,a) * sqrt(n::R) udft(v) == #v ~= n => error err1 udftMatrix() * v iudft(v) == #v ~= n => error err1 udftInvMatrix() * v \end{spad} Test flavours \begin{axiom} --)co DFT -- -- https://en.wikipedia.org/w/index.php?title= -- Discrete_Fourier_transform&action=edit§ion=3 -- X ==> EXPR COMPLEX INT a:=convert(zerosOf(z^4-1).2)$ROU(4,X) F ==> DFT(X,4,inv a) -- Note: inv(a) gives another result than a !!! m:=dftMatrix()$F v:=vector [1::X,2-%i,-%i,-1+2*%i] w:=dft(v)$F test(w=m*v) test(idft(w)$F=v) determinant(m/2) -- Examples ----------- R:=IntegerMod 5 a:=2::R n:=4 DFTZ5==>DFT(R,n,a) dftMatrix()$DFTZ5 dftInvMatrix()$DFTZ5 dft([1,2,3,4])$DFTZ5 idft([0,4,3,2])$DFTZ5 -- R:=Expression Complex Integer -- n:=3 -- a:=exp(-2*%i*%pi/n) -- ... no! M:=udftMatrix()$UDFT(5) d:=determinant(M) test(d=1) -- test false does not necessarily mean d<>1 test(d^2=1) -- true r:=retract(d) -- Type: AlgebraicNumber test(r=1) -- true :) )set mess time on N:=20 -- 20 is quick, 50 ~12sec, 100 ?long M:=udftMatrix()$UDFT(N) IM:=udftInvMatrix()$UDFT(N) w:=M*vector([i for i in 1..N]) IM*w E:=M*IM -- needs some time! trace(E) -- N --)synonym lispversion )lisp (lisp-implementation-version) map(x+->x::Complex Float,m) fN:=complexNumeric(exp(-2*%pi*%i/N)) MF:=map(x+->subst(x,%z5=fN),M) trace(MF) -- 1+i determinant(MF) -- -i \end{axiom}
)abbrev package TYPEPKG TypePackage TypePackage (T : Type) : Exports == Implementation where Exports == with typeof : T -> Type Implementation == add typeof(x:T) == T
)abbrev package CMPTYPE CompareTypes CompareTypes(T:Type,S:Type):Exports == Implementation where SEX ==> SExpression Exports == with sameType? : (T, S) -> Boolean Implementation == add
import TypePackage(T) import TypePackage(S)
sameType?(x:T,y:S):Boolean == test ( typeof(x)::SEX = typeof(y):: SEX )
)abbrev domain ROU RootOfUnity ++ Author: Kurt Pagani ++ Date Created: Fri Jun 01 17:24:19 CEST 2018 ++ License: BSD ++ References: ++ https://en.wikipedia.org/wiki/Root_of_unity ++ https://en.wikipedia.org/wiki/Principal_root_of_unity ++ Description: ++ The nth roots of unity are,by definition, the roots of the polynomial ++ $P(z)=z^n−1$, and are therefore algebraic numbers. As the polynomial $P$ ++ is not irreducible - unless $n=1$, the primitive nth roots of unity are ++ roots of an irreducible polynomial of lower degree, called the cyclotomic ++ polynomial. ++ ++ Group of nth roots of unity ++ The product and the multiplicative inverse of two nth roots of unity ++ are also nth roots of unity. Therefore, the nth roots of unity form ++ a group under multiplication. ++ ++ Notes ++ Any algebraically closed field has exactly $n$ nth roots of unity if ++ $n$ is not divisible by the characteristic of the field. ++ ++ The significance of a root of unity being principal is that it is a ++ necessary condition for the theory of the discrete Fourier transform ++ to work out correctly. ++ ++ Usage and Examples ++ X ==> Expression Complex Integer ++ R ==> RootOfUnity(5, X) ++ z:X ++ r:=rootsOf(z^5-1) or zerosOf(z^5-1) or solve(z^5=1, 'z) ++ q:=[convert(t)$R for t in r] ++ [primitive?(t) for t in q] ++ [principal?(t) for t in q] ++ RootOfUnity(n, R) : Exports == Implementation where
n:PositiveInteger R:Ring
CTOF ==> CoercibleTo OutputForm
Exports == Join(Group,CTOF) with
convert : R -> % ++ Convert r:R to a n-th root of unity if r^n=1$R. retract : % -> R ++ Retract a n-th root of unity to a member of R. 1 : () -> % ++ The ring unit. primitive? : % -> Boolean ++ An nth root of unity is primitive if it is not a kth root of unity ++ for some smaller k. principal? : % -> Boolean ++ A principal n-th root of unity of a ring is an element a:R ++ satisfying the equations a^n=1$R,sum(a^(j*k), j=0..n-1)=0 ++ for all 1<=k<n. coerce : % -> OutputForm ++ Coerce to output form.
if R has ExpressionSpace then ExpressionSpace
Implementation == R add
Rep := R
convert(x) == x^n = 1$R => x error "Probably not a root of unity."
retract(x:%):R == x@Rep
primitive?(x:%):Boolean == b:List Boolean:=[test(x^m=1$R) for m in 1..n-1] not reduce(_and,b)
summ(a:R,m:PositiveInteger):R == s:List R:=[a^j for j in 0..m] reduce(_+, s)
principal?(x:%):Boolean == n=1 => false a:R:=x@Rep nn:PositiveInteger:=(n-1)::PositiveInteger b:List Boolean:=[test(summ(a^k,nn)=0$R) for k in 1..nn] reduce(_and, b)
)abbrev package DFT DiscreteFourierTransform ++ Author: Kurt Pagani ++ Date Created: Wed Mar 02 23:13:52 CET 2016 ++ License: BSD ++ References: ++ https://en.wikipedia.org/wiki/Discrete_Fourier_transform ++ https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general) ++ ++ Description: ++ Discrete Fourier transform (DFT) over any ring,commonly called ++ a number-theoretic transform (NTT) in the case of finite fields. ++ Given a ring R and a principal n-th root of unity a, the package ++ designator DFT(R, n, a) is used to compute various objects: ++ dftMatrix()$DFT(R, n, a) ++ for instance, will compute and display a matrix M (the DFT matrix) ++ such the discrete Fourier transform of the vector v is given by M*v. ++ ++ Usage and Examples ++ X ==> EXPR COMPLEX INT ++ a:=convert(zerosOf(z^5-1).2)$ROU(5, X) ++ F ==> DFT(X, 5, a) ++ m:=dftMatrix()$F ++ v:=vector [1::X, 2, 3, 4, 5] ++ w:=dft(v)$F ++ test(w=m*v) ++ test(idft(w)$F=v) ++ more? see test_dft.input ++ DiscreteFourierTransform(R, n, a) : Exports == Implementation where
R:Ring n:PositiveInteger a:RootOfUnity(n,R)
MATR ==> Matrix R
Exports == with
dftMatrix : () -> MATR ++ dftMatrix() computes and returns the DFT(R,n, a) matrix. dftInvMatrix : () -> MATR ++ dftInvMatrix() computes and returns the inverse DFT(R, n, a) matrix. dft : Vector R -> Vector R ++ dft(v) performs the DFT of the vector v. idft : Vector R -> Vector R ++ idft(v) performs the inverse DFT of the vector v.
Implementation == add
err1 := "The n-th root of unity provided is not principal." err2 := "Not invertible." err3 := "Wrong dimension of vector."
not principal?(a) => error err1
dftMatrix() == w:R:=retract(a) m:MATR:=matrix [[w^(i*j) for i in 0..n-1] for j in 0..n-1]
dftInvMatrix() == w:R:=retract(a) b:=recip(w) c:=recip(n*1::R) if b case R and c case R then m:MATR:=c * matrix [[b^(i*j) for i in 0..n-1] for j in 0..n-1] else error err2
dft(v) == #v ~= n => error err3 dftMatrix() * v
idft(v) == #v ~= n => error err3 dftInvMatrix() * v
)abbrev package UDFT UnitaryDiscreteFourierTransform ++ Description: ++ With unitary normalization constants 1/sqrt(n),the DFT becomes a ++ unitary transformation, defined by a unitary matrix: ++ UDFT = DFT/sqrt(n) ++ Since IDFT = DFT*/N, we have IUDFT = sqrt(n) * IDFT. ++ Moreover, |det(UDFT)|=1, IUDFT=UDFT*. ++ ++ Usage and Examples ++ M:=udftMatrix()$UDFT(100) ++ IM:=udftInvMatrix()$UDFT(100) ++ w:=M*vector([i for i in 1..100]) ++ IM*w ++ E:=M*IM -- needs some time! ++ trace(E) -- 100 ++ ++ Note that udft/iudft also compute udftMatrix internally each time ++ the functions are called, so it is certainly better to compute the ++ matrices once, if several transformations are to be performed. ++ UnitaryDiscreteFourierTransform(n) : Exports == Implementation where
n:PositiveInteger R ==> Expression Complex Integer DFT ==> DiscreteFourierTransform MATR ==> Matrix R
Exports == with
principalNthRootOfUnity : () -> RootOfUnity(n,R) ++ principalNthRootOfUnity() returns the principal n-th root ++ of unity which generates the cyclic group of the primitive ++ roots of unity as well as the Vandermonde matrix UDFT for ++ $\mathbb{C}^n$ ~ Expression Complex Integer. udftMatrix : () -> MATR ++ udftMatrix() computes and returns the UDFT(n) matrix. udftInvMatrix : () -> MATR ++ udftInvMatrix() computes and returns the inverse UDFT(n) matrix. udft : Vector R -> Vector R ++ udft(v) performs the UDFT of the vector v. iudft : Vector R -> Vector R ++ iudft(v) performs the inverse UDFT of the vector v.
Implementation == add
err1 := "Wrong dimension of vector."
principalNthRootOfUnity():RootOfUnity(n,R) == z:R:='z::R p:R:=z^n-1$R zop:List R:=zerosOf(p) a:R:=zop.2 convert(inv a)$RootOfUnity(n, R)
a:RootOfUnity(n,R):=principalNthRootOfUnity()
udftMatrix() == dftMatrix()$DFT(R,n, a) / sqrt(n::R)
udftInvMatrix() == dftInvMatrix()$DFT(R,n, a) * sqrt(n::R)
udft(v) == #v ~= n => error err1 udftMatrix() * v
iudft(v) == #v ~= n => error err1 udftInvMatrix() * v
Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/6294884381899805856-25px001.spad using old system compiler. TYPEPKG abbreviates package TypePackage ------------------------------------------------------------------------ initializing NRLIB TYPEPKG for TypePackage compiling into NRLIB TYPEPKG compiling exported typeof : T$ -> Type Time: 0 SEC.
(time taken in buildFunctor: 0)
;;; *** |TypePackage| REDEFINED
;;; *** |TypePackage| REDEFINED Time: 0.01 SEC.
Cumulative Statistics for Constructor TypePackage Time: 0.01 seconds
finalizing NRLIB TYPEPKG Processing TypePackage for Browser database: --->-->TypePackage(constructor): Not documented!!!! --->-->TypePackage((typeof ((Type) T$))): Not documented!!!! --->-->TypePackage(): Missing Description ; compiling file "/var/aw/var/LatexWiki/TYPEPKG.NRLIB/TYPEPKG.lsp" (written 16 JUL 2018 08:34:29 PM):
; /var/aw/var/LatexWiki/TYPEPKG.NRLIB/TYPEPKG.fasl written ; compilation finished in 0:00:00.011 ------------------------------------------------------------------------ TypePackage is now explicitly exposed in frame initial TypePackage will be automatically loaded when needed from /var/aw/var/LatexWiki/TYPEPKG.NRLIB/TYPEPKG
CMPTYPE abbreviates package CompareTypes ------------------------------------------------------------------------ initializing NRLIB CMPTYPE for CompareTypes compiling into NRLIB CMPTYPE importing TypePackage T$ importing TypePackage S compiling exported sameType? : (T$,S) -> Boolean Time: 0 SEC.
(time taken in buildFunctor: 0)
;;; *** |CompareTypes| REDEFINED
;;; *** |CompareTypes| REDEFINED Time: 0 SEC.
Cumulative Statistics for Constructor CompareTypes Time: 0 seconds
finalizing NRLIB CMPTYPE Processing CompareTypes for Browser database: --->-->CompareTypes(constructor): Not documented!!!! --->-->CompareTypes((sameType? ((Boolean) T$ S))): Not documented!!!! --->-->CompareTypes(): Missing Description ; compiling file "/var/aw/var/LatexWiki/CMPTYPE.NRLIB/CMPTYPE.lsp" (written 16 JUL 2018 08:34:29 PM):
; /var/aw/var/LatexWiki/CMPTYPE.NRLIB/CMPTYPE.fasl written ; compilation finished in 0:00:00.012 ------------------------------------------------------------------------ CompareTypes is now explicitly exposed in frame initial CompareTypes will be automatically loaded when needed from /var/aw/var/LatexWiki/CMPTYPE.NRLIB/CMPTYPE
ROU abbreviates domain RootOfUnity ------------------------------------------------------------------------ initializing NRLIB ROU for RootOfUnity compiling into NRLIB ROU compiling exported convert : R -> $ Time: 0.01 SEC.
compiling exported retract : $ -> R ROU;retract;$R;2 is replaced by x Time: 0 SEC.
compiling exported primitive? : $ -> Boolean Time: 0.02 SEC.
compiling local summ : (R,PositiveInteger) -> R Time: 0.01 SEC.
compiling exported principal? : $ -> Boolean Time: 0 SEC.
****** Domain: $ already in scope augmenting $: (RetractableTo (Integer)) ****** Domain: R already in scope augmenting R: (ExpressionSpace) ****** Domain: $ already in scope augmenting $: (Ring) ****** Domain: R already in scope augmenting R: (ExpressionSpace) ****** Domain: R already in scope augmenting R: (ExpressionSpace) (time taken in buildFunctor: 0)
;;; *** |RootOfUnity| REDEFINED
;;; *** |RootOfUnity| REDEFINED Time: 0.01 SEC.
Warnings: [1] not known that (Comparable) is of mode (CATEGORY domain (SIGNATURE convert ($ R)) (SIGNATURE retract (R $)) (SIGNATURE (One) ($)) (SIGNATURE primitive? ((Boolean) $)) (SIGNATURE principal? ((Boolean) $)) (SIGNATURE coerce ((OutputForm) $)) (IF (has R (ExpressionSpace)) (ATTRIBUTE (ExpressionSpace)) noBranch))
Cumulative Statistics for Constructor RootOfUnity Time: 0.05 seconds
finalizing NRLIB ROU Processing RootOfUnity for Browser database: --------constructor--------- --------(convert (% R))--------- --->-->RootOfUnity((convert (% R))): Improper first word in comments: Convert "Convert \\spad{r:R} to a \\spad{n}-th root of unity if \\spad{r^n=1}\\$\\spad{R}." --------(retract (R %))--------- --->-->RootOfUnity((retract (R %))): Improper first word in comments: Retract "Retract a \\spad{n}-th root of unity to a member of \\spad{R}." --------((One) (%))--------- --->-->RootOfUnity(((One) (%))): Improper first word in comments: The "The ring unit." --------(primitive? ((Boolean) %))--------- --->-->RootOfUnity((primitive? ((Boolean) %))): Improper first word in comments: An "An \\spad{n}th root of unity is primitive if it is not a \\spad{k}th root of unity for some smaller \\spad{k}." --------(principal? ((Boolean) %))--------- --->-->RootOfUnity((principal? ((Boolean) %))): Improper first word in comments: A "A principal \\spad{n}-th root of unity of a ring is an element a:R satisfying the equations \\spad{a^n=1}\\$\\spad{R},{} sum(a^(\\spad{j*k}), {}\\spad{j=0}..\\spad{n}-1)\\spad{=0} for all 1<=k<n." --------(coerce ((OutputForm) %))--------- --->-->RootOfUnity((coerce ((OutputForm) %))): Improper first word in comments: Coerce "Coerce to output form." ; compiling file "/var/aw/var/LatexWiki/ROU.NRLIB/ROU.lsp" (written 16 JUL 2018 08:34:30 PM):
; /var/aw/var/LatexWiki/ROU.NRLIB/ROU.fasl written ; compilation finished in 0:00:00.040 ------------------------------------------------------------------------ RootOfUnity is now explicitly exposed in frame initial RootOfUnity will be automatically loaded when needed from /var/aw/var/LatexWiki/ROU.NRLIB/ROU
DFT abbreviates package DiscreteFourierTransform ------------------------------------------------------------------------ initializing NRLIB DFT for DiscreteFourierTransform compiling into NRLIB DFT compiling exported dftMatrix : () -> Matrix R Time: 0.02 SEC.
compiling exported dftInvMatrix : () -> Matrix R Time: 0.01 SEC.
compiling exported dft : Vector R -> Vector R Time: 0.01 SEC.
compiling exported idft : Vector R -> Vector R Time: 0 SEC.
(time taken in buildFunctor: 0)
;;; *** |DiscreteFourierTransform| REDEFINED
;;; *** |DiscreteFourierTransform| REDEFINED Time: 0 SEC.
Warnings: [1] unknown Functor code (error (QREFELT $ 9))
Cumulative Statistics for Constructor DiscreteFourierTransform Time: 0.04 seconds
finalizing NRLIB DFT Processing DiscreteFourierTransform for Browser database: --------constructor--------- --------(dftMatrix ((Matrix R)))--------- --------(dftInvMatrix ((Matrix R)))--------- --------(dft ((Vector R) (Vector R)))--------- --------(idft ((Vector R) (Vector R)))--------- ; compiling file "/var/aw/var/LatexWiki/DFT.NRLIB/DFT.lsp" (written 16 JUL 2018 08:34:30 PM):
; /var/aw/var/LatexWiki/DFT.NRLIB/DFT.fasl written ; compilation finished in 0:00:00.044 ------------------------------------------------------------------------ DiscreteFourierTransform is now explicitly exposed in frame initial DiscreteFourierTransform will be automatically loaded when needed from /var/aw/var/LatexWiki/DFT.NRLIB/DFT
UDFT abbreviates package UnitaryDiscreteFourierTransform ------------------------------------------------------------------------ initializing NRLIB UDFT for UnitaryDiscreteFourierTransform compiling into NRLIB UDFT Local variable err1 type redefined: (Wrong dimension of vector.) to (The n-th root of unity provided is not principal.) compiling exported principalNthRootOfUnity : () -> RootOfUnity(n,Expression Complex Integer) Time: 0.06 SEC.
compiling exported udftMatrix : () -> Matrix Expression Complex Integer Time: 0.03 SEC.
compiling exported udftInvMatrix : () -> Matrix Expression Complex Integer Time: 0.02 SEC.
compiling exported udft : Vector Expression Complex Integer -> Vector Expression Complex Integer Time: 0.01 SEC.
compiling exported iudft : Vector Expression Complex Integer -> Vector Expression Complex Integer Time: 0.02 SEC.
(time taken in buildFunctor: 0)
;;; *** |UnitaryDiscreteFourierTransform| REDEFINED
;;; *** |UnitaryDiscreteFourierTransform| REDEFINED Time: 0 SEC.
Cumulative Statistics for Constructor UnitaryDiscreteFourierTransform Time: 0.14 seconds
finalizing NRLIB UDFT Processing UnitaryDiscreteFourierTransform for Browser database: --------constructor--------- --------(principalNthRootOfUnity ((RootOfUnity n (Expression (Complex (Integer))))))--------- --->-->UnitaryDiscreteFourierTransform((principalNthRootOfUnity ((RootOfUnity n (Expression (Complex (Integer))))))): Unexpected HT command: \mathbb "\\spad{principalNthRootOfUnity()} returns the principal \\spad{n}-th root of unity which generates the cyclic group of the primitive roots of unity as well as the Vandermonde matrix UDFT for \\$\\mathbb{\\spad{C}}\\spad{^n}\\$ ~ Expression Complex Integer." --------(udftMatrix ((Matrix (Expression (Complex (Integer))))))--------- --------(udftInvMatrix ((Matrix (Expression (Complex (Integer))))))--------- --------(udft ((Vector (Expression (Complex (Integer)))) (Vector (Expression (Complex (Integer))))))--------- --------(iudft ((Vector (Expression (Complex (Integer)))) (Vector (Expression (Complex (Integer))))))--------- ; compiling file "/var/aw/var/LatexWiki/UDFT.NRLIB/UDFT.lsp" (written 16 JUL 2018 08:34:30 PM):
; /var/aw/var/LatexWiki/UDFT.NRLIB/UDFT.fasl written ; compilation finished in 0:00:00.025 ------------------------------------------------------------------------ UnitaryDiscreteFourierTransform is now explicitly exposed in frame initial UnitaryDiscreteFourierTransform will be automatically loaded when needed from /var/aw/var/LatexWiki/UDFT.NRLIB/UDFT
Test flavours
--)co DFT
a:=convert(zerosOf(z^4-1).2)$ROU(4,X)
(1) |
F ==> DFT(X,4, inv a) -- Note: inv(a) gives another result than a !!!
m:=dftMatrix()$F
(2) |
v:=vector [1::X,2-%i, -%i, -1+2*%i]
(3) |
w:=dft(v)$F
(4) |
test(w=m*v)
(5) |
test(idft(w)$F=v)
(6) |
determinant(m/2)
(7) |
-- Examples ----------- R:=IntegerMod 5
(8) |
a:=2::R
(9) |
n:=4
(10) |
DFTZ5==>DFT(R,n, a)
dftMatrix()$DFTZ5
(11) |
dftInvMatrix()$DFTZ5
(12) |
dft([1,2, 3, 4])$DFTZ5
(13) |
idft([0,4, 3, 2])$DFTZ5
(14) |
-- R:=Expression Complex Integer -- n:=3 -- a:=exp(-2*%i*%pi/n) -- ... no!
M:=udftMatrix()$UDFT(5)
(15) |
d:=determinant(M)
(16) |
test(d=1) -- test false does not necessarily mean d<>1
(17) |
test(d^2=1) -- true
(18) |
r:=retract(d) --
(19) |
test(r=1) -- true :)
(20) |
)set mess time on
N:=20 -- 20 is quick,50 ~12sec, 100 ?long
(21) |
Time: 0 sec M:=udftMatrix()$UDFT(N)
(22) |
Time: 0.12 (EV) + 0.01 (OT) = 0.13 sec IM:=udftInvMatrix()$UDFT(N)
(23) |
Time: 0.06 (EV) = 0.06 sec w:=M*vector([i for i in 1..N])
(24) |
Time: 0.01 (IN) + 0.02 (EV) = 0.03 sec IM*w
(25) |
Time: 0.03 (EV) = 0.03 sec E:=M*IM -- needs some time!
(26) |
Time: 0.46 (EV) = 0.46 sec trace(E) -- N
(27) |
Time: 0 sec
--)synonym lispversion )lisp (lisp-implementation-version) map(x+->x::Complex Float,m)
(28) |
Time: 0.01 (EV) = 0.01 sec
fN:=complexNumeric(exp(-2*%pi*%i/N))
(29) |
Time: 0.02 (OT) = 0.02 sec MF:=map(x+->subst(x,%z5=fN), M)
(30) |
Time: 0.68 (IN) + 0.02 (EV) + 0.07 (OT) = 0.77 sec trace(MF) -- 1+i
(31) |
Time: 0.01 (IN) = 0.01 sec determinant(MF) -- -i
(32) |
Time: 0.02 (EV) = 0.02 sec