login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

# Edit detail for SandBoxDFT revision 1 of 1

 1 Editor: pagani Time: 2018/07/16 20:34:29 GMT+0 Note:

changed:
-
)abbrev package TYPEPKG TypePackage
TypePackage (T : Type) : Exports == Implementation where
Exports == with
typeof : T -> Type
typeof(x:T) == T

)abbrev package CMPTYPE CompareTypes
CompareTypes(T:Type,S:Type):Exports == Implementation  where
SEX ==> SExpression
Exports == with
sameType? : (T, S) -> Boolean

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
++ 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
++ 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&section=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
typeof(x:T) == T
)abbrev package CMPTYPE CompareTypes
CompareTypes(T:Type,S:Type):Exports == Implementation  where
SEX ==> SExpression
Exports == with
sameType? : (T, S) -> Boolean
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
++ 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
++ 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 spad  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:  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:  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

fricas
--)co DFT
-- https://en.wikipedia.org/w/index.php?title= -- Discrete_Fourier_transform&action=edit&section=3 -- X ==> EXPR COMPLEX INT
Type: Void
fricas
a:=convert(zerosOf(z^4-1).2)$ROU(4,X) (1) Type: RootOfUnity?(4,Expression(Complex(Integer))) fricas F ==> DFT(X,4,inv a) -- Note: inv(a) gives another result than a !!! Type: Void fricas m:=dftMatrix()$F (2)
Type: Matrix(Expression(Complex(Integer)))
fricas
v:=vector [1::X,2-%i,-%i,-1+2*%i] (3)
Type: Vector(Expression(Complex(Integer)))
fricas
w:=dft(v)$F (4) Type: Vector(Expression(Complex(Integer))) fricas test(w=m*v) (5) Type: Boolean fricas test(idft(w)$F=v) (6)
Type: Boolean
fricas
determinant(m/2) (7)
Type: Expression(Complex(Integer))
fricas
-- Examples
-----------
R:=IntegerMod 5 (8)
Type: Type
fricas
a:=2::R (9)
Type: IntegerMod?(5)
fricas
n:=4 (10)
Type: PositiveInteger?
fricas
DFTZ5==>DFT(R,n,a)
Type: Void
fricas
dftMatrix()$DFTZ5 (11) Type: Matrix(IntegerMod?(5)) fricas dftInvMatrix()$DFTZ5 (12)
Type: Matrix(IntegerMod?(5))
fricas
dft([1,2,3,4])$DFTZ5 (13) Type: Vector(IntegerMod?(5)) fricas idft([0,4,3,2])$DFTZ5 (14)
Type: Vector(IntegerMod?(5))
fricas
-- R:=Expression Complex Integer
-- n:=3
-- a:=exp(-2*%i*%pi/n)
-- ... no!
M:=udftMatrix()$UDFT(5) (15) Type: Matrix(Expression(Complex(Integer))) fricas d:=determinant(M) (16) Type: Expression(Complex(Integer)) fricas test(d=1) -- test false does not necessarily mean d<>1 (17) Type: Boolean fricas test(d^2=1) -- true (18) Type: Boolean fricas r:=retract(d) -- Type: AlgebraicNumber? (19) Type: AlgebraicNumber? fricas test(r=1) -- true :) (20) Type: Boolean fricas )set mess time on N:=20 -- 20 is quick, 50 ~12sec, 100 ?long (21) Type: PositiveInteger? fricas Time: 0 sec M:=udftMatrix()$UDFT(N) (22)
Type: Matrix(Expression(Complex(Integer)))
fricas
Time: 0.12 (EV) + 0.01 (OT) = 0.13 sec
IM:=udftInvMatrix()\$UDFT(N) (23)
Type: Matrix(Expression(Complex(Integer)))
fricas
Time: 0.06 (EV) = 0.06 sec
w:=M*vector([i for i in 1..N]) (24)
Type: Vector(Expression(Complex(Integer)))
fricas
Time: 0.01 (IN) + 0.02 (EV) = 0.03 sec
IM*w (25)
Type: Vector(Expression(Complex(Integer)))
fricas
Time: 0.03 (EV) = 0.03 sec
E:=M*IM -- needs some time! (26)
Type: Matrix(Expression(Complex(Integer)))
fricas
Time: 0.46 (EV) = 0.46 sec
trace(E) -- N (27)
Type: PositiveInteger?
fricas
Time: 0 sec
--)synonym lispversion )lisp (lisp-implementation-version)
map(x+->x::Complex Float,m) (28)
Type: Matrix(Complex(Float))
fricas
Time: 0.01 (EV) = 0.01 sec
fN:=complexNumeric(exp(-2*%pi*%i/N)) (29)
Type: Complex(Float)
fricas
Time: 0.02 (OT) = 0.02 sec
MF:=map(x+->subst(x,%z5=fN),M) (30)
Type: Matrix(Expression(Complex(Float)))
fricas
Time: 0.68 (IN) + 0.02 (EV) + 0.07 (OT) = 0.77 sec
trace(MF) -- 1+i (31)
Type: Complex(Float)
fricas
Time: 0.01 (IN) = 0.01 sec
determinant(MF) -- -i (32)
Type: Expression(Complex(Float))
fricas
Time: 0.02 (EV) = 0.02 sec