12 ms·

Fast Fourier transform in x86 assembly

158 points by gozzoo 10y ago

jlebar 10y ago
> gcc -O1 -o fft-test fft-test.c fft-portable.c -lm -O1 instructs the compiler to leave some optimizations on the table. If you want the fastest code the compiler can generate, you want -O2 or -O3. This is true for clang and gcc.
- user5994461 10y ago
  And I believe he should also give flags to target the specific hardware and instruction set that he has. In extreme cases, GCC/Intel/MSVC are able to optimize pure C code with auto-vectorization and SSE/AVX instructions, possibly better than hand crafted ASM... but the compiler will never do that if not given the flags to allow it.
  - astrange 10y ago
    Autovectorization is okay but not that good. Its advantage is just that it can transform any part of the program, when your limited dev time won't touch it all. But it can't focus on the hot loops or know about cache efficiency etc. even if you use PGO. And it's sensitive and unstable wrt code changes which is never good.
    - aarongolliver 10y ago
      Yeah, and when autovectorization fails it can be really hard to tell why. The diagnostic messages aren't the best (or, at least, they weren't when I was doing this ~4 years ago) EG: ICC has a AVX version of atan2(), which you can use through intrinsics[0]. But it refused to vectorize even the simplest loop with the function call in it. Even with every possible fast-math flag or every #pragma(please vectorize), it wouldn't do it. I assume due to reduced precision? Nothing I could find online had the answer (again, 4 years ago) [0] https://software.intel.com/en-us/node/524364 https://software.intel.com/en-us/node/524364
- benchaney 10y ago
  To be honest, this oversight makes the rest of his analysis somewhat suspect.
  - dagss 10y ago
    I would assume it meant that it made things slower, after all this was a hand-tuned numerical core. Optimizations is for what a human did not already optimize.
    - benchaney 10y ago
      If by this you mean higher levels of optimization made things slower, than I must disagree. He says "Compiling C code with optimizations enabled is obviously a good idea" he would have mentioned it here if -O2 or -O3 were slower. Also, the optimizer doesn't optimize hand written assembly. Only c code. What I actually think happened is he didn't understand that there were higher optimization levels that -O1, which gives his code an advantage when compared to the c code implementation.
      - dagss 10y ago
        He compared with baseline C code with -O1? That is really bad.
        benchaney 10y ago
        Well, it is hard to actually tell but that is what it looks like.
llukas 10y ago
Any comparison to fftw or mkl?
- nitrogen 10y ago
  Also kissfft numbers would be interesting for comparison. Even though it's slower, it has a very easy license to work with.
- etrain 10y ago
  I would also be interested in seeing this.
kken 10y ago
"The hand-coded x86-64 AVX code is from 1.0× to 2.4× faster than the best C code (which is the AVX model code). This speedup was not as impressive as hoped." Curious. So assembly optimization does not really help that much with modern CPUs?
- iammyIP 10y ago
  the usual speedup with sse/avx is in about that 2.5 range from what i've observed. If it's lower, then there are likely other effects that stall the execution, like cache misses and other data logistic ineffieciency. The fft shown here uses a precalced lookup table for sin & cos. For a fft size of 2048 this takes 2 * 8 * 2048 bytes, around 32 kb, which is the already size of the usual L1 cache. The bigger the fft size the bigger the lookup table and the less all the asm will be effective. The asm benefits more of a pure algorithmic approach, might try some sin cos approximation instead of a lookup table, then the speed gain might be more consistent.
- user5994461 10y ago
  Whatever can be done with assembly can be done the same with C (or C++). So strictly speaking, assembly optimization can never help you [compared to not using assembly].
  - benchaney 10y ago
    This isn't true at all, in fact, the opposite is true. Compilers are often better than naive micro optimizations, but they do not always have every optimization programmed for every architecture. An expert can absolutely beat the compiler in some cases.
    - Veedrac 10y ago
      Although I agree the compiler rarely produces great code automatically, it's normally possible to handhold the compiler into producing code roughly equivalent to hand-written ASM. Doing so is rarely easier, but it is at least more portable!
mathetic 10y ago
The page doesn't mention compiling code with -O2 or -O3. Does anyone know the reason for this? In particular -O3 is meant to enable vectorisation optimisations, which seems like it might subsume the author's hand optimisations. Also I believe Intel C Compiler is better at optimising numerical code, so results with that would also be illuminating.
- raverbashing 10y ago
  You can enable vectorization manually without going for -O2 or -O3
  - loeg 10y ago
    It isn't clear the author did so, though.
    - mathetic 10y ago
      Actually, it is clear the author didn't because he gives the shell commands that are used to compile.
andrepd 10y ago
A few questions: Is this just vectorization? Because FFT libraries often implement this using either asm or intrinsics. Why compile with -O1? Why not let the gcc optimize with -O2 or -O3? In recent versions I believe this enables vectorization attempts which could produce similar speedups to your handwritten asm.
mpreda 10y ago
The GIMPS (Mersenne numbers search) project has a really extremely fast AVX2 assembler implementation of double FFT http://www.mersenne.org/download/ http://www.mersenne.org/download/ .
- adito 10y ago
  So, what it's licensed under? Can I use it in other project? From quick glance, it seems that they have their own license.
  - mwfunk 10y ago
    The terms and conditions are right here, you can evaluate it for yourself (you are the only person who knows if it suits your needs or not): http://www.mersenne.org/legal/#TCU http://www.mersenne.org/legal/#TCU
