C++ Fast Fourier transform This is a very simple FFT, I am wondering what I can do to make this faster and more memory efficient from the programming side (better data types, and maybe some tricks like unrolling loops or using the preprocessor if that is useful here), and not by using a more efficient mathematical algorithm. Obviously I would appreciate advice on best practices as well.
#include &lt;stdio.h&gt;
#include &lt;vector&gt;
#include &lt;iostream&gt;
#include &lt;complex&gt;
#include &lt;cmath&gt;
#include &lt;algorithm&gt;

#define N 1048576
#define PI 3.14159265358979323846

/*
Creating the table of all N'th roots of unity.
We use notation omega_k = e^(-2 pi i / n).
*/
template&lt;typename U&gt;
std::vector&lt; std::complex&lt;U&gt; &gt; rootsOfUnityCalculator() {
    std::vector&lt; std::complex&lt;U&gt; &gt; table;

    for (size_t k = 0; k &lt; N; k++) {
        std::complex&lt;U&gt; kthRootOfUnity(std::cos(-2.0 * PI * k / N), std::sin(-2.0 * PI * k / N));
        table.push_back(kthRootOfUnity);
    }

    return table;
}

/*
Fast Fourier transform, T is the precision level, so float or double.
table is a look up table of the roots of unity. Overwrites the input.
For now only works for N a power of 2.
*/
template&lt;typename T&gt;
void FFT(std::complex&lt;T&gt;* input, const std::vector&lt; std::complex&lt;T&gt; &gt;&amp; table, size_t n) {

    if (n % 2 == 0) {
        // Split up the input in even and odd components
        std::complex&lt;T&gt;* evenComponents = new std::complex&lt;T&gt;[n/2];
        std::complex&lt;T&gt;* oddComponents = new std::complex&lt;T&gt;[n/2];

        for (size_t k = 0; k &lt; n/2; k++) {
            evenComponents[k] = input[2 * k];
            oddComponents[k] = input[2 * k + 1];
        }

        // Transform the even and odd input
        FFT(evenComponents, table, n/2);
        FFT(oddComponents, table, n/2);

        // Use the algorithm from Danielson and Lanczos
        for (size_t k = 0; k &lt; n/2; k++) {
            std::complex&lt;T&gt; plusMinus = table[N / n * k] * oddComponents[k]; // omega_n^k = (omega_N^(N/n))^k = omega_N^(Nk/n)
            input[k] = evenComponents[k] + plusMinus;
            input[k + n/2] = evenComponents[k] - plusMinus;
        }

        delete[] evenComponents;
        delete[] oddComponents;

    } else {
        // The Fourier transform on one element does not do anything, so
        // nothing needed here.
    }
}

int main() {
    std::complex&lt;double&gt;* input = new std::complex&lt;double&gt;[N];

    for (size_t k = 0; k &lt; N; k++) {
        input[k] = k;
    }

    const std::vector&lt; std::complex&lt;double&gt; &gt; table = rootsOfUnityCalculator&lt;double&gt;();

    // Overwrites the input with its Fourier transform
    FFT&lt;double&gt;(input, table, N);

    delete[] input;

    return 0;
}
Prefer C++ headers (such as &lt;cstdio&gt;) to the C compatibility headers (&lt;stdio.h&gt;).  The C++ headers define identifiers in the std namespace where we want them.
It seems that this header isn't even used here, so we can omit it completely.

std::size_t is misspellt throughout as size_t.  This commonly happens when writing on a platform that declares it in the global namespace as well as in std, but that's not a portable assumption.

We only use PI as a double, so we can declare it as a constant rather than as a macro:
constexpr double pi = 3.14159265358979323846;

Or to make it more obviously correct:
const double pi = std::acos(-1.0);


We can avoid naked new[] and delete[] using std::vector.  If we use std::vector::reserve before adding any elements, there should be no performance impact, only improved code safety.
Using a vector saves having to pass the array size explicitly, too.

A helpful technique to determine whether a number is a power of two is that bitwise-and of the number and its arithmetic predecessor is zero for powers of two:
constexpr bool is_power_of_two(unsigned n)
{
    return n &amp;&amp; (n &amp; (n-1) == 0);
}


A single value of N is probably not very useful in general.  We could make it a template parameter, or a simple function argument, to make our code useful for other sizes of input.

Modified version
I split FFT() into a simple public function that ensures arguments are valid and the internal recursive implementation.  I've used a requires to constrain the complex type we accept (integer types aren't useful here), but that can be omitted for plain C++17.
#include &lt;complex&gt;
#include &lt;vector&gt;

namespace internal
{
    /*
      Creating the table of all N'th roots of unity.
      We use notation ω_i = e^(-2πi/n).
    */
    template&lt;typename U&gt;
    constexpr auto rootsOfUnity(std::size_t n)
    {
        constexpr auto pi = std::acos(U{-1.0});
        std::vector&lt;std::complex&lt;U&gt;&gt; table;
        table.reserve(n);

        for (std::size_t i = 0;  i &lt; n;  ++i) {
            table.emplace_back(std::cos(-2.0 * pi * i / n), std::sin(-2.0 * pi * i / n));
        }

        return table;
    }

    template&lt;typename V&gt;        // V is vector of complex
    void FFT(V&amp; input, const V&amp; table, std::size_t total_size)
    {
        if (input.size() &lt;= 1) {
            // The Fourier transform on one element does not do
            // anything, so nothing needed here.
            return;
        }

        const auto n2 = input.size() / 2;

        // Split up the input in even and odd components
        V evenComponents;  evenComponents.reserve(n2);
        V oddComponents;   oddComponents.reserve(n2);

        for (std::size_t i = 0;  i &lt; input.size();  ) {
            evenComponents.push_back(input[i++]);
            oddComponents.push_back(input[i++]);
        }

        // Transform the even and odd input
        FFT(evenComponents, table, total_size);
        FFT(oddComponents, table, total_size);

        // Use the algorithm from Danielson and Lanczos
        for (std::size_t i = 0;  i &lt; n2;  ++i) {
             // ω_n^i = (ω_N^(N/n))^i = ω_N^(Nk/n)
            auto const plusMinus = table[total_size / n2 * i] * oddComponents[i];
            input[i] = evenComponents[i] + plusMinus;
            input[i + n2] = evenComponents[i] - plusMinus;
        }
    }
}

/*
  Fast Fourier transform
  T is the precision level: float, double, long double.
  Overwrites the input.
  Only works for input sizes that are a power of 2.
*/
template&lt;typename T&gt;
    requires std::is_floating_point&lt;T&gt;
void FFT(std::vector&lt;std::complex&lt;T&gt;&gt;&amp; input)
{
    if (input.size() == 0) {
        throw std::invalid_argument(&quot;input is empty&quot;);
    }
    if (input.size() &amp; (input.size() - 1)) {
        throw std::invalid_argument(&quot;input length is not a power of two&quot;);
    }

    auto const table = internal::rootsOfUnity&lt;T&gt;(input.size());
    internal::FFT(input, table, input.size());
}


// Test program
#include &lt;iostream&gt;

int main()
{
    static constexpr size_t n = 1048576;
    std::vector&lt;std::complex&lt;double&gt;&gt; input;
    for (std::size_t i = 0;  i &lt; n;  ++i) {
        input.emplace_back(i);
    }

    // Overwrite input with its Fourier transform
    FFT&lt;double&gt;(input);

    // use the result, to avoid optimising out
    std::cout &lt;&lt; input[0] &lt;&lt; '\n';
}