Как перенос вычислений из времени выполнения во время компиляции замедлил мой код?
Я пишу функцию для вычисления быстрого преобразования Фурье с использованием алгоритма Стокхэма БПФ и обнаружил, что все "факторы поворота" для вычисления могут быть предварительно рассчитаны во время компиляции, если длина БПФ равна степени двойки.
В вычислениях БПФ вычисления коэффициента вращения обычно занимают большую часть общего времени, поэтому теоретически выполнение этого должно привести к резкому повышению производительности алгоритма.
Вчера я потратил весь день на повторную реализацию своего алгоритма на новом компиляторе (gcc 10), чтобы я мог использовать C++20. consteval
функция для предварительного расчета всех факторов вращения во время компиляции. Мне удалось это сделать, но в итоге код, который предварительно вычисляет все факторы поворота во время компиляции, действительно работает медленнее!
Вот код, который выполняет все вычисления во время выполнения:
#include <algorithm>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <iostream>
#include <vector>
using namespace std;
static vector<complex<double>> StockhamFFT(const vector<complex<double>>& x);
constexpr bool IsPowerOf2(const size_t value)
{
return value && (!(value & (value - 1)));
}
vector<complex<double>> FFT(const vector<double>& x)
{
const auto N = x.size();
assert(IsPowerOf2(x.size()));
const auto NOver2 = N/2;
vector<complex<double>> x_p(N);
transform(x.begin(), x.end(), x_p.begin(), [](const double value) {
return complex<double>(value);
});
return StockhamFFT(x_p);
}
// C++ implementation of the Stockam FFT algorithm
static vector<complex<double>> StockhamFFT(const vector<complex<double>>& x)
{
const auto N = x.size();
assert(IsPowerOf2(N));
const auto NOver2 = N/2;
// Pre-calculate the twiddle factors (at runtime)
vector<complex<double>> W(NOver2);
const auto omega = 2.0 * M_PI / N;
for (size_t n = 0; n < NOver2; ++n)
{
W[n] = complex{ cos(-omega * n), sin(-omega * n) };
}
// The Stockham algorithm requires one vector for input/output data and
// another as a temporary workspace
vector<complex<double>> a(x);
vector<complex<double>> b(N);
// Set the spacing between twiddle factors used at the first stage
size_t WStride = N/2;
// Loop through each stage of the FFT
for (size_t stride = 1; stride < N; stride *= 2)
{
// Loop through the individual FFTs of each stage
for (size_t m = 0; m < NOver2; m += stride)
{
const auto mTimes2 = m*2;
// Perform each individual FFT
for (size_t n = 0; n < stride; ++n)
{
// Calculate the input indexes
const auto aIndex1 = n + m;
const auto aIndex2 = aIndex1 + NOver2;
// Calculate the output indexes
const auto bIndex1 = n + mTimes2;
const auto bIndex2 = bIndex1 + stride;
// Perform the FFT
const auto tmp1 = a[aIndex1];
const auto tmp2 = W[n*WStride]*a[aIndex2];
// Sum the results
b[bIndex1] = tmp1 + tmp2;
b[bIndex2] = tmp1 - tmp2; // (>*.*)> symmetry! <(*.*<)
}
}
// Spacing between twiddle factors is half for the next stage
WStride /= 2;
// Swap the data (output of this stage is input of the next)
a.swap(b);
}
return a;
}
int main()
{
size_t N = pow(2, 18);
vector<double> x(N);
int f_s = 1000;
double t_s = 1.0 / f_s;
for (size_t n = 0; n < N; ++n)
{
x[n] = sin(2 * M_PI * 100 * n * t_s);
}
auto sum = 0;
for (int i = 1; i < 100; ++i)
{
auto start = chrono::high_resolution_clock::now();
auto X = FFT(x);
auto stop = chrono::high_resolution_clock::now();
auto duration = chrono::duration_cast<chrono::microseconds>(stop - start);
sum += duration.count();
}
auto average = sum / 100;
std::cout << "duration " << average << " microseconds." << std::endl;
}
Используя это в качестве отправной точки, я смог извлечь расчеты коэффициента вращения из StockhamFFT
функции и выполнять их во время компиляции с помощью consteval
функция. Вот как выглядит код после:
#include <algorithm>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <iostream>
#include <vector>
using namespace std;
static vector<complex<double>> StockhamFFT(const vector<complex<double>>& x);
constexpr bool IsPowerOf2(const size_t value)
{
return value && (!(value & (value - 1)));
}
// Evaluates FFT twiddle factors at compile time!
template <size_t N>
static consteval array<complex<double>, N/2> CalculateTwiddleFactors()
{
static_assert(IsPowerOf2(N), "N must be a power of 2.");
array<complex<double>, N/2> W;
const auto omega = 2.0*M_PI/N;
for (size_t n = 0; n < N/2; ++n)
{
W[n] = complex{cos(-omega*n), sin(-omega*n)};
}
return W;
}
// Calculate the twiddle factors (>*O*)> AT COMPILE TIME <(*O*<)
constexpr auto W = CalculateTwiddleFactors<static_cast<size_t>(pow(2,18))>();
vector<complex<double>> FFT(const vector<double>& x)
{
const auto N = x.size();
assert(IsPowerOf2(x.size()));
const auto NOver2 = N/2;
vector<complex<double>> x_p(N);
transform(x.begin(), x.end(), x_p.begin(), [](const double value) {
return complex<double>(value);
});
return StockhamFFT(x_p);
}
// C++ implementation of the Stockam FFT algorithm
static vector<complex<double>> StockhamFFT(const vector<complex<double>>& x)
{
const auto N = x.size();
assert(IsPowerOf2(N));
const auto NOver2 = N/2;
//***********************************************************************
// Twiddle factors are already calculated at compile time!!!
//***********************************************************************
// The Stockham algorithm requires one vector for input/output data and
// another as a temporary workspace
vector<complex<double>> a(x);
vector<complex<double>> b(N);
// Set the spacing between twiddle factors used at the first stage
size_t WStride = N/2;
// Loop through each stage of the FFT
for (size_t stride = 1; stride < N; stride *= 2)
{
// Loop through the individual FFTs of each stage
for (size_t m = 0; m < NOver2; m += stride)
{
const auto mTimes2 = m*2;
// Perform each individual FFT
for (size_t n = 0; n < stride; ++n)
{
// Calculate the input indexes
const auto aIndex1 = n + m;
const auto aIndex2 = aIndex1 + NOver2;
// Calculate the output indexes
const auto bIndex1 = n + mTimes2;
const auto bIndex2 = bIndex1 + stride;
// Perform the FFT
const auto tmp1 = a[aIndex1];
const auto tmp2 = W[n*WStride]*a[aIndex2];
// Sum the results
b[bIndex1] = tmp1 + tmp2;
b[bIndex2] = tmp1 - tmp2; // (>*.*)> symmetry! <(*.*<)
}
}
// Spacing between twiddle factors is half for the next stage
WStride /= 2;
// Swap the data (output of this stage is input of the next)
a.swap(b);
}
return a;
}
int main()
{
size_t N = pow(2, 18);
vector<double> x(N);
int f_s = 1000;
double t_s = 1.0 / f_s;
for (size_t n = 0; n < N; ++n)
{
x[n] = sin(2 * M_PI * 100 * n * t_s);
}
auto sum = 0;
for (int i = 1; i < 100; ++i)
{
auto start = chrono::high_resolution_clock::now();
auto X = FFT(x);
auto stop = chrono::high_resolution_clock::now();
auto duration = chrono::duration_cast<chrono::microseconds>(stop - start);
sum += duration.count();
}
auto average = sum / 100;
std::cout << "duration " << average << " microseconds." << std::endl;
}
Обе версии были скомпилированы на Ubuntu 19.10 с gcc 10.0.1:
g++ -std=c++2a -o main main.cpp
Обратите внимание, что компилятор gcc требуется специально, потому что это единственный компилятор, поддерживающий constexpr
версии sin
а также cos
Пример "времени выполнения" дает следующее:
длительность 292854 микросекунды.
Пример "времени компиляции" дает следующее:
длительность 295230 микросекунд.
Компиляция версии во время компиляции заняла заметно больше времени, но все равно почему-то дольше запускается, хотя большая часть вычислений выполняется еще до запуска программы! Как это возможно? Я упустил что-то важное?