rjtobin 10y ago
Here are some performances with higher optimization parameters and with the Intel C Compiler. All run with a Haswell Xeon (E5-2698 v3). Using icc, for the larger problem sizes you get performance speedup that very similar to the speedup from the hard-coded asm. With same arguments as blog post: gcc -O1 -o fft-test fft-test.c fft-portable.c -lm: ... 262144 min=15721397 mean=15748901 sd=0.12% 524288 min=34474976 mean=34546878 sd=0.18% 1048576 min=99814826 mean=100391926 sd=0.26% ... Adding -march=native and -O2: ... 262144 min=14668178 mean=14701671 sd=0.14% 524288 min=31297081 mean=31395468 sd=0.14% 1048576 min=94027867 mean=94196746 sd=0.13% ... With Intel C Compiler (-O2 -xHost): ... 262144 min=14258167 mean=14297606 sd=0.16% 524288 min=31126755 mean=31212153 sd=0.23% 1048576 min=92893600 mean=93491419 sd=0.46% ... Can provide more numbers if people want them. EDIT: looks like I was using the less optimized C code, will re-run with the better performing C code when I get the chance.
- inlineint 10y ago
  Did you try GCC with -O3 and -ffast-math? It seems like ICC has its analogue to ffast-math enabled by default: http://scicomp.stackexchange.com/questions/20815/performance-of-icc-main-cpp-g-ffast-math-main-cpp http://scicomp.stackexchange.com/questions/20815/performance...
  - conistonwater 10y ago
    Good point. GCC has a whole bunch of options here: https://gcc.gnu.org/wiki/FloatingPointMath https://gcc.gnu.org/wiki/FloatingPointMath, and icc's default setting is some, possibly unknown, combination of them.
- simcop2387 10y ago
  I just started a run with the better C code from a sibling commentor. I've also got it doing -ffast-math, and each of the -O1,2,3. I'll post the raw data up to https://gist.github.com/simcop2387/32c8fb863403ee714227afe0ee3dd637 https://gist.github.com/simcop2387/32c8fb863403ee714227afe0e...
