# Too Fast, Too Fourier: FFT Code Golf

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)])

Edit:

• the code should implement the standard forward Fast Fourier Transform, the form of which can be seen in equation (3) of this Wolfram article,

• 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.
• This is underspecified. At the very least you need to define the normalisation factors, and you also ought to be aware that any ambiguity will be wilfully misinterpreted. E.g. is "Implement" satisfied by the answer "FFT (3 chars): it's in the standard library"? Some test cases would be good too. – Peter Taylor Aug 30 '13 at 7:24
• Does it matter about the order of the output elements, i.e. do we need to implement bit reversed unscrambling or can we leave the output in scrambled order ? – Paul R Aug 30 '13 at 8:22
• See the edits to the rules. The output should be a list/array with values ordered according to the indices in the standard DFT expression, referenced above. – jakevdp Aug 30 '13 at 16:05
• VTC as no IO requirements mentioned. – user80551 Apr 15 '14 at 11:14
• Can you post some example inputs and outputs so we can test our implementations? – FUZxxl Mar 8 '15 at 11:57

# Mathematica, 95 bytes

Another implementation of the Cooley–Tukey FFT with help from @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]]

• I think #[[;;;;2]]==#[[1;;N;;2]] and [[2;;;;2]]==[[2;;N;;2]]. – chyanog Sep 8 '13 at 16:07
• 101 characters：With[{L=Length@#},If[L>1,Join[+##,#-#2]&[#0@#[[;;;;2]],#0@#[[2;;;;2]]E^(-2I*Pi(Range[L/2]-1)/L)],#]]& – chyanog Sep 8 '13 at 16:20
• Nice, you can condense another anonymous function inside it without conflicting with the recursive one. Also learned that Part fills in missing indices. We can take it further using Unicode. – miles Sep 8 '13 at 17:30

# Python, 166151 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

• 2 things: it's typically best to use from x import*, and sum(([x for x in y] for y in z),[]) is longer than [x for y in z for x in y]. – boothby Sep 1 '13 at 1:40
• Thanks - that saves 15 characters! 11 more and it's a tweet. – jakevdp Sep 4 '13 at 3:25
• Oh, that's definitely possible. Often when you find one improvement, an old one becomes a stumbling block. – boothby Sep 4 '13 at 6:11

# 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. Try it online! ## 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


### R: 14213399 95 bytes

Thanks to @Giuseppe 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

• Four years later, and we get it down to 101 bytes. Not 100% sure why you used Recall as it's a named function already, but hey, it's easy to golf in hindsight! :) +1, very nice. – Giuseppe Sep 10 '17 at 22:05
• Yes Recall is now unnecessary, indeed. I noticed that a few months ago but was too lazy to change it :) I'll modify it. – plannapus Sep 11 '17 at 9:18
• Very nice! I squeezed another 4 bytes out!, putting this on par with Mathematica. – Giuseppe Sep 11 '17 at 11:22
• Thanks! I thought about putting y up there but didn't notice it could be use for the exp(...) part as well. – plannapus Sep 11 '17 at 11:53

## 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: lambdaize, 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])))]

• I think that this is no longer a fast fourier transform, though... by "killing t" you added in some unnecessary calculations which move it from O[N log(N)] to O[N^2] – jakevdp Apr 2 '14 at 20:38
• It appears that I cannot downvote my own post. You are correct, I exchanged the loop order and killed the performance. I will leave this up for now, in case I find a way to fix it. – boothby Apr 9 '14 at 3:52
• 101 bytes with f=lambda x:x*(len(x)<2)or[u+v/1j**(4*i/len(x))for i,(u,v)in enumerate(zip(f(x[::2])*2,f(x[1::2])*2))] – miles Aug 12 '17 at 14:26
• You can replace for s in(1,-1)for with for s in 1,-1for or even, if order does not matter, for s in-1,1for. – Jonathan Frech Sep 10 '17 at 14:23

# 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.

• You can remove some more whitespace to get a lower amout of characters, for example the step < n can be changed into step<n and the step * 2 can be changed into step*2. – ProgramFOX Apr 15 '14 at 11:32
• all the variables and functions and typedefs should have one-letter names to save a lot of chars – user16402 Apr 15 '14 at 11:43
• You had someone suggest a bunch of improvements for this. Take a look at them here: codegolf.stackexchange.com/review/suggested-edits/17119 – Justin Nov 24 '14 at 4:35
• You can remove all newlines, newlines are useless in C – TuxCrafting Jun 18 '16 at 13:36
• @TùxCräftîñg Not all newlines are useless. They are needed for directives like #include, #define, #if, etc. – Nayuki Aug 12 '17 at 1:50

## Python 3: 140134 113 characters

Short version - short and sweet, fits in a tweet (with thanks to 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
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)]

• Shortened to 113 bytes using e^(-2j * pi * i / n) = (-1)^(2 * i / n) = (1j)^(4 * i / n) – miles Aug 11 '17 at 12:04
• @miles Very impressive observation, thank you! Having repeatedly implemented DFTs for over a decade, I got obsessed with sin/cos/exp and forgot that simple powers of i can be used. I edited my answer to incorporate the new insight and credit you. – Nayuki Aug 12 '17 at 1:46

# Matlab, 12811810710210194 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

• Tip: you can combine the two d= lines into one: m=n/2;d=f(Y(2:2:n)).*exp(-pi*i*(0:m-1)/m).';. Furthermore, consider changing y=f(Y) to Y=f(Y) and remove line 3 (and promise you'll never do that outside of code-golf) – Sanchises Mar 8 '15 at 9:12
• Oh thanks! Does function Y = f(Y) have any disadvanteages other than unreadability? – flawr Mar 8 '15 at 10:13
• Well, MATLAB will never complain about a return value, even if Y is never changed. It is slightly faster though, so I guess it isn't so bad after all for some purposes (i.e. a function that almost never changes the input variable) – Sanchises Mar 8 '15 at 10:48
• Now, for shaving off more: m=n/2 could be removed, and instead m replaced by n/2 and n*2 respectively. And then, I strongly believe, the program is as short as could be in MATLAB. – Sanchises Mar 8 '15 at 21:05
• And then, I strongly believe, the program is as short as could be in MATLAB. – Sanchises Mar 8 '15 at 21:05 Famous last words... – Sanchises Aug 14 '17 at 8:33

# Jelly, 313028 26 bytes, non-competing

LḶ÷$N-*×,N$+ḷF
s2Zß€ç/µ¹Ṗ?


Jelly was created after this challenge so it is non-competing.

This uses the Cooley-Tukey radix-2 recursive algorithm. For an un-golfed version, see my answer in Mathematica.

## 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


# C (gcc), 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!

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));
}


# 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]

• Isn't this the naive DFT? (ie theta(N^2)) – miles Aug 31 '13 at 23:59

# 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.