fast fourier transform using only integers on HVM

FFT.hvm

(Rot (V z))   = (V (- 0.0 z))
(Rot (G x y)) = (G (Rot y) x)

(Add (V z)   (V w))   = (V (+ z w))
(Add (G x y) (G w z)) = (G (Add x w) (Add y z))

(Get (V x)   f) = (f x)
(Get (G x y) f) = (f x y)

Nil     = λm λx (m x)
(Skp f) = λm (f λa (Get a λx λy (G (m x) (m y))))
(Twx f) = λm (f λa (Get a λx λy (Rot (G (m x) (m y)))))

(Mix (L x) (L y) twid zero) =
  let mins = (Add x (twid λa(a)     y))
  let maxs = (Add x (twid λa(Rot a) y))
  (B (L (G mins zero)) (L (G maxs zero)))
(Mix (B x y) (B w z) twid zero) =
  let lfts = (Mix x w (Skp twid) (G zero zero))
  let rgts = (Mix y z (Twx twid) (G zero zero))
  (B lfts rgts)

(FFT (L x)) = (L (V x))
(FFT (B x y)) =
  let eves = (FFT x)
  let odds = (FFT y)
  (Mix eves odds Nil (V 0.0))

// Conversions
// -----------

// List

(List.concat List.nil         ys) = ys
(List.concat (List.cons x xs) ys) = (List.cons x (List.concat xs ys))

// Complex

PI = 3.141592653589793

(CSca (C ar ai) s)         = (C (* ar s) (* ai s))
(CAdd (C ar ai) (C br bi)) = (C (+ ar br) (+ ai bi))
(CMul (C ar ai) (C br bi)) = (C (- (* ar br) (* ai bi)) (+ (* ar bi) (* ai br)))
(CPol ang)                 = (C (& ang 0.0) (& (/ PI 2.0) ang))

(Com num) = (Com.val num E)
  (Com.val (V x)     path) = (CSca (Com.twi path (/ PI 2.0) (C 1.0 0.0)) x)
  (Com.val (G a0 a1) path) = (CAdd (Com.val a0 (O path)) (Com.val a1 (I path)))
  (Com.twi E        ang num) = num
  (Com.twi (O path) ang num) = (Com.twi path (/ ang 2.0) num)
  (Com.twi (I path) ang num) = (Com.twi path (/ ang 2.0) (CMul num (CPol ang)))

// Tree

(Map f (L x))   = (L (f x))
(Map f (B x y)) = (B (Map f x) (Map f y))

(Flt (L x))   = [x]
(Flt (B x y)) = (List.concat (Flt x) (Flt y))

(Idx E     (L x))   = x
(Idx (O p) (B x y)) = (Idx p x)
(Idx (I p) (B x y)) = (Idx p y)

(Dep (L x))   = 0
(Dep (B x y)) = (+ 1 (Dep x))

(Inv tree) = (Inv.go (Dep tree) E tree)
  (Inv.go 0 path tree) = (L (Idx path tree))
  (Inv.go n path tree) = (B (Inv.go (- n 1) (O path) tree) (Inv.go (- n 1) (I path) tree))

(See x) = (Flt (Inv (Map λx(Com x) x)))

// Main

Main =
  let c0 = 1.0
  let c1 = 2.0
  let c2 = 3.0
  let c3 = 4.0
  let c4 = 5.0
  let c5 = 6.0
  let c6 = 7.0
  let c7 = 8.0
  (See (FFT (B (B (B (L c0) (L c4)) (B (L c2) (L c6))) (B (B (L c1) (L c5)) (B (L c3) (L c7))))))

Improved:

BE     = λx(x)
(B0 f) = λx(f (f x))
(B1 f) = λx(Rot (f (f x)))

(Rot (V z))   = (V (- 0.0 z))
(Rot (G x y)) = (G (Rot y) x)

(Add (V z)   (V w))   = (V (+ z w))
(Add (G x y) (G w z)) = (G (Add x w) (Add y z))

(Sub (V z)   (V w))   = (V (- z w))
(Sub (G x y) (G w z)) = (G (Sub x w) (Sub y z))

(Get (G x y) f) = (f x y)

(Mix ang (B ax ay) (B bx by)) =
  let ptl = (Mix λx(ang (B0 x)) ax bx)
  let ptr = (Mix λx(ang (B1 x)) ay by)
  (B ptl ptr)
