[Question]
  [
Implement the Fast Fourier Transform in the fewest possible characters.
Rules:
* Shortest solution wins
* It can be assumed that the input is a 1D array whose length is a power of two.
* You may use the algorithm of your choice, but the solution must actually be a Fast Fourier Transform, not just a naive Discrete Fourier Transform (that is, it must have asymptotic computation cost of \$O(N \log N)\$)
* the code should implement the standard forward Fast Fourier Transform, the form of which can be seen in equation (3) of [this Wolfram article](http://mathworld.wolfram.com/DiscreteFourierTransform.html),
$$F\_n = \sum\_{k=0}^{N-1}f\_ke^{-2\pi ink/N}$$
* Using an FFT function from a pre-existing standard library or statistics package is not allowed. The challenge here is to succinctly *implement* the FFT algorithm itself.
  
[Answer]
# Mathematica, 95 bytes
Another implementation of the Cooley–Tukey FFT with help from @[chyaong](https://codegolf.stackexchange.com/users/7254/chyaong).
```
{n=Length@#}~With~If[n>1,Join[+##,#-#2]&[#0@#[[;;;;2]],#0@#[[2;;;;2]]I^Array[-4#/n&,n/2,0]],#]&
```
### Ungolfed
```
FFT[x_] := With[{N = Length[x]},
  If[N > 1,
    With[{a = FFT[ x[[1 ;; N ;; 2]] ], 
          b = FFT[ x[[2 ;; N ;; 2]] ] * Table[E^(-2*I*Pi*k/N), {k, 0, N/2 - 1}]},
      Join[a + b, a - b]],
    x]]
```
[Answer]
# J, 37 bytes
```
_2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#
```
An improvement after a few years. Still uses the Cooley-Tukey FFT algorithm.
Saved 4 bytes using *e**πi* = -1, thanks to @[Leaky Nun](https://codegolf.stackexchange.com/users/48934/leaky-nun).
[Try it online!](https://tio.run/##HYu7DoJAFAV7vuJE5BXxuu9lN2qwsrewEimMRG0saCjEX1/BnEwmU5xX6LHzYJgIrUhzluercl00ScuvT6rjJC7Spd/UH99cii/fxqGIFoSsm28ZSowefTTMJTASDqfjObrfHm90KS0o3WMIgeO/yQISKmgSyggDSY5ZNdsKKTU0WaWEgyLrKomKpOKWT6mMrTQYWcO4dD8)
## Usage
```
   f =: _2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#
   f 1 1 1 1
4 0 0 0
   f 1 2 3 4
10 _2j2 _2 _2j_2
   f 5.24626 3.90746 3.72335 5.74429 4.7983 8.34171 4.46785 0.760139
36.9894 _6.21186j0.355661 1.85336j_5.74474 7.10778j_1.13334 _0.517839 7.10778j1.13334 1.85336j5.74474 _6.21186j_0.355661
```
## Explanation
```
_2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#  Input: array A
                                    #  Length
                                  1<   Greater than one?
_2&(                            )~     Execute this if true, else return A
_2                            ]\         Get non-overlapping sublists of size 2
    0                       |:           Move axis 0 to the end, equivalent to transpose
                          /@             Reduce [even-indexed, odd-indexed]
                       &$:               Call recursively on each 
                   #                     Get the length of the odd list
                i.@                      Range from 0 to that length exclusive
                    %#                   Divide each by the odd length
             _1^                         Compute (-1)^x for each x
           ]                             Get the odd list
            %                            Divide each in that by the previous
       +                                 Add the even values and modified odd values
         -                               Subtract the even values and modified odd values
        ,                                Join the two lists and return
```
[Answer]
# Python, 166 151 150 characters
This uses the radix-2 Cooley-Tukey FFT algorithm
```
from math import*
def F(x):N=len(x);t=N<2or(F(x[::2]),F(x[1::2]));return N<2and x or[
a+s*b/e**(2j*pi*n/N)for s in[1,-1]for(n,a,b)in zip(range(N),*t)]
```
Testing the result
```
>>> import numpy as np
>>> x = np.random.random(512)
>>> np.allclose(F(x), np.fft.fft(x))
True
```
[Answer]
## Python 3: 140 134 113 characters
Short version - short and sweet, fits in a tweet (with thanks to [miles](https://codegolf.stackexchange.com/users/6710/miles)):
```
from math import*
def f(v):
 n=len(v)
 if n<2:return v
 a,b=f(v[::2])*2,f(v[1::2])*2;return[a[i]+b[i]/1j**(i*4/n)for i in range(n)]
```
(In Python 2, `/` is truncating division when both sides are integers. So we replace `(i*4/n)` by `(i*4.0/n)`, which bumps the length to 115 chars.)
Long version - more clarity into the internals of the classic Cooley-Tukey FFT:
```
import cmath
def transform_radix2(vector):
    n = len(vector)
    if n <= 1:  # Base case
        return vector
    elif n % 2 != 0:
        raise ValueError("Length is not a power of 2")
    else:
        k = n // 2
        even = transform_radix2(vector[0 : : 2])
        odd  = transform_radix2(vector[1 : : 2])
        return [even[i % k] + odd[i % k] * cmath.exp(i * -2j * cmath.pi / n) for i in range(n)]
```
[Answer]
### R: 142 ~~133~~ ~~99~~ 95 bytes
Thanks to [@Giuseppe](https://codegolf.stackexchange.com/users/67312) for helping me shaving down ~~32~~ 36 bytes!
```
f=function(x,n=sum(x|1),y=1:(n/2)*2)`if`(n>1,f(x[-y])+c(b<-f(x[y]),-b)*exp(-2i*(y/2-1)*pi/n),x)
```
An additional trick here is to use the main function default arguments to instantiate some variables.  
Usage is still the same:
```
x = c(1,1,1,1)
f(x)
[1] 4+0i 0+0i 0+0i 0+0i
```
**4-year old version at 133 bytes:**
```
f=function(x){n=length(x);if(n>1){a=Recall(x[seq(1,n,2)]);b=Recall(x[seq(2,n,2)]);t=exp(-2i*(1:(n/2)-1)*pi/n);c(a+b*t,a-b*t)}else{x}}
```
With indentations:
```
f=function(x){
    n=length(x)
    if(n>1){
        a=Recall(x[seq(1,n,2)])
        b=Recall(x[seq(2,n,2)])
        t=exp(-2i*(1:(n/2)-1)*pi/n)
        c(a+b*t,a-b*t)
        }else{x}
    }
```
It uses also Cooley-Tukey algorithm. The only tricks here are the use of function `Recall` that allows recursivity and the use of R vectorization that shorten greatly the actual computation.
Usage:
```
x = c(1,1,1,1)
f(x)
[1] 4+0i 0+0i 0+0i 0+0i
```
[Answer]
## Python, 134
This borrows heavily from jakevdp's solution, so I've set this one to a community wiki.
```
from math import*
F=lambda x:x*(len(x)<2)or[a+s*b/e**(2j*pi*n/len(x))for s in(1,-1)for n,(a,b)in
enumerate(zip(F(x[::2]),F(x[1::2])))]
```
Changes:
-12 chars: kill `t`.
```
def F(x):N=len(x);t=N<2or(F(x[::2]),F(x[1::2]));return ... in zip(range(N),*t)]
def F(x):N=len(x);return ... in zip(range(N),F(x[::2]),F(x[1::2]))]
```
-1 char: exponent trick, `x*y**-z == x/y**z` (this could help some others)
```
...[a+s*b*e**(-2j*pi*n/N)...
...[a+s*b/e**(2j*pi*n/N)...
```
-2 char: replace `and` with `*`
```
...return N<2and x or[
...return x*(N<2)or[
```
+1 char: `lambda`ize, killing `N`
```
def F(x):N=len(x);return x*(N<2)or[a+s*b/e**(2j*pi*n/N) ... zip(range(N) ...
F=lambda x:x*(len(x)<2)or[a+s*b/e**(2j*pi*n/len(x)) ... zip(range(len(x)) ...
```
-2 char: use `enumerate` instead of `zip(range(len(`
```
...for(n,a,b)in zip(range(len(x)),F(x[::2]),F(x[1::2]))]
...for n,(a,b)in enumerate(zip(F(x[::2]),F(x[1::2])))]
```
[Answer]
# C, 259
```
typedef double complex cplx;
void fft(cplx buf[],cplx out[],int n,int step){
if(step < n){
fft(out, buf,n, step * 2);
fft(out+step,buf+step,n,step*2);
for(int i=0;i<n;i+=2*step){
cplx t=cexp(-I*M_PI*i/n)*out[i+step];
buf[i/2]=out[i]+t;
buf[(i+n)/2]=out[i]-t;
}}}
```
The problem is, such implementations are useless, and straightforward algorithm is MUCH faster.
[Answer]
# Matlab, ~~128~~ ~~118~~ ~~107~~ ~~102~~ ~~101~~ ~~94~~ 93 bytes
EDIT6: thanks @algmyr for another byte!
```
function Y=f(Y);
n=numel(Y);
k=2:2:n;
if k;
   c=f(Y(k-1));
   d=f(Y(k)).*i.^(2*(2-k)/n);
   Y=[c+d;c-d];
end
```
EDIT5: Still getting shorter:) thanks to @sanchises
```
function Y=f(Y)
n=numel(Y);
k=2:2:n;
if k;
   c=f(Y(k-1));
   d=f(Y(k)).*(-1).^((2-k)/n);
   Y=[c+d;c-d];
end
```
EDIT4: Yay, -1 character more (could aslo have done without the `k`):
```
function Y=f(Y)
n=numel(Y);
if n>1;
   k=2:2:n;
   c=f(Y(k-1));
   d=f(Y(k)).*(-1).^((k/2-1)*2/n)';
   Y=[c+d;c-d];
end
```
EDIT2/3: Thanks for @sanchises for further improvements!
```
function Y=f(Y)
n=numel(Y);  
if n>1;
   c=f(Y(1:2:n));
   d=f(Y(2:2:n)).*(-1).^(-(0:n/2-1)*2/n).';
   Y=[c+d;c-d]; 
end
```
EDIT: Could make some improvements, and noticed that the scaling constant is not required.
This is the expanded version, character count is valid if you remove the newlines/spaces. (Works only for column vectors.)
```
function y=f(Y)
n=numel(Y);  
y=Y;
if n>1;
   c=f(Y(1:2:n));
   d=f(Y(2:2:n));
   n=n/2;
   d=d.*exp(-pi*i*(0:n-1)/n).';
   y=[c+d;c-d]; 
end
```
[Answer]
# [C (gcc)](https://gcc.gnu.org/), ~~188 186 184~~ 183 bytes
```
#define d(a,b,c)f(a,b,c,1,0)
f(a,b,c,n,k)_Complex*a,*b;{_Complex z[c];*b=*a;if(n<c)for(f(a,z,c,n*2),f(a+n,z+n,c,n*2);k<c;k+=n*2)b[k+c>>1]=z[k]*2-(b[k/2]=z[k]+z[k+n]/cpow(1i,2.*k/c));}
```
[Try it online!](https://tio.run/##rVDBasMwDL33K0RHwY7dNUnDKLjuZZ@RmWG7cTFJnNB0bKT015fZzVo22NihO0jPT5b0JOn5TuthuNsWxroCtkhSRTU2I9KExnhyIY6W@PmxqduqeIskjRQ7Xij0uRYsUjySzBrk1r5Fs0ehsg@VUYqpJ8TR3tsYYOVas5Lw8FZ5SfRmkwje56WI0jnykUU6UuIdcWKh2@YVJZam91G50Biz03DVl3kSCza5cjVy6w5QS@sQPk7AJ8WC@8zg0uCWgifMf4RZIeRa4BAzD2vIPBACGEIlQLv3/wbBdGbIzFiYUtD7QlZI5lZgT2wtdyMBHHqevLUvhw5Np2e@RSApKAoZvlVSfZVUP0k@uVH0c1m/Z8zO2yfsvHrqIRN86eFB8Cyk1kXdFQcUJowpdLYvGgPq11lX/3ce@H6f1c2af98nXOc0vGtTyV03zKv6Aw "C (gcc) – Try It Online")
Slightly golfed less
```
#define d(a,b,c)f(a,b,c,1,0)
f(a,b,c,n,k)_Complex*a,*b;{
  _Complex z[c];
  *b=*a;
  if(n<c)
    for(f(a,z,c,n*2),f(a+n,z+n,c,n*2);k<c;k+=n*2)
      b[k+c>>1]=z[k]*2-(b[k/2]=z[k]+z[k+n]/cpow(1i,2.*k/c));
}
```
[Answer]
# Jelly, ~~31~~ ~~30~~ ~~28~~ 26 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
```
LḶ÷$N-*×,N$+ḷF
s2Zß€ç/µ¹Ṗ?
```
This uses the Cooley-Tukey radix-2 recursive algorithm. For an un-golfed version, see my [answer](https://codegolf.stackexchange.com/a/12429/6710) in Mathematica.
[Try it online](http://jelly.tryitonline.net/#code=TOG4tsO3JE4tKsOXLE4kK-G4t0YKczJaw5_igqzDpy_CtcK54bmWPw&input=&args=WzEsIDIsIDMsIDRd) or [Verify multiple test cases](http://jelly.tryitonline.net/#code=TOG4tsO3JE4tKsOXLE4kK-G4t0YKczJaw5_igqzDpy_CtcK54bmWPwrDh8WS4bmYJOKCrGrigbc&input=&args=W1sxLCAxLCAxLCAxXSwKIFsxLCAyLCAzLCA0XSwKIFs1LjI0NjI2LCAzLjkwNzQ2LCAzLjcyMzM1LCA1Ljc0NDI5LCA0Ljc5ODMsIDguMzQxNzEsIDQuNDY3ODUsIDAuNzYwMTM5XV0).
## Explanation
```
LḶ÷$N-*×,N$+ḷF  Helper link. Input: lists A and B
L               Get the length of A
   $            Operate on that length
 Ḷ                Make a range [0, 1, ..., length-1]
  ÷               Divide each by length
    N           Negate each
     -          The constant -1
      *         Compute -1^(x) for each x in that range
       ×        Multiply elementwise between that range and B, call it B'  
          $     Operate on that B'
         N        Negate each
        ,         Make a list [B', -B']
            ḷ   Get A
           +    Add vectorized, [B', -B'] + A = [A+B', A-B']
             F  Flatten that and return
s2Zß€ç/µ¹Ṗ?  Main link. Input: list X
         Ṗ   Curtail - Make a copy of X with the last value removed
          ?  If that list is truthy (empty lists are falsey)
       µ       Parse to the left as a monad
s2             Split X into sublists of length 2
  Z            Transpose them to get [even-index, odd-index]
   ß€          Call the main link recursively on each sublist
     ç/        Call the helper link as a dyad on the sublists and return
             Else
        ¹      Identity function on X and return
```
[Answer]
# Pari/GP, 76 characters
```
X(v)=my(t=-2*Pi*I/#v,s);vector(#v,k,s=t*(k-1);sum(n=0,#v-1,v[n+1]*exp(s*n)))
```
Usage
```
X([1,1,1,1])
%2 = [4.000000000000000000000000000, 0.E-27 + 0.E-28*I, 0.E-28 + 0.E-27*I, 0.E-27 + 0.E-28*I]
```
[Answer]
# [Octave](https://www.gnu.org/software/octave/), ~~109 103 101~~ 100 bytes
```
f(f=@(f)@(x,n=rows(x)){@(d=f(f)(x(k=2:2:n)).*i.^((k*2-4)/n)')[d+(c=f(f)(x(k-1)));c-d],x}{1+(n<2)}())
```
[Try it online!](https://tio.run/##TY3BCsIwGIPvPkVv@7N1da3zMi30PUYF2SzIoAMVLYw9ey0diAS@QBLIPLyu71t0WggRHTltyMFQ4F4/5s@TArAYGnWqQIEmrTrVeUCUd3EhmkpVt9h7FOjHiobfrpYATkM9Wh7WRVbkzworAemkt9glHjeTrPmXLbaUs4YzlXnIbHku4xc "Octave – Try It Online")
Ooooo do my eyes bleed from this ~~recursive~~ accursed lambda. Large parts of this were lifted from @flawr's answer.
```
f(                                          % lambda function
  f=@(f)                                    % defined in its own argument list, 
                                            % accepts itself as parameter (for recursion)
    @(x,n=rows(x)){                         % calls another lambda,
                                            % 2nd parameter is just to define a variable
      @(d=f(f)(x(k=2:2:n)).*i.^((k*2-4)/n)')% 1/4 of FFT (argument just defines a variable)
        [d+(c=f(f)(x(k-1)));                % 2/4 of FFT
         c-d                                % 4/4 of FFT
        ],                                  % This is in a @()[] to inhibit evaluation
                                            % unless actually called
      x                                     % FFT of length 1
    }{1+(n<2)}                              % if len(x)==1, return x
                                            % else return [d+c;c-d]
  ()                                        % this is probably important too
)
```
[Answer]
# APL(NARS), 58 chars, 116 bytes
```
{1≥k←≢⍵:⍵⋄(∇⍵[y∼⍨⍳k])(+,-)(∇⍵[y←2×⍳t])×0J1*t÷⍨2-2×⍳t←⌊k÷2}
```
test
```
  f←{1≥k←≢⍵:⍵⋄(∇⍵[y∼⍨⍳k])(+,-)(∇⍵[y←2×⍳t])×0J1*t÷⍨2-2×⍳t←⌊k÷2}
  f 1 1 1 1
4J0 0J0 0J0 0J0 
  f 1 2 3 4
10J0 ¯2J2 ¯2J0 ¯2J¯2 
  f 1J1 2 ¯2J1  9
10J2 3J7 ¯12J2 3J¯7 
  f 5.24626,3.90746,3.72335,5.74429,4.7983,8.34171,4.46785,0.760139
36.989359J0 ¯6.211855215J0.3556612739 1.85336J¯5.744741 7.107775215J¯1.133338726 ¯0.517839J0 
  7.107775215J1.133338726 1.85336J5.744741 ¯6.211855215J¯0.3556612739 
```
[Answer]
# Axiom,  ~~259~~, ~~193~~, ~~181~~, 179 bytes
```
L(g,n,f)==>[g for i in 1..n|f]
h(a)==(n:=#a;n=1=>a;c:=h(L(a.i,n,odd? i));d:=h(L(a.i,n,even? i));n:=n/2;t:=1>0;v:=L(d.i*%i^(-2*(i-1)/n),n,t);append(L(c.i+v.i,n,t),L(c.i-v.i,n,t)))
```
Even if h(a) could pass all the test and would be ok as entry for this 'competition'
one has to call h() or hlp() thru fft() below, for *checking arguments*.
I don't know if this software can be ok because i only had seen what other wrote, and search
the way it could run in Axiom for return some possible right result. Below ungolfed code with few comments:
```
-- L(g,n,f)==>[g for i in 1..n|f]
-- this macro L, build one List from other list, where in g, there is the generic element of index i
-- (as a.i, or a.i*b.i or a.i*4), n build 1..n that is the range of i, f is the condition 
-- for insert the element in the list result.
hlp(a)==
    n:=#a;n=1=>a
    -- L(a.i,n,odd? i)  it means build a list getting "even indices i of a.i as starting from index 0" [so even is odd and odd is even]
    -- L(a.i,n,even? i) it means build a list getting "odd  indices i of a.i as starting from index 0"
    c:=hlp(L(a.i,n,odd? i));d:=hlp(L(a.i,n,even? i))
    n:=n/2;t:=1>0
    v:=L(d.i*%i^(-2*(i-1)/n),n,t)
    append(L(c.i+v.i,n,t),L(c.i-v.i,n,t))
-- Return Fast Fourier transform of list a, in the case #a=2^n
fft(a)==(n:=#a;n=0 or gcd(n,2^30)~=n=>[];hlp(a))
(5) -> h([1,1,1,1])
   (5)  [4,0,0,0]
                                    Type: List Expression Complex Integer
(6) -> h([1,2,3,4])
   (6)  [10,- 2 + 2%i,- 2,- 2 - 2%i]
                                    Type: List Expression Complex Integer
(7) -> h([5.24626,3.90746,3.72335,5.74429,4.7983,8.34171,4.46785,0.760139])
   (7)
   [36.989359, - 6.2118552150 341603904 + 0.3556612739 187363298 %i,
    1.85336 - 5.744741 %i, 7.1077752150 341603904 - 1.1333387260 812636702 %i,
    - 0.517839, 7.1077752150 341603904 + 1.1333387260 812636702 %i,
    1.85336 + 5.744741 %i,
    - 6.2118552150 341603904 - 0.3556612739 187363298 %i]
                                      Type: List Expression Complex Float
(8) -> h([%i+1,2,%i-2,9])
   (8)  [10 + 2%i,3 + 7%i,- 12 + 2%i,3 - 7%i]
                                    Type: List Expression Complex Integer
```
in the few i had seen h() or fft() would return exact solution, but if the simplification
is not good as in:
```
(13) -> h([1,2,3,4,5,6,7,8])
   (13)
                    +--+                                   +--+
        (- 4 + 4%i)\|%i  - 4 + 4%i             (- 4 - 4%i)\|%i  - 4 + 4%i
   [36, --------------------------, - 4 + 4%i, --------------------------, - 4,
                    +--+                                   +--+
                   \|%i                                   \|%i
            +--+                                   +--+
    (- 4 + 4%i)\|%i  + 4 - 4%i             (- 4 - 4%i)\|%i  + 4 - 4%i
    --------------------------, - 4 - 4%i, --------------------------]
                +--+                                   +--+
               \|%i                                   \|%i
                                    Type: List Expression Complex Integer
```
than it is enought change the type of only one element of list, as in below writing 8. (Float)
for find the approximate solution:
```
(14) -> h([1,2,3,4,5,6,7,8.])
   (14)
   [36.0, - 4.0000000000 000000001 + 9.6568542494 923801953 %i, - 4.0 + 4.0 %i,
    - 4.0 + 1.6568542494 92380195 %i, - 4.0, - 4.0 - 1.6568542494 92380195 %i,
    - 4.0 - 4.0 %i, - 4.0 - 9.6568542494 923801953 %i]
                                      Type: List Expression Complex Float
```
I wrote it, seen all other answers because in the link, the page it was too much difficult so I don't know if this code can be right. I'm not one fft expert so all this can (it is probable) be wrong.
]