)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

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