(Mix ang (L pt0) (L pt1)) =
  let pt2 = (ang λx(x) pt1)
  let ptl = (L (Add pt0 pt2))
  let ptr = (L (Sub pt0 pt2))
  (B ptl ptr)

(FFT ang (L x)) = (L x)
(FFT ang (B x y)) =
  let pt0 = (FFT λx(ang (B0 x)) x)
  let pt1 = (FFT λx(ang (B0 x)) y)
  (Mix ang pt0 pt1)

// Conversions
// -----------

// List

(List.concat List.nil         ys) = ys
(List.concat (List.cons x xs) ys) = (List.cons x (List.concat xs ys))

// Complex

PI = 3.141592653589793

(CSca (C ar ai) s)         = (C (* ar s) (* ai s))
(CAdd (C ar ai) (C br bi)) = (C (+ ar br) (+ ai bi))
(CMul (C ar ai) (C br bi)) = (C (- (* ar br) (* ai bi)) (+ (* ar bi) (* ai br)))
(CPol ang)                 = (C (& ang 0.0) (& (/ PI 2.0) ang))

(Zer 0) = (V 0.0)
(Zer n) = let z = (Zer (- n 1)); (G z z)

(Neo 0 x) = (V x)
(Neo n x) = (G (Neo (- n 1) x) (Zer (- n 1)))

(Com num) = (Com.val num E)
  (Com.val (V x)     path) = (CSca (Com.twi path (/ PI 2.0) (C 1.0 0.0)) x)
  (Com.val (G a0 a1) path) = (CAdd (Com.val a0 (O path)) (Com.val a1 (I path)))
  (Com.twi E        ang num) = num
  (Com.twi (O path) ang num) = (Com.twi path (/ ang 2.0) num)
  (Com.twi (I path) ang num) = (Com.twi path (/ ang 2.0) (CMul num (CPol ang)))

// Tree

(Map f (L x))   = (L (f x))
(Map f (B x y)) = (B (Map f x) (Map f y))

(Flt (L x))   = [x]
(Flt (B x y)) = (List.concat (Flt x) (Flt y))

(Idx E     (L x))   = x
(Idx (O p) (B x y)) = (Idx p x)
(Idx (I p) (B x y)) = (Idx p y)

(Dep (L x))   = 0
(Dep (B x y)) = (+ 1 (Dep x))

(Inv tree) = (Inv.go (Dep tree) E tree)
  (Inv.go 0 path tree) = (L (Idx path tree))
  (Inv.go n path tree) = (B (Inv.go (- n 1) (O path) tree) (Inv.go (- n 1) (I path) tree))

(See x) = (Flt (Inv (Map λx(Com x) x)))

// Main

(Rol E     x)       = x
(Rol (O p) (G a b)) = (G (Rol p a) (Rol p b))
(Rol (I p) (G a b)) = (Rot (G (Rol p a) (Rol p b)))

Main =
  let c0 = (Neo 3 1.0)
  let c1 = (Neo 3 2.0)
  let c2 = (Neo 3 3.0)
  let c3 = (Neo 3 4.0)
  let c4 = (Neo 3 5.0)
  let c5 = (Neo 3 6.0)
  let c6 = (Neo 3 7.0)
  let c7 = (Neo 3 8.0)
  (See (FFT λx(B0 x) (B (B (B (L c0) (L c4)) (B (L c2) (L c6))) (B (B (L c1) (L c5)) (B (L c3) (L c7))))))

// [(C 36.0 0.0), (C -4.0 -9.656854249492369), (C -4.0 -4.0), (C -4.0 -1.6568542494923832), (C -4.0 0.0), (C -3.999999999999993 1.6568542494923832), (C -3.999999999999993 4.0), (C -3.999999999999986 9.656854249492369)]
// [(C 36.0 0.0), (C -4.0 -9.656854249492369), (C -4.0 -4.0), (C -4.0 -1.6568542494923832), (C -4.0 0.0), (C -3.999999999999993 1.6568542494923832), (C -3.999999999999993 4.0), (C -3.999999999999986 9.656854249492369)]
// [(C 36.0 0.0), (C -4.0 -9.656854249492369), (C -4.0 -4.0), (C -4.0 -1.6568542494923832), (C -4.0 0.0), (C -3.999999999999993 1.6568542494923832), (C -3.999999999999993 4.0), (C -3.999999999999986 9.656854249492369)]