%%HP: T(3)F(,);
@ by Jurjen NE Bos
DIR
  CST { { "FFT"
    \<< 1 FFT
    \>> } { "-FFT"
    \<< -1 FFT
    \>> } }
  FFT   @ FFT routine.
        @ Input:        2: Vector 1: 1 for FFT, -1 for inverse
        @ Ouput:        1: transformed vector
    \<< OVER SIZE 1 GET
      IF OVER 0 <       @ Divide by length if inverse FFT
      THEN ROT OVER / ROT ROT
      END (0;6,28318530718) ROT * OVER / SWAP DUP
      WHILE DUP 2 MOD NOT       @ Divide out factors of 2
      REPEAT 2 /
      END SWAP OVER / \-> \Gr r p       @ p: twopower, r: odd rest
                                        @ rho: root of unity
      \<< { r p } RDM
        IF r 1 \=/      @ Compute regular r-point FFT p times simultaneously
        THEN A\->L \Gr p * EXP 1 \-> m \Grp \Grn
          \<< 1 r
            START m 1 GET 1 2 r
              FOR k \Grn * m k GET OVER * ROT + SWAP
              NEXT DROP OBJ\-> DROP \Grp '\Grn' STO*
            NEXT
          \>> { r p } \->ARRY
        END TRN CONJ    @ TRN does transpose & conjugate, so conjugate back
        WHILE p 1 \=/   @ combine r-point FFTs to 2r-point FFTs
        REPEAT OBJ\-> DROP 'p' 2 STO/ { p r } \->ARRY TRN A\->L
        @ from here on, we work with conjugated numbers!
          \Gr NEG p * EXP 1 \-> m \Grp \Grn
          \<< { p r } \->ARRY TRN 1 r
            FOR k m k GET \Grn * OBJ\-> DROP \Grp '\Grn' STO*
            NEXT
          \>> { r p } \->ARRY + LASTARG - \-> m
          \<< OBJ\-> DROP m
          \>> OBJ\-> DROP 'r' 2 STO* { r p } \->ARRY TRN
          @ numbers are normal again
        END { r } RDM
      \>>
    \>>
  A\->L @ convert matrix to list of rows.  This doesn't use sigma+, because
        @ this trick doesn't work with complex numbers :-(
    \<< OBJ\-> OBJ\-> DROP { } \-> r c u
      \<< 1 r
        START c \->ARRY 'u' STO+
        NEXT u
      \>>
    \>>
  DFT   @ The regular Fourier Transform.  Extra short version.
    \<< SWAP DUP SIZE 1 GET (0;6,28318530718) OVER / 4 PICK * \-> d x s q
      \<< 0 s 1 -
        FOR k '\GS(n=0;s-1;x(n+1)*EXP(q*k*n))' \->NUM
        NEXT s \->ARRY
        IF d -1 ==
        THEN s /
        END
      \>>
    \>>
END