tsomctl 10y ago
How does it compare to the two giants, FFTW and Intel IPP? (http://www.fftw.org/speed/CoreDuo-3.0GHz-icc64/ http://www.fftw.org/speed/CoreDuo-3.0GHz-icc64/). In particular, I don't see any mention of cache awareness. I don't know much about FFTW, but all the parts of Intel IPP that I've used are ridiculously fast, and I'm sure Intel is using all the tricks they can to speed up the cache (ie proper layout, cache prefetch). For example, Intel's fft forward transform on complex floats is about as fast as memcpy with a cold cache.
- dsharlet 10y ago
  It's not going to compare well. The assembly is a pretty simple port of a basic C implementation of radix 2 Cooley-Tukey FFTs. The author gets ~2x speedup from writing this in assembly. In contrast, highly optimized FFT implementations that use higher radix transforms and other FFT algorithm optimizations, along with SIMD and good C (not necessarily assembly) will get speedups of 10-40x over basic C implementations with a radix 2 implementation: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.153.6089&rep=rep1&type=pdf http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.153... (Implementing FFTs in practice by Frigo and Johnson, of FFTW) FWIW, it is possible to significantly beat FFTW, if you highly optimize for particular cases: https://github.com/halide/Halide/pull/977 https://github.com/halide/Halide/pull/977. I don't have data handy, but I've seen similar results comparing against IPP.
- gsteinb88 10y ago
  FFTW is one of the seminal examples of the importance of cache awareness, if I remember my history correctly. As I understand it, the reason the fftw authors (Matteo Frigo and Steven G. Johnson) won the Wilkinson prize was essentially because fftw exemplified the importance of cache-aware algorithms, as opposed to blindly optimizing number of operations. I don't know, but I would guess that intel's implementation post-dates fftw?
- alfalfasprout 10y ago
  Probably a massive difference. The generic X86 IPP implementation of 1 dimensional FFT is around 8x faster than the naive FFT variants I've tried. Intel also claims an effective 5.5x improvement over that when targeting AVX2 (so Haswell and later). Essentially, a 40x is not unusual. When it comes to insanely optimized code, most of Intel's libraries (MKL, VSL, IPP, DAAL) are hard to beat.
ziikutv 10y ago
Great work, just finished a course in DSP too!
throwaway_26758 10y ago
I always cringe when developers compare compiler-optimized C to 'hand rolled' assembly. In the dawning age of SAT-based program synthesis and vectorization-aware superoptimizers, it's extremely anachronistic. It only shows a lack of awareness that in reality your job as a programmer is not to give instructions to a machine, but to provide enough help that a machine can create its own instructions better than you can. For example, let's take one critical function here: bit reversal. The C code provided (renamed variables so it's easy to consistently compare) has an impressive speedup just going from -O1 to higher optimization levels. And then, on top of that speedup with a decent optimizer, we can again halve the amount of time spent in this function by providing better intentions to the compiler which allow it to perform automatic vectorization. Note this is all in standard C, and makes the code simpler to understand (and maintain) by humans as well. uint64_t reverse_slow(uint64_t n, uint64_t k) { uint64_t r, i; for (r = 0, i = 0; i < k; ++i, n >>= 1) r = (r << 1) | (n & 1); return r; } uint64_t reverse_fast(const uint64_t n, const uint64_t k) { uint64_t r, i; for (r = 0, i = 0; i < k; ++i) r |= ((n >> i) & 1) << (k - i - 1); return r; } Lesson: do your job as a human and provide guidance, let the machines do their job at taking your advice and running wild with creativity. It helps you both become more effective.
- conistonwater 10y ago
  While what you write has a point, I think it's not quite right. I remember stumbling on this HN conversation recently: https://news.ycombinator.com/item?id=9202858 https://news.ycombinator.com/item?id=9202858 — the two realistic examples mentioned there are video codecs and language interpreters. > In the dawning age of SAT-based program synthesis and vectorization-aware superoptimizers, it's extremely anachronistic. I find this particular argument unconvincing, because those aren't real yet. (Consider also that if it was easy and feasible, they probably would have been real already.) In particular, with respect to automatic vectorization (not even superoptimizers!), there is one curious thing that should be noted: clang, with its automatic vectorizer, gives you special options to check if the vectorizer is working for your code (http://llvm.org/docs/Vectorizers.html#diagnostics http://llvm.org/docs/Vectorizers.html#diagnostics). This is necessary: from the point of view of the programmer, performance-sensitive code that is not auto-vectorized is miscompiled because vectorizing it and the consequent performance gain is part of its correctness. One group of examples of this is all the linear algebra libraries. > do your job as a human and provide guidance This advice is only good in the general setting, where it is indeed good advice. The people whose job it is to write all that performance-sensitive code absolutely must get it right, and saying "a compiler will eventually do it for you" is simply not good enough. FWIW, I think you're not addressing your argument at the right crowd.
edblarney 10y ago
Performance? Man, I wish someone would release a decent tutorial on how to actually use the FFT in real life. There are some decent explanations of FT, but often they fail at the very end of FFT: what the output means, explicitly. Particularly, this concern of 'the same freq. appearing in many buckets' due to issues of the phase of the wave? So much theory, so much Matlab ... so little actual applied Engineering. Making 'stuff work' in the real world is basket of caveats ... looking forward to that article: 'FFT in the wild' :)
- ww520 10y ago
  I was baffled by FFT, too, when I first encountered in math class. It seemed useless. It's not until I started working with sound that it really dawned on me how it's being used. I think the best illustrated usage of DFT (FFT is just the fast algorithm of DFT) is the sound equalizer. Sound comes in. DFT converts it into the strength over frequency on the equalizer. You can visually see the strength bars on each frequency slots (from the high, the mid range, to low frequency slots). This let you adjust the frequency strength in each frequency slot individually. You can bump up the high frequency slot or depress the mid range slot. The adjusted frequencies are combined to produce the new distorted sound. When a sound comes in straight, you can only adjust its amplitude (volume). Now with DFT, you can adjust the decomposed frequency components of the sound individually, giving you finer control over the sound's pitch and quality. You can amplify the high frequencies to make it higher pitch. Or amplify the lower frequencies to give more booming bass. You can chop off the low frequency completely to remove background noise. You can chop off both the low and high frequencies to make the sound "clearer." DFT/FFT basically let you control the different dimension of the sound (or any time-based sine wave data).
  - edblarney 10y ago
    That I get. But the FFT results show all sorts of harmonics and results that effectively need to be further decomposed. It does not produce a the results as you described, in reality. Using FFT to decompose someones speech, I'm getting peaks at 300Hz, 600Hz, 900Hz etc.. My crude understanding is that those are all related to one another, though I can't quite fathom why. I fully grasp the example you give ... but the actual results are something bewildering altogether ...
paulsutter 10y ago
Halide is an interesting direction for generating optimized code. By isolating algorithm from schedule, it is better able to optimize for vectorization, locality, pipelining, threading, tiling, and recomputation vs. storage. It then uses LLVM to do all the ordinary optimization. Example code: // The algorithm - no storage or order blur_x(x, y) = (input(x-1, y) + input(x, y) + input(x+1, y))/3; blur_y(x, y) = (blur_x(x, y-1) + blur_x(x, y) + blur_x(x, y+1))/3; // The schedule - defines order, locality; implies storage blur_y.tile(x, y, xi, yi, 256, 32) .vectorize(xi, 8).parallel(y); blur_x.compute_at(blur_y, x).vectorize(x, 8); Halide can auto-optimize to create schedules for specific architectures, which can match or exceed performance of C code written by experts. It's intended for image processing, and has an unnatural affinity for rectangles, but it's an interesting direction for optimization. Darkroom does similar optimization for FPGAs Website: http://halide-lang.org/ http://halide-lang.org/ Video: https://www.youtube.com/watch?v=3uiEyEKji0M https://www.youtube.com/watch?v=3uiEyEKji0M Paper: http://graphics.cs.cmu.edu/projects/halidesched/ http://graphics.cs.cmu.edu/projects/halidesched/
- throwaway_26758 10y ago
  Thanks for bringing up specialized optimization tools such as Halide. Recently there's been a proliferation of domain-specific languages targeting this automatic optimization-from-intentions paradigm. While it's great for testing new ideas isolated from the complexity of modern day compiling suites (academic excursions), actually devoting your codebase to these tools is a huge burden for future maintainability, and eliminates the possibility of future improvements provided outside the small development effort of a single academic group. Many of the techniques given here (and elsewhere in the glut of DSLs) are quietly folded into mainstream C compilers within a short timeframe -- first as requiring explicit hints just like the DSLs, then secondly simplified to allow the compiler to automatically infer hints from well-structured code relying on transformation guarantees of standard C idioms. So if you stick with strict C and come to understand the implicit assumptions/guarantees the standard dictates while conveying structure in your code, in time the compiler suites will march forward in improvements that come with each update 'for free.' And most importantly, update your testing infrastructure to be formally verified each refactoring (such as with CBMC) so that there's no regression problems caused by compiler updates with free reign!
babbage12 10y ago
I love finding blogs like this that I did not know about before. Thanks HN! If anybody knows any other cool blogs in the same vein, please tell me!
deleted 10y ago
[deleted]
nayuki 10y ago
Author here. I am embarrassed and surprised to find this article on the front page of HN. I think it isn't well-written, and I never intended it to appeal to a wide audience. The only fact this article really proves is that I can write an x86 AVX implementation of FFT, such that it has reasonable performance and reasonably readable code. I threw together the C code quickly - the portable version is basically the result of "what if I implement FFT in a relatively naive way?", and the AVX model version is the result of "what if I simulate the computation order and memory access pattern of the AVX code in pure C?". You wouldn't be wrong to say that I created strawmen C code as comparisons to the AVX code. I view the article as a non-thorough wrapper around a proof-of-concept AVX code. I failed on a lot of fronts - comparing performance against well-known high-performance libraries like FFTW, explaining why I designed my AVX code the way it is, elaborating on the memory access patterns, et cetera. For these I apologize to the reader who had their hopes up.
- nkurz 10y ago
  No need to be embarrassed. I think it's useful article, both as a pedagogical "proof of concept" and as a way of evaluating the ability of current compilers to vectorize code. I used a shell script to compile your code with a variety of compilers (gcc-4.8, gcc-4.9, gcc-5, gcc-6, icc-16, icc-17, and clang-4.0) each with a variety of compilation options (-O0, -O1, -O2, -O3, -Ofast, -march=native, -flto/-ipo) to see how they fared on a recent Intel Skylake machine. Across all combinations, I got results that mostly matched what you found. Some takeaways (highly influenced by my own opinions): Avoid -O0 for production code. It no longer means "without optimization", and should be understood as "put everything on the stack for easy debugging". The difference between -O1, -O2, -O3, and -Ofast can be large, but in this case the spread was small (about 10%). Higher numbers usually mean faster performance, but not always. If it really matters, measure. If you have to guess, use -Ofast (which for GCC is approximately -O3 plus -ffast-math). The flag that makes the biggest positive difference here is "-march=native", which tells the compiler it can use all the instructions supported by the machine you are compiling on. If you are running the software yourself rather than distributing it to others, this should probably be your default. If you don't specify it (or another architecture baseline) your code will run on much older machines, but at the cost of missing many of the vectorization improvements made over the last decade. Link time optimization (-flto for gcc and clang) or interprocedural optimization (-ipo in icc) lets the compiler optimize across compilation units rather than only within them. Approximately, it makes everything a candidate for inlining, thus giving many more opportunities for optimization. It doesn't make a big difference here (since almost all the work is done in one function) but you should likely be using it more than you are. Overall, the "portable" version is slower than the "model", which in turn is slower than the "assembly". One read of this is that we don't yet have a "sufficiently clever compiler": the "naive" C version is slower than the "explicitly vectorizable" one, and both are beaten by "naive assembly". This is definitely interesting, and often true, and counters the refrain of "don't bother trying to beat the compiler", but I think it's worth digging deeper and asking why the simplistic assembly is faster than the complicated compiler. I think the main reason here is that the C code is not sufficiently precise, and doesn't actually reflect the programmer's intent. The compiler fails to vectorize because the code as written isn't safe to vectorize: legally, there could be "aliasing" between the different arrays. The assembly is written to take advantage of the absence of aliasing, but the compiler wasn't given enough information to do this. Some of this information can be provided by a liberal sprinkling of "restrict" keywords: "void fft_transform(const void * tables, double * restrict real, double * restrict imag)". But I don't think it's fair to put the blame on the programmer for skipping these. Rather than being safe and assuming that the arrays for "real" and "imag" might overlap, in this case it would clearly be more helpful for the compiler to ask the programmer for more guidance rather than silently giving up. Intel's been doing a fairly good job of giving icc (and icpc) the ability to explain why loops cannot be vectorized. "-qopt-report" and "-guide-vec" give useful hints about what to be changed to allow vectorization. But it still requires significant expertise to figure out what the hints mean, which pragmas to add, and whether it's safe to add them. In this case, the hand written assembly is faster not because the compiler is generating bad code, but because it incorporates the programmer's silent assumptions. It would be nice if compilers were able to give better feedback to the programmer, as a means of making these assumptions (or their absence) more explicit. Do gcc and clang have options to provide this guidance? Are there outside tools that can? compilation script: https://gist.github.com/nkurz/ff4f6a495428c5c0eaffb3521ad8439b https://gist.github.com/nkurz/ff4f6a495428c5c0eaffb3521ad843... portable results: https://gist.github.com/nkurz/a44d57538d4530362951beb1ea2c4cad https://gist.github.com/nkurz/a44d57538d4530362951beb1ea2c4c... model results: https://gist.github.com/nkurz/2f362185b7fd42889b8bb533ced5f7d5 https://gist.github.com/nkurz/2f362185b7fd42889b8bb533ced5f7... assembly results: https://gist.github.com/nkurz/6e797e8b635931d2a0e4f32af6a12d1a https://gist.github.com/nkurz/6e797e8b635931d2a0e4f32af6a12d...
Too 10y ago
The real and imag-pointers should be marked with "restrict"-keyword to enable aliasing optimizations. This goes for all other pointers also. And as others have already said, leaving optimzation at -O1 and not targeting your architecture (-march) can also have big impact. real[j + halfsize] = real[j] - tpre; imag[j + halfsize] = imag[j] - tpim; // 1. Reading imag[j] real[j] += tpre; // 2. Writing real[j], which could be same data as imag[j] imag[j] += tpim; // 3. Must re-read imag[j] even though we read it in #1, because real[j] might alias it
dharma1 10y ago
This is a pretty nice CUDA implementation of sFFT http://hgpu.org/?p=16805 http://hgpu.org/?p=16805

158 points by gozzoo 10y ago

> gcc -O1 -o fft-test fft-test.c fft-portable.c -lm -O1 instructs the compiler to leave some optimizations on the table. If you want the fastest code the compiler can generate, you want -O2 or -O3. This is true for clang and gcc.

user5994461 10y ago
And I believe he should also give flags to target the specific hardware and instruction set that he has. In extreme cases, GCC/Intel/MSVC are able to optimize pure C code with auto-vectorization and SSE/AVX instructions, possibly better than hand crafted ASM... but the compiler will never do that if not given the flags to allow it.
- astrange 10y ago
  Autovectorization is okay but not that good. Its advantage is just that it can transform any part of the program, when your limited dev time won't touch it all. But it can't focus on the hot loops or know about cache efficiency etc. even if you use PGO. And it's sensitive and unstable wrt code changes which is never good.
  - aarongolliver 10y ago
    Yeah, and when autovectorization fails it can be really hard to tell why. The diagnostic messages aren't the best (or, at least, they weren't when I was doing this ~4 years ago) EG: ICC has a AVX version of atan2(), which you can use through intrinsics[0]. But it refused to vectorize even the simplest loop with the function call in it. Even with every possible fast-math flag or every #pragma(please vectorize), it wouldn't do it. I assume due to reduced precision? Nothing I could find online had the answer (again, 4 years ago) [0] https://software.intel.com/en-us/node/524364 https://software.intel.com/en-us/node/524364
