fricas
(1) -> <spad>
fricas
)abbrev package TYPEPKG TypePackage
TypePackage (T : Type) : Exports == Implementation where
Exports == with
typeof : T -> Type
Implementation == add
typeof(x:T) == T
fricas
)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 )
fricas
)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)
fricas
)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
fricas
)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>
fricas
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 SEC.
Cumulative Statistics for Constructor TypePackage
Time: 0 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 03 DEC 2024 07:19:37 PM):
; wrote /var/aw/var/LatexWiki/TYPEPKG.NRLIB/TYPEPKG.fasl
; compilation finished in 0:00:00.008
------------------------------------------------------------------------
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 03 DEC 2024 07:19:37 PM):
; wrote /var/aw/var/LatexWiki/CMPTYPE.NRLIB/CMPTYPE.fasl
; compilation finished in 0:00:00.004
------------------------------------------------------------------------
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 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 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: 1316)
;;; *** |RootOfUnity| REDEFINED
;;; *** |RootOfUnity| REDEFINED
Time: 0 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.03 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 03 DEC 2024 07:19:37 PM):
; wrote /var/aw/var/LatexWiki/ROU.NRLIB/ROU.fasl
; compilation finished in 0:00:00.012
------------------------------------------------------------------------
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 SEC.
compiling exported dftInvMatrix : () -> Matrix R
Time: 0 SEC.
compiling exported dft : Vector R -> Vector R
Time: 0 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.02 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 03 DEC 2024 07:19:37 PM):
; wrote /var/aw/var/LatexWiki/DFT.NRLIB/DFT.fasl
; compilation finished in 0:00:00.020
------------------------------------------------------------------------
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.02 SEC.
compiling exported udftMatrix : () -> Matrix Expression Complex Integer
Time: 0.01 SEC.
compiling exported udftInvMatrix : () -> Matrix Expression Complex Integer
Time: 0 SEC.
compiling exported udft : Vector Expression Complex Integer -> Vector Expression Complex Integer
Time: 0 SEC.
compiling exported iudft : Vector Expression Complex Integer -> Vector Expression Complex Integer
Time: 0 SEC.
(time taken in buildFunctor: 0)
;;; *** |UnitaryDiscreteFourierTransform| REDEFINED
;;; *** |UnitaryDiscreteFourierTransform| REDEFINED
Time: 0 SEC.
Cumulative Statistics for Constructor UnitaryDiscreteFourierTransform
Time: 0.06 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 03 DEC 2024 07:19:37 PM):
; wrote /var/aw/var/LatexWiki/UDFT.NRLIB/UDFT.fasl
; compilation finished in 0:00:00.008
------------------------------------------------------------------------
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§ion=3
--
X ==> EXPR COMPLEX INT
Type: Void
fricas
a:=convert(zerosOf(z^4-1).2)$ROU(4,X)
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
Type: Matrix(Expression(Complex(Integer)))
fricas
v:=vector [1::X,2-%i,-%i,-1+2*%i]
Type: Vector(Expression(Complex(Integer)))
fricas
w:=dft(v)$F
Type: Vector(Expression(Complex(Integer)))
fricas
test(w=m*v)
Type: Boolean
fricas
test(idft(w)$F=v)
Type: Boolean
fricas
determinant(m/2)
Type: Expression(Complex(Integer))
fricas
-- Examples
-----------
R:=IntegerMod 5
Type: Type
fricas
a:=2::R
fricas
n:=4
fricas
DFTZ5==>DFT(R,n,a)
Type: Void
fricas
dftMatrix()$DFTZ5
Type: Matrix(IntegerMod
?(5))
fricas
dftInvMatrix()$DFTZ5
Type: Matrix(IntegerMod
?(5))
fricas
dft([1,2,3,4])$DFTZ5
Type: Vector(IntegerMod
?(5))
fricas
idft([0,4,3,2])$DFTZ5
Type: Vector(IntegerMod
?(5))
fricas
-- R:=Expression Complex Integer
-- n:=3
-- a:=exp(-2*%i*%pi/n)
-- ... no!
M:=udftMatrix()$UDFT(5)
Type: Matrix(Expression(Complex(Integer)))
fricas
d:=determinant(M)
Type: Expression(Complex(Integer))
fricas
test(d=1) -- test false does not necessarily mean d<>1
Type: Boolean
fricas
test(d^2=1) -- true
Type: Boolean
fricas
r:=retract(d) --
fricas
test(r=1) -- true :)
Type: Boolean
fricas
)set mess time on
N:=20 -- 20 is quick, 50 ~12sec, 100 ?long
fricas
Time: 0 sec
M:=udftMatrix()$UDFT(N)
Type: Matrix(Expression(Complex(Integer)))
fricas
Time: 0.08 (EV) + 0.01 (OT) = 0.10 sec
IM:=udftInvMatrix()$UDFT(N)
Type: Matrix(Expression(Complex(Integer)))
fricas
Time: 0.04 (EV) + 0.01 (OT) = 0.05 sec
w:=M*vector([i for i in 1..N])
Type: Vector(Expression(Complex(Integer)))
fricas
Time: 0.01 (EV) = 0.02 sec
IM*w
Type: Vector(Expression(Complex(Integer)))
fricas
Time: 0.02 (EV) = 0.02 sec
E:=M*IM -- needs some time!
Type: Matrix(Expression(Complex(Integer)))
fricas
Time: 0.23 (EV) = 0.23 sec
trace(E) -- N
fricas
Time: 0 sec
--)synonym lispversion )lisp (lisp-implementation-version)
map(x+->x::Complex Float,m)
Type: Matrix(Complex(Float))
fricas
Time: 0 sec
fN:=complexNumeric(exp(-2*%pi*%i/N))
Type: Complex(Float)
fricas
Time: 0.02 sec
MF:=map(x+->subst(x,%z5=fN),M)
Type: Matrix(Expression(Complex(Float)))
fricas
Time: 0.16 (IN) + 0.03 (OT) = 0.20 sec
trace(MF) -- 1+i
Type: Complex(Float)
fricas
Time: 0 sec
determinant(MF) -- -i
Type: Expression(Complex(Float))
fricas
Time: 0 sec