benchaney 10y ago
To be honest, this oversight makes the rest of his analysis somewhat suspect.
- dagss 10y ago
  I would assume it meant that it made things slower, after all this was a hand-tuned numerical core. Optimizations is for what a human did not already optimize.
  - benchaney 10y ago
    If by this you mean higher levels of optimization made things slower, than I must disagree. He says "Compiling C code with optimizations enabled is obviously a good idea" he would have mentioned it here if -O2 or -O3 were slower. Also, the optimizer doesn't optimize hand written assembly. Only c code. What I actually think happened is he didn't understand that there were higher optimization levels that -O1, which gives his code an advantage when compared to the c code implementation.
    - dagss 10y ago
      He compared with baseline C code with -O1? That is really bad.
      - benchaney 10y ago
        Well, it is hard to actually tell but that is what it looks like.

llukas 10y ago

Any comparison to fftw or mkl?

nitrogen 10y ago
Also kissfft numbers would be interesting for comparison. Even though it's slower, it has a very easy license to work with.
etrain 10y ago
I would also be interested in seeing this.

kken 10y ago

"The hand-coded x86-64 AVX code is from 1.0× to 2.4× faster than the best C code (which is the AVX model code). This speedup was not as impressive as hoped." Curious. So assembly optimization does not really help that much with modern CPUs?

iammyIP 10y ago
the usual speedup with sse/avx is in about that 2.5 range from what i've observed. If it's lower, then there are likely other effects that stall the execution, like cache misses and other data logistic ineffieciency. The fft shown here uses a precalced lookup table for sin & cos. For a fft size of 2048 this takes 2 * 8 * 2048 bytes, around 32 kb, which is the already size of the usual L1 cache. The bigger the fft size the bigger the lookup table and the less all the asm will be effective. The asm benefits more of a pure algorithmic approach, might try some sin cos approximation instead of a lookup table, then the speed gain might be more consistent.
user5994461 10y ago
Whatever can be done with assembly can be done the same with C (or C++). So strictly speaking, assembly optimization can never help you [compared to not using assembly].
- benchaney 10y ago
  This isn't true at all, in fact, the opposite is true. Compilers are often better than naive micro optimizations, but they do not always have every optimization programmed for every architecture. An expert can absolutely beat the compiler in some cases.
  - Veedrac 10y ago
    Although I agree the compiler rarely produces great code automatically, it's normally possible to handhold the compiler into producing code roughly equivalent to hand-written ASM. Doing so is rarely easier, but it is at least more portable!

mathetic 10y ago

The page doesn't mention compiling code with -O2 or -O3. Does anyone know the reason for this? In particular -O3 is meant to enable vectorisation optimisations, which seems like it might subsume the author's hand optimisations. Also I believe Intel C Compiler is better at optimising numerical code, so results with that would also be illuminating.

raverbashing 10y ago
You can enable vectorization manually without going for -O2 or -O3
- loeg 10y ago
  It isn't clear the author did so, though.
  - mathetic 10y ago
    Actually, it is clear the author didn't because he gives the shell commands that are used to compile.

andrepd 10y ago

A few questions: Is this just vectorization? Because FFT libraries often implement this using either asm or intrinsics. Why compile with -O1? Why not let the gcc optimize with -O2 or -O3? In recent versions I believe this enables vectorization attempts which could produce similar speedups to your handwritten asm.

mpreda 10y ago

The GIMPS (Mersenne numbers search) project has a really extremely fast AVX2 assembler implementation of double FFT http://www.mersenne.org/download/ http://www.mersenne.org/download/ .

adito 10y ago
So, what it's licensed under? Can I use it in other project? From quick glance, it seems that they have their own license.
- mwfunk 10y ago
  The terms and conditions are right here, you can evaluate it for yourself (you are the only person who knows if it suits your needs or not): http://www.mersenne.org/legal/#TCU http://www.mersenne.org/legal/#TCU

rjtobin 10y ago

Here are some performances with higher optimization parameters and with the Intel C Compiler. All run with a Haswell Xeon (E5-2698 v3). Using icc, for the larger problem sizes you get performance speedup that very similar to the speedup from the hard-coded asm. With same arguments as blog post: gcc -O1 -o fft-test fft-test.c fft-portable.c -lm: ... 262144 min=15721397 mean=15748901 sd=0.12% 524288 min=34474976 mean=34546878 sd=0.18% 1048576 min=99814826 mean=100391926 sd=0.26% ... Adding -march=native and -O2: ... 262144 min=14668178 mean=14701671 sd=0.14% 524288 min=31297081 mean=31395468 sd=0.14% 1048576 min=94027867 mean=94196746 sd=0.13% ... With Intel C Compiler (-O2 -xHost): ... 262144 min=14258167 mean=14297606 sd=0.16% 524288 min=31126755 mean=31212153 sd=0.23% 1048576 min=92893600 mean=93491419 sd=0.46% ... Can provide more numbers if people want them. EDIT: looks like I was using the less optimized C code, will re-run with the better performing C code when I get the chance.

inlineint 10y ago
Did you try GCC with -O3 and -ffast-math? It seems like ICC has its analogue to ffast-math enabled by default: http://scicomp.stackexchange.com/questions/20815/performance-of-icc-main-cpp-g-ffast-math-main-cpp http://scicomp.stackexchange.com/questions/20815/performance...
- conistonwater 10y ago
  Good point. GCC has a whole bunch of options here: https://gcc.gnu.org/wiki/FloatingPointMath https://gcc.gnu.org/wiki/FloatingPointMath, and icc's default setting is some, possibly unknown, combination of them.
simcop2387 10y ago
I just started a run with the better C code from a sibling commentor. I've also got it doing -ffast-math, and each of the -O1,2,3. I'll post the raw data up to https://gist.github.com/simcop2387/32c8fb863403ee714227afe0ee3dd637 https://gist.github.com/simcop2387/32c8fb863403ee714227afe0e...

tsomctl 10y ago

How does it compare to the two giants, FFTW and Intel IPP? (http://www.fftw.org/speed/CoreDuo-3.0GHz-icc64/ http://www.fftw.org/speed/CoreDuo-3.0GHz-icc64/). In particular, I don't see any mention of cache awareness. I don't know much about FFTW, but all the parts of Intel IPP that I've used are ridiculously fast, and I'm sure Intel is using all the tricks they can to speed up the cache (ie proper layout, cache prefetch). For example, Intel's fft forward transform on complex floats is about as fast as memcpy with a cold cache.

dsharlet 10y ago
It's not going to compare well. The assembly is a pretty simple port of a basic C implementation of radix 2 Cooley-Tukey FFTs. The author gets ~2x speedup from writing this in assembly. In contrast, highly optimized FFT implementations that use higher radix transforms and other FFT algorithm optimizations, along with SIMD and good C (not necessarily assembly) will get speedups of 10-40x over basic C implementations with a radix 2 implementation: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.153.6089&rep=rep1&type=pdf http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.153... (Implementing FFTs in practice by Frigo and Johnson, of FFTW) FWIW, it is possible to significantly beat FFTW, if you highly optimize for particular cases: https://github.com/halide/Halide/pull/977 https://github.com/halide/Halide/pull/977. I don't have data handy, but I've seen similar results comparing against IPP.
gsteinb88 10y ago
FFTW is one of the seminal examples of the importance of cache awareness, if I remember my history correctly. As I understand it, the reason the fftw authors (Matteo Frigo and Steven G. Johnson) won the Wilkinson prize was essentially because fftw exemplified the importance of cache-aware algorithms, as opposed to blindly optimizing number of operations. I don't know, but I would guess that intel's implementation post-dates fftw?
alfalfasprout 10y ago
Probably a massive difference. The generic X86 IPP implementation of 1 dimensional FFT is around 8x faster than the naive FFT variants I've tried. Intel also claims an effective 5.5x improvement over that when targeting AVX2 (so Haswell and later). Essentially, a 40x is not unusual. When it comes to insanely optimized code, most of Intel's libraries (MKL, VSL, IPP, DAAL) are hard to beat.

ziikutv 10y ago

Great work, just finished a course in DSP too!

throwaway_26758 10y ago

I always cringe when developers compare compiler-optimized C to 'hand rolled' assembly. In the dawning age of SAT-based program synthesis and vectorization-aware superoptimizers, it's extremely anachronistic. It only shows a lack of awareness that in reality your job as a programmer is not to give instructions to a machine, but to provide enough help that a machine can create its own instructions better than you can. For example, let's take one critical function here: bit reversal. The C code provided (renamed variables so it's easy to consistently compare) has an impressive speedup just going from -O1 to higher optimization levels. And then, on top of that speedup with a decent optimizer, we can again halve the amount of time spent in this function by providing better intentions to the compiler which allow it to perform automatic vectorization. Note this is all in standard C, and makes the code simpler to understand (and maintain) by humans as well. uint64_t reverse_slow(uint64_t n, uint64_t k) { uint64_t r, i; for (r = 0, i = 0; i < k; ++i, n >>= 1) r = (r << 1) | (n & 1); return r; } uint64_t reverse_fast(const uint64_t n, const uint64_t k) { uint64_t r, i; for (r = 0, i = 0; i < k; ++i) r |= ((n >> i) & 1) << (k - i - 1); return r; } Lesson: do your job as a human and provide guidance, let the machines do their job at taking your advice and running wild with creativity. It helps you both become more effective.

conistonwater 10y ago
While what you write has a point, I think it's not quite right. I remember stumbling on this HN conversation recently: https://news.ycombinator.com/item?id=9202858 https://news.ycombinator.com/item?id=9202858 — the two realistic examples mentioned there are video codecs and language interpreters. > In the dawning age of SAT-based program synthesis and vectorization-aware superoptimizers, it's extremely anachronistic. I find this particular argument unconvincing, because those aren't real yet. (Consider also that if it was easy and feasible, they probably would have been real already.) In particular, with respect to automatic vectorization (not even superoptimizers!), there is one curious thing that should be noted: clang, with its automatic vectorizer, gives you special options to check if the vectorizer is working for your code (http://llvm.org/docs/Vectorizers.html#diagnostics http://llvm.org/docs/Vectorizers.html#diagnostics). This is necessary: from the point of view of the programmer, performance-sensitive code that is not auto-vectorized is miscompiled because vectorizing it and the consequent performance gain is part of its correctness. One group of examples of this is all the linear algebra libraries. > do your job as a human and provide guidance This advice is only good in the general setting, where it is indeed good advice. The people whose job it is to write all that performance-sensitive code absolutely must get it right, and saying "a compiler will eventually do it for you" is simply not good enough. FWIW, I think you're not addressing your argument at the right crowd.

edblarney 10y ago

Performance? Man, I wish someone would release a decent tutorial on how to actually use the FFT in real life. There are some decent explanations of FT, but often they fail at the very end of FFT: what the output means, explicitly. Particularly, this concern of 'the same freq. appearing in many buckets' due to issues of the phase of the wave? So much theory, so much Matlab ... so little actual applied Engineering. Making 'stuff work' in the real world is basket of caveats ... looking forward to that article: 'FFT in the wild' :)

ww520 10y ago
I was baffled by FFT, too, when I first encountered in math class. It seemed useless. It's not until I started working with sound that it really dawned on me how it's being used. I think the best illustrated usage of DFT (FFT is just the fast algorithm of DFT) is the sound equalizer. Sound comes in. DFT converts it into the strength over frequency on the equalizer. You can visually see the strength bars on each frequency slots (from the high, the mid range, to low frequency slots). This let you adjust the frequency strength in each frequency slot individually. You can bump up the high frequency slot or depress the mid range slot. The adjusted frequencies are combined to produce the new distorted sound. When a sound comes in straight, you can only adjust its amplitude (volume). Now with DFT, you can adjust the decomposed frequency components of the sound individually, giving you finer control over the sound's pitch and quality. You can amplify the high frequencies to make it higher pitch. Or amplify the lower frequencies to give more booming bass. You can chop off the low frequency completely to remove background noise. You can chop off both the low and high frequencies to make the sound "clearer." DFT/FFT basically let you control the different dimension of the sound (or any time-based sine wave data).
- edblarney 10y ago
  That I get. But the FFT results show all sorts of harmonics and results that effectively need to be further decomposed. It does not produce a the results as you described, in reality. Using FFT to decompose someones speech, I'm getting peaks at 300Hz, 600Hz, 900Hz etc.. My crude understanding is that those are all related to one another, though I can't quite fathom why. I fully grasp the example you give ... but the actual results are something bewildering altogether ...

paulsutter 10y ago

Halide is an interesting direction for generating optimized code. By isolating algorithm from schedule, it is better able to optimize for vectorization, locality, pipelining, threading, tiling, and recomputation vs. storage. It then uses LLVM to do all the ordinary optimization. Example code: // The algorithm - no storage or order blur_x(x, y) = (input(x-1, y) + input(x, y) + input(x+1, y))/3; blur_y(x, y) = (blur_x(x, y-1) + blur_x(x, y) + blur_x(x, y+1))/3; // The schedule - defines order, locality; implies storage blur_y.tile(x, y, xi, yi, 256, 32) .vectorize(xi, 8).parallel(y); blur_x.compute_at(blur_y, x).vectorize(x, 8); Halide can auto-optimize to create schedules for specific architectures, which can match or exceed performance of C code written by experts. It's intended for image processing, and has an unnatural affinity for rectangles, but it's an interesting direction for optimization. Darkroom does similar optimization for FPGAs Website: http://halide-lang.org/ http://halide-lang.org/ Video: https://www.youtube.com/watch?v=3uiEyEKji0M https://www.youtube.com/watch?v=3uiEyEKji0M Paper: http://graphics.cs.cmu.edu/projects/halidesched/ http://graphics.cs.cmu.edu/projects/halidesched/

throwaway_26758 10y ago
Thanks for bringing up specialized optimization tools such as Halide. Recently there's been a proliferation of domain-specific languages targeting this automatic optimization-from-intentions paradigm. While it's great for testing new ideas isolated from the complexity of modern day compiling suites (academic excursions), actually devoting your codebase to these tools is a huge burden for future maintainability, and eliminates the possibility of future improvements provided outside the small development effort of a single academic group. Many of the techniques given here (and elsewhere in the glut of DSLs) are quietly folded into mainstream C compilers within a short timeframe -- first as requiring explicit hints just like the DSLs, then secondly simplified to allow the compiler to automatically infer hints from well-structured code relying on transformation guarantees of standard C idioms. So if you stick with strict C and come to understand the implicit assumptions/guarantees the standard dictates while conveying structure in your code, in time the compiler suites will march forward in improvements that come with each update 'for free.' And most importantly, update your testing infrastructure to be formally verified each refactoring (such as with CBMC) so that there's no regression problems caused by compiler updates with free reign!

babbage12 10y ago

I love finding blogs like this that I did not know about before. Thanks HN! If anybody knows any other cool blogs in the same vein, please tell me!

deleted 10y ago

[deleted]

nayuki 10y ago

Author here. I am embarrassed and surprised to find this article on the front page of HN. I think it isn't well-written, and I never intended it to appeal to a wide audience. The only fact this article really proves is that I can write an x86 AVX implementation of FFT, such that it has reasonable performance and reasonably readable code. I threw together the C code quickly - the portable version is basically the result of "what if I implement FFT in a relatively naive way?", and the AVX model version is the result of "what if I simulate the computation order and memory access pattern of the AVX code in pure C?". You wouldn't be wrong to say that I created strawmen C code as comparisons to the AVX code. I view the article as a non-thorough wrapper around a proof-of-concept AVX code. I failed on a lot of fronts - comparing performance against well-known high-performance libraries like FFTW, explaining why I designed my AVX code the way it is, elaborating on the memory access patterns, et cetera. For these I apologize to the reader who had their hopes up.

nkurz 10y ago
No need to be embarrassed. I think it's useful article, both as a pedagogical "proof of concept" and as a way of evaluating the ability of current compilers to vectorize code. I used a shell script to compile your code with a variety of compilers (gcc-4.8, gcc-4.9, gcc-5, gcc-6, icc-16, icc-17, and clang-4.0) each with a variety of compilation options (-O0, -O1, -O2, -O3, -Ofast, -march=native, -flto/-ipo) to see how they fared on a recent Intel Skylake machine. Across all combinations, I got results that mostly matched what you found. Some takeaways (highly influenced by my own opinions): Avoid -O0 for production code. It no longer means "without optimization", and should be understood as "put everything on the stack for easy debugging". The difference between -O1, -O2, -O3, and -Ofast can be large, but in this case the spread was small (about 10%). Higher numbers usually mean faster performance, but not always. If it really matters, measure. If you have to guess, use -Ofast (which for GCC is approximately -O3 plus -ffast-math). The flag that makes the biggest positive difference here is "-march=native", which tells the compiler it can use all the instructions supported by the machine you are compiling on. If you are running the software yourself rather than distributing it to others, this should probably be your default. If you don't specify it (or another architecture baseline) your code will run on much older machines, but at the cost of missing many of the vectorization improvements made over the last decade. Link time optimization (-flto for gcc and clang) or interprocedural optimization (-ipo in icc) lets the compiler optimize across compilation units rather than only within them. Approximately, it makes everything a candidate for inlining, thus giving many more opportunities for optimization. It doesn't make a big difference here (since almost all the work is done in one function) but you should likely be using it more than you are. Overall, the "portable" version is slower than the "model", which in turn is slower than the "assembly". One read of this is that we don't yet have a "sufficiently clever compiler": the "naive" C version is slower than the "explicitly vectorizable" one, and both are beaten by "naive assembly". This is definitely interesting, and often true, and counters the refrain of "don't bother trying to beat the compiler", but I think it's worth digging deeper and asking why the simplistic assembly is faster than the complicated compiler. I think the main reason here is that the C code is not sufficiently precise, and doesn't actually reflect the programmer's intent. The compiler fails to vectorize because the code as written isn't safe to vectorize: legally, there could be "aliasing" between the different arrays. The assembly is written to take advantage of the absence of aliasing, but the compiler wasn't given enough information to do this. Some of this information can be provided by a liberal sprinkling of "restrict" keywords: "void fft_transform(const void * tables, double * restrict real, double * restrict imag)". But I don't think it's fair to put the blame on the programmer for skipping these. Rather than being safe and assuming that the arrays for "real" and "imag" might overlap, in this case it would clearly be more helpful for the compiler to ask the programmer for more guidance rather than silently giving up. Intel's been doing a fairly good job of giving icc (and icpc) the ability to explain why loops cannot be vectorized. "-qopt-report" and "-guide-vec" give useful hints about what to be changed to allow vectorization. But it still requires significant expertise to figure out what the hints mean, which pragmas to add, and whether it's safe to add them. In this case, the hand written assembly is faster not because the compiler is generating bad code, but because it incorporates the programmer's silent assumptions. It would be nice if compilers were able to give better feedback to the programmer, as a means of making these assumptions (or their absence) more explicit. Do gcc and clang have options to provide this guidance? Are there outside tools that can? compilation script: https://gist.github.com/nkurz/ff4f6a495428c5c0eaffb3521ad8439b https://gist.github.com/nkurz/ff4f6a495428c5c0eaffb3521ad843... portable results: https://gist.github.com/nkurz/a44d57538d4530362951beb1ea2c4cad https://gist.github.com/nkurz/a44d57538d4530362951beb1ea2c4c... model results: https://gist.github.com/nkurz/2f362185b7fd42889b8bb533ced5f7d5 https://gist.github.com/nkurz/2f362185b7fd42889b8bb533ced5f7... assembly results: https://gist.github.com/nkurz/6e797e8b635931d2a0e4f32af6a12d1a https://gist.github.com/nkurz/6e797e8b635931d2a0e4f32af6a12d...

Too 10y ago

The real and imag-pointers should be marked with "restrict"-keyword to enable aliasing optimizations. This goes for all other pointers also. And as others have already said, leaving optimzation at -O1 and not targeting your architecture (-march) can also have big impact. real[j + halfsize] = real[j] - tpre; imag[j + halfsize] = imag[j] - tpim; // 1. Reading imag[j] real[j] += tpre; // 2. Writing real[j], which could be same data as imag[j] imag[j] += tpim; // 3. Must re-read imag[j] even though we read it in #1, because real[j] might alias it

dharma1 10y ago

This is a pretty nice CUDA implementation of sFFT http://hgpu.org/?p=16805 http://hgpu.org/?